Table Of Content

6 0 Hopf bifurcation analysis for a 0 2 mathematical model of P53-MDM2 n a interaction J 9 1 Mihaela Neam¸tua∗, Raul Florin Horhatb, Dumitru Opri¸sc ] S D . h t aDepartmentofEconomicInformatics,Mathematics andStatistics,FacultyofEconomics, a m WestUniversityofTimi¸soara,str. Pestalozzi,nr. 16A,300115,Timi¸soara,Romania, E-mail:[email protected], [ b DepartmentofBiophysicsandMedicalInformatics, 1 v UniversityofMedicineandPharmacy,PiataEftimieMurgu,nr. 3,300041, Timi¸soara,Romania, 1 E-mail: [email protected] 8 4 c DepartmentofAppliedMathematics, FacultyofMathematics, 1 WestUniversityofTimi¸soara,Bd. V.Parvan,nr. 4,300223, Timi¸soara,Romania, 0 E-mail: [email protected] 6 0 / h Abstract t a m Inthispaperweanalyzeasimplemathematicalmodelwhichdescribes : v the interaction between the proteins p53 and Mdm2. For the statio- i X nary state we discuss the local stability and the existence of the Hopf r bifurcation. Choosing the delay as a bifurcation parameter we study a the direction and stability of the bifurcating periodic solutions. Some numerical examples and the conclusions are finally made. Keywords: delay differential equation, stability, Hopf bifurcation, p53, Mdm2. 2000AMSMathematicsSubjectClassification:34C23, 34C25, 37G05, 37G15, 92D10 . ∗ Corresponding author 1 1. Introduction In every normal cell, there is a protective mechanism against tumoral degeneration. This mechanism isbased ontheP53 network. P53, alsoknown as ”the guardian of the genome”, is a gene that codes a protein in order to regulate the cell cycle. The name is due to its molecular mass: it is a 53 kilo Dalton fraction of cell proteins. Mdm2 gene plays a very important role in P53network. Itregulatesthelevels ofintracellular P53proteinconcentration through a feedback loop. Under normal conditions the P53 levels are kept very low. When there is DNA damage the levels of P53 protein rise and if there is a prolonged elevation the cell shifts to apoptosis, and if there is only a short elevation the cell cycle is arrested and the repair process is begun. The first pathway protects the cell from tumoral transformation when there is a massive DNA damage that cannot be repaired, and the second pathway protects a number of important cells (neurons, myocardic cells) from death after DNA damage. In these cells first pathway, apoptosis, is not an option because they do not divide in adult life and their importance is obvious. Due toitsmajorimplicationincancer prevention andduetotheactionsdescribed above, P53 has been intensively studied in the last two decades. During the years, several models which describe the interaction between P53andMdm2havebeenstudied. Wementionsomeoftheminreferences[2], [3], [4], [6], [7], [8], [9], [10]. This paper gives a mathematical approach to the model described in [10]. The authors of paper [10] make a molecular energy calculation based on the classical force fields, and they also use chemical reactions constants from literature. Their results obtained by simulations in accordance with experimental behavior of the P53-Mdm2 complex, but they lack the mathematical part which we develop in this paper. We analyze the Hopf bifurcation with time delay as a bifurcation parameter using the methods from [1], [5], [7]. The paper is organized as follows. We present the mathematical model in section 2 and the existence of the stationary state is study. In section 3, we discuss the local stability for the stationary state of system (2) and we investigate the existence of the Hopf bifurcation for system (2) using time delay as the bifurcation parameter. In section 4, the direction of Hopf bifurcation is analyzed according to the normal form theory and the center manifold theorem introduced by Hassard [5]. Numerical simulations for confirming the theoretical results are illustrated in section 5. Finally, some conclusions 2 are made. 2. The mathematical model and the stationary state The state variables are: y (t),y (t) the totalnumber of p53 molecules and 1 2 the total number of Mdm2 proteins. The interaction function between P53 and MDM2 is f : IR2 → IR given + by [10]: 1 f(y ,y ) = (y +y +k − (y +y +k)2 −4y y ). (1) 1 2 1 2 1 2 1 2 2 q The parameters of the model are: s the production rate of p53, a the degradation rate of p53 (through ubiquitin pathway), and also the rate at which Mdm2 re-enters the loop, b the spontaneous decay rate of p53, d the decay rate of the protein rate Mdm2, k the dissociation constant of the 1 complex P53-Mdm2, c the constant of proportionality of the production rate of Mdm2 gene with the probability that the complex P53-Mdm2 is build. These parameters are positive numbers. The mathematical model is described by the following differential system with time delay [10]: y˙ (t) = s−af(y (t),y (t))−by (t), 1 1 2 1 (2) y˙ (t) = cg(y (t−τ),y (t−τ))−dy (t), 2 1 2 2 where f is given by (1) and g : IR2 → IR, is + y −f(y ,y ) 1 1 2 g(y ,y ) = . (3) 1 2 k +y −f(y ,y ) 1 1 1 2 For the study of the model (1) we consider the following initial values: y (θ) = ϕ (θ),y (θ) = ϕ (θ),θ ∈ [−τ,0], 1 1 2 2 with ϕ ,ϕ : [−τ,0] → IR are differentiable functions. 1 2 + In the second equation of (2) there is delay, because the trancription and translation of Mdm2 last for some time after that p53 was bound to the gene. The stationary state (y ,y ) ∈ IR2 is given by the solution of the system 10 20 + of equations: 3 s−af(y ,y )−by = 0, 1 2 1 (4) cg(y ,y )−dy = 0. 1 2 2 From (1), (3) and (4) we deduce that the stationary state can be found through the intersection of the curves: y = f (y ), y = f (y ), (5) 2 1 1 1 2 1 where f ,f : IR → IR are given by: 1 2 + c(y (a+b)−s) 1 f (y ) = , 1 1 d(k a+(a+b)y −s) 1 1 (6) (s−by )(ka+(a+b)y −s) 1 1 f (y ) = . 2 2 a((a+b)y −s) 1 Proposition 1. The function f ,f from (6) have the following proper- 1 2 ties: s s (i) f is strictly increasing, f is strictly decreasing on , ; 1 2 a+b b s s (cid:18) (cid:19) (ii) There is a unique value y ∈ , so that f (y ) = f (y ), 10 1 10 2 10 a+b b (cid:18) (cid:19) s s where y is the solution of the equation ϕ(x) = 0, from , where 10 a+b b (cid:18) (cid:19) ϕ(x)=ac((a+b)x−s)2−d(k a+(a+b)x−s)(ka+(a+b)x−s)(s−bx). (7) 1 Proof. (i) From (6) we have: (a+b)ack f′(y ) = 1 , 1 1 d(k a+(a+b)y −s)2 1 1 (8) b bk k(a+b)(s−by ) f′(y ) = − − − 1 . 2 2 a (a+b)y −s ((a+b)y −s)2 1 1 s s For y ∈ , , the relations (8) lead to f′(y ) > 0, f′(y ) < 0, so that 1 a+b b 1 1 2 1 (cid:18) (cid:19) f is strictly increasing and f is strictly deceasing. 1 2 s s (ii) By (i) there is y ∈ , so that f (y ) = f (y ). From (6) 10 1 10 2 10 a+b b (cid:18) (cid:19) and (7) it results that: ϕ(y ) = f (y )h(y ), 1 12 1 1 4 where h(y ) = ad((a+b)y −s)((a+b)y +k a−s) 1 1 1 1 f (y ) = f (y )−f (y ). 12 1 1 1 2 1 s s On , , the functions f and h are strictly increasing and conse- 12 a+b b (cid:18) (cid:19) s a3dkk s 1 quently ϕ is strictly increasing. Because ϕ( ) = − < 0 and a+b a+b s a3cs2 ϕ( ) = > 0 we can conclude that equation ϕ(x) = 0 has a unique b b2 s s solution y on , . 10 a+b b (cid:18) (cid:19) 3. The analysis of the stationary state and the existence of the Hopf bifurcation. We consider the following translation: y = x +y ,y = x +y 1 1 10 2 2 20 and system (2) can be expressed as: x˙ (t) = s−af(x (t)+y ,x (t)+y )−b(y (t)+y ), 1 1 10 2 20 1 10 (9) x˙ (t) = cg(x (t−τ)+y ,x (t−τ)+y )−d(x (t)+y ). 2 1 10 2 20 2 20 System (9) has a unique stationary state (0,0). To investigate the local stability of the equilibrium state we linearize system (9). Let u (t) and u (t) 1 2 be the linearized system variables. Then (9) is rewritten as: U˙(t) = AU(t)+BU(t−τ), (10) where −(b+aρ ) −aρ 0 0 10 01 A= ,B= (11)  0 −d  cγ cγ  10 01         with U(t) = (u (t),u (t))T and ρ , ρ , γ , γ are the values of the first 1 2 10 01 10 01 order derivatives for the functions: 1 f(x,y) = (x+y +k − (x+y +k)2 −4xy) 2 q 5 and x−f(x,y) g(x,y) = k +x−f(x,y) 1 evaluated at (y ,y ). 10 20 The characteristic equation corresponding to system (10) is ∆(λ,τ) = det(λI −A−e−λτB) = 0 which leads to: λ2 +p λ+p −(q λ+q )e−λτ = 0, (12) 1 0 1 0 where p = b+d+aρ , p = d(b+aρ ), q = cγ , 1 10 0 10 1 01 q = cγ (b+aρ )−acρ γ . 0 01 10 01 10 If there is no delay, the characteristic equation (12) becomes: ∆(λ,0) = λ2 +(p −q )λ+p −q . (13) 1 1 0 0 Then, the stationary state (0,0) is locally asymptotically stable if p −q > 0, p −q > 0. (14) 1 1 0 0 When τ > 0, the stationary state is asymptotically stable if and only if all roots of equation (12) have a negative real part. We are determining the interval [0,τ ) so that the stationary state remain asymptotically stable. 0 InwhatfollowswestudytheexistenceoftheHopfbifurcationforequation (10) choosing τ as the bifurcation parameter. We are looking for the values τ of τ so that the stationary state (0,0) changes from local asymptotic 0 stability to instability or vice versa. We need the pure imaginary solutions of equation (12). Let λ = ±iω be these solutions and without loss of generality 0 we assume ω > 0. Replacing λ = iω and τ = τ in (12) we obtain: 0 0 0 q cosω τ +q ω sinω τ = p −ω2 0 0 0 1 0 0 0 0 0 q sinω τ −q ω cosω τ = −ω p , 0 0 0 1 0 0 0 0 1 which implies that 1 p ω q ω 1 0 1 0 τ = (2k+1)π+arcsin +arcsin 0 ω0 (p −ω2)2+ω2p2 (p −ω2)2+ω2p2 0 0 0 1 0 0 0 1  q q (15) 6 where ω is a solution of the equation: 0 ω4 +(−p2 −2p +q2)ω2 +p2 −q2 = 0. 1 0 1 0 0 dλ Now we have to calculate Re evaluated at λ = iω and τ = τ . We 0 0 dτ ! have: dλ | = M +iN dτ λ=iω0,τ=τ0 where q2ω6 +2q2ω4+(p2q2 −p2q2 −2p q2)ω2 M = 1 0 0 0 1 0 0 1 0 0 0 (16) l2 +l2 1 2 and −q2τ ω7 +ω5(q q −p q2+τ (2p q2 −p2q2 −q2)) N = 1 0 0 0 0 1 1 1 0 0 1 1 1 0 + l2 +l2 1 2 ω3(−p q2−2p q q +p2q q −p p q2+τ (−p2q2−q2p2+2p q2)) + 0 1 0 0 0 1 1 0 1 0 1 1 0 1 0 1 0 0 0 + l2 +l2 1 2 ω (−p p q2 +p2q q −τ p2q2) + 0 0 1 0 0 0 1 0 0 0 l2 +l2 1 2 (17) with l = −q ω2+p q −q p +τ (−q p ω2−q ω2+q p ), 1 1 0 1 0 1 0 0 1 1 0 0 0 0 0 l = 2ω q +τ (−q ω3+p q ω +p q ω ). 2 0 0 0 1 0 0 1 0 1 0 0 We conclude with: Theorem 1. If there is no delay, under condition (14) system (10) has an asymtotically stable stationary state. If τ > 0 and p2q2−q2p2−2p q2 > 0 1 0 1 0 0 0 dλ then there is τ = τ given by (15) so that Re > 0 and therefor 0 dτ ! λ=iω0,τ=τ0 a Hopf bifurcation occurs at (y ,y ). 10 20 4. Direction and stability of the Hopf bifurcation In this section we describe the direction, stability and the period of the bifurcating periodic solutions of system (2). The method we use is based on the normal form theory and the center manifold theorem introduced by Hassard [5]. Taking into account the previous section, if τ = τ then all roots 0 7 of equation (11) other than ±iω have negative real parts, and any root of 0 equation (11) of the form λ(τ) = α(τ)±iω(τ) satisfies α(τ ) = 0, ω(τ ) = ω 0 0 0 dα(τ ) 0 and 6= 0. For notational convenience let τ = τ + µ,µ ∈ IR. Then 0 dτ µ = 0 is the Hopf bifurcation value for equations (2). The Taylor expansion at (y ,y ) of the right members from (9) until the 10 20 third order leads to: x˙(t) = Ax(t)+Bx(t−τ)+F(x(t),x(t−τ)), (18) where x(t) = (x (t),x (t))T, A, B are given by (11) and 1 2 F(x(t),x(t−τ)) = (F1(x(t)),F2(x(t−τ)))T, (19) where a F1(x (t),x (t)) = − [ρ x2(t)+2ρ x (t)x (t)+ρ x2(t)]− 1 2 2 20 1 11 1 2 02 2 a − [ρ x3(t)+3ρ x2(t)x (t)+3ρ x (t)x2(t)+ρ x3(t)] 6 30 1 21 1 2 12 1 2 03 2 c F2(x (t−τ),x (t−τ)) = [γ x2(t−τ)+2γ x (t−τ)x (t−τ)+ 1 2 2 20 1 11 1 2 + γ x2(t−τ)]+ 02 2 c + [γ x3(t−τ)+3γ x2(t−τ)x (t−τ)+ 6 30 1 21 1 2 + 3γ x (t−τ)x2(t−τ)+γ x3(t−τ)] 12 1 2 03 2 and ρ , γ , i,j = 0,1,2,3 are the values of the second and third order ij ij derivatives for the functions f(x,y) and g(x,y). Definethespaceofcontinuousreal-valuedfunctionsasC = C([−τ ,0],IR4). 0 In τ = τ + µ,µ ∈ IR, we regard µ as the bifurcation parameter. For 0 Φ ∈ C we define a linear operator: L(µ)Φ = AΦ(0)+BΦ(−τ) where Aand Bare given by (11) anda nonlinear operatorF(µ,Φ) = F(Φ(0), Φ(−τ)), where F(Φ(0),Φ(−τ)) is given by (19). According to the Riesz representation theorem, there is a matrix whose components are bounded variation functions, η(θ,µ) with θ ∈ [−τ ,0] so that: 0 0 L(µ)Φ = dη(θ,µ)φ(θ), θ ∈ [−τ ,0]. 0 −Zτ0 8 For Φ ∈ C1([−τ ,0],IR4) we define: 0 dΦ(θ) , θ ∈ [−τ ,0) 0 dθ A(µ)Φ(θ) =   0 dη(t,µ)φ(t), θ = 0, −τ0  R 0, θ ∈ [−τ ,0) 0 R(µ)Φ(θ) =  F(µ,Φ), θ = 0.  We can rewrite (18) in the following vector form:  u˙ = A(µ)u +R(µ)u (20) t t t where u = (u ,u )T, u = u(t+θ) for θ ∈ [−τ ,0]. 1 2 t 0 For Ψ ∈ C1([0,τ ],IR∗4), we define the adjunct operator A∗ of A by: 0 dΨ(s) − , s ∈ (0,τ ] 0 ds A∗Ψ(s) =   0 dηT(t,0)ψ(−t), s = 0. −τ0  R We define the following bilinear form: 0 θ < Ψ(θ),Φ(θ) >= Ψ¯T(0)Φ(0)− Ψ¯T(ξ −θ)dη(θ)φ(ξ)dξ, Z−τ0 Zξ=0 where η(θ) = η(θ,0). We assume that ±iω are eigenvalues of A(0). Thus, they are also eigen- 0 ∗ values of A . We can easily obtain: Φ(θ) = veλ1θ, θ ∈ [−τ ,0] (21) 0 where v = (v ,v )T, 1 2 v = −aρ ,v = λ +b+aρ , 1 01 2 1 10 is the eigenvector of A(0) corresponding to λ = iω and 1 0 Ψ(s) = weλ2s, s ∈ [0,∞) 9 where w = (w ,w ), 1 2 w cγ e−λ1τ 1 2 10 w = ,w = , 1 2 b+λ +aρ η¯ 1 10 cγ e−λ2τ cγ v +cγ v η = v 10 +v − 10 1 01 2(−τ λ e−λ1τ0 −1+e−λ1τ0)] 1b+λ +aρ 2 λ2 0 1 2 10 1 is the eigenvector of A(0) corresponding to λ = −iω . 2 0 Wecanverifythat: < Ψ(s),Φ(s) >= 1,< Ψ(s),Φ¯(s) >=< Ψ¯(s),Φ(s) >= 0, < Ψ¯(s),Φ¯(s) >= 1. Using the approach of Hassard [5], we next compute the coordinates to describe the center manifold Ω at µ = 0. Let u = u (t+θ),θ ∈ [−τ ,0) be 0 t t 0 the solution of equation (20) when µ = 0 and z(t) =< Ψ,u >, w(t,θ) = u (θ)−2Re{z(t)Φ(θ)}. t t On the center manifold Ω , we have: 0 w(t,θ) = w(z(t),z¯(t),θ) where z2 z¯2 z3 w(z,z¯,θ) = w (θ) +w (θ)zz¯+w (θ) +w (θ) +... 20 11 02 30 2 2 6 in which z and z¯ are local coordinates for the center manifold Ω in the 0 direction of Ψ and Ψ¯ and w (θ) = w¯ (θ). Note that w and u are real. 02 20 t For solution u ∈ Ω of equation (20), as long as µ = 0, we have: t 0 z˙(t) = λ z(t)+g(z,z¯) (22) 1 where g(z,z¯) = Ψ¯(0)F(w(z(t),z¯(t),0)+2Re(z(t)Φ(0))) = z(t)2 z¯(t)2 z(t)2z¯(t) = g +g z(t)z¯(t)+g +g +... 20 11 02 21 2 2 2 where g = F1 w¯ +F2 w¯ ,g = F1 w¯ +F2 w¯ ,g = F1 w¯ +F2w¯ , (23) 20 20 1 20 2 11 11 1 11 2 02 02 1 02 2 10