Table Of Content6
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 inves-
tigate 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
