Table Of Content6
0 A dynamic p53-mdm2 model with
0
2
distributed time delay
n
a
J
4 M. NEAMT¸Ua∗, D. OPRIS¸b, R.F. HORHATc
2
]
S
D aDepartmentofEconomicInformaticsandStatistics,FacultyofEconomics,
. WestUniversityofTimi¸soara,str. Pestalozzi,nr. 16A,300115,Timi¸soara,Romania,
h
t E-mail:[email protected],
a
m b DepartmentofAppliedMathematics,FacultyofMathematics,
[ WestUniversityofTimi¸soara,Bd. V.Parvan,nr. 4,300223, Timi¸soara,Romania,
1 E-mail: [email protected]
v c Department ofBiophysicsandMedicalInformatics,
1
UniversityofMedicineandPharmacy,PiataEftimieMurgu,nr. 3,300041, Timi¸soara,Romania,
9
5 E-mail: [email protected]
1
0
6
Abstract
0
/
h
t Theobjective of thispaperistoinvestigate thestability oflimitcycles
a
m of a mathematical model with a distributed delay which describes the
interaction between p53 and mdm2. Choosing the delay as a bifurca-
:
v
tion parameter we study the direction and stability of the bifurcating
i
X periodic solutions using the normal form and the center manifold the-
r orem. Some numerical examples are finally made in order to confirm
a
the theoretical results.
Keywords: delay differential equation, stability, Hopf bifurcation, P53, MDM2.
2000 AMS Mathematics Subject Classification: 34C23, 34C25, 37G05, 37G15,
92D10 .
∗
Corresponding author
1
1. Introduction
The tumor suppresser gene p53 plays a key role in oncogenesis and its
anomalies are almost universal in tumoral cells [4]. To make the reading
easier we briefly present the role of p53 gene in living cells. Normally the
activity of p53 is inhibited and the concentration of p53 protein are kept at
very low levels. The activation occurs when there are DNA damages [2,3].
Depending on these damages there are two outcomes: one is the cell cycle
arrest induced by a low level or a brief elevation of p53 protein, and the
other is the apoptosis induced by a high level or a prolonged elevation of p53
protein [5]. Knowing these 2 different outcomes, it is now clear that the level
of p53 protein should be kept tight control. This control is achieved with the
help of mdm2 gene with which p53 makes a feedback loop [5, 10]. Recently it
had been discovered that this loop is not quite straightforward because there
are two isoforms of the mdm2 protein, i.e. p76mdm2 and p90mdm2 which
have different features and roles [7].
This new model is based on [8, 9]. Here we achieve a smoother modeling
of the phenomenon, i.e. the interaction p53-mdm2. The production of p53
protein is continuous, so is the binding between p53 and the promoter of the
mdm2. The difference from the previous model [9] lies in the introduction of
the integral form in the third equation, which is the natural way of modeling
a continuous process. To have a clear picture of this aspect, we explain the
needfortheintroductionoftheintegralform,bythefactthatthesynthesis of
p53 protein and the binding p53-mdm2 are continuous and there is no buffer
to store p53 and then, at a given moment, to release it, and after that all the
quantity of p53 to bind mdm2 promoter. As a matter of fact, the molecules
of p53 protein that are bound with mdm2 promoter were not synthesized at
the moment t, but they were synthesized at a previous different moments,
they were bound at different moments.
The variables of the model are: x , x mRNA concentrations and y , y
1 2 1 2
the protein concentrations.
2
The mathematical system which describes our model is:
x˙ (t) = 1 b x (t),
1 1 1
−
y˙ (t) = x (t) (a +a y (t))y (t),
1 1 1 12 2 1
−
1 τ (1)
x˙ (t) = (αf(y (t))+(1 α)f(y (t s)))ds b x (t),
2 τ Z0 1 − 1 − − 2 2
y˙ (t) = x (t) (a +a y (t))y (t)
2 2 2 12 1 2
−
where: b ,b are the rates for mRNA degradation, a ,a ,a are the rates for
1 2 1 2 12
protein degradation. The function f : IR IR, is the Hill function, given
+
→
by:
xn
f(x) = (2)
a+xn
∗
with n IN ,a > 0. The parameters a , a , b , b , a of the model are
1 2 1 2 12
∈
assumed to be positive numbers less or equal to 1, α [0,1] and τ 0.
∈ ≥
For α = 0 in the present model, we obtain the model from [9], which
suggests that there is an oscillatory behavior based on using only numerical
simulations.
For the study of the model (1) we consider the following initial values:
x (0) = x¯ ,y (θ) = ϕ (θ),θ [ τ,0],x (0) = x¯ ,y (0) = y¯ ,
1 1 1 1 2 2 2 2
∈ −
with x¯ 0, x¯ 0, y¯ 0, ϕ (θ) 0, for all θ [ τ,0] and ϕ is a
1 2 2 1 1
≥ ≥ ≥ ≥ ∈ −
differentiable function.
The paper is organized as follows. In section 2, we discuss the local sta-
bility for the equilibrium state of system (1). We investigate the existence
of the Hopf bifurcation for system (1) using time delay as the bifurcation
parameter. In section 3, the direction of Hopf bifurcation is analyzed by the
normal form theory and the center manifold theorem introduced by Hassard
[6]. Numerical simulations for justifying the theoretical results are illustrated
in section 4. Finally, some conclusions are made and further research direc-
tions are presented.
2. Local stability and the existence of the Hopf bifurcation.
The equilibrium points of system (1) are given by the solutions of the
following system of equations:
3
1 b x (t) = 0,
1 1
−
x (t) (a +a y (t))y (t) = 0,
1 − 1 12 2 1 (3)
f(y (t)) b x (t) = 0,
1 2 2
−
x (t) (a +a y (t))y (t) = 0.
2 2 12 1 2
−
Let g : (0, ) IR be the function given by:
∞ →
(a +a x)(1 a b x)
2 12 1 1
g(x) = − . (4)
b a x
1 12
From (3) and (4) it results that one solution of system (3) is:
1 (a +a y )( a b y +1) a b y +1
2 12 1 1 1 1 1 1 1
x = , y , x = − , y = − , (5)
1 1 2 2
b b a y b a y
1 1 12 1 1 12 1
where y is the solution of the equation:
1
f(x) g(x) = 0. (6)
−
Because f(x) is an increasing function with lim f(x) = 1 then g(x) is an
x→∞
1
increasing function too, for x < . It results that equation (6) has a
a b
1 1
unique solution.
Proposition 1. If y is a real solution of (6) then the equilibrium point
10
of system (1) is:
1 1 a b y (a +a y )(1 a b y )
1 1 10 2 12 10 1 1 10
x = , y = − , x = − . (7)
10 20 20
b b a y b a y
1 1 12 10 1 12 10
We consider the following translation:
x = u +x ,y = u +y ,x = u +x ,y = u +y . (8)
1 1 10 1 2 10 2 3 20 2 4 20
With respect to (8), system (1) can be expressed as:
u˙ (t) = b u (t),
1 1 1
−
u˙ (t) = u (t) (a +a y )u (t) a y u (t) a u (t)u (t),
2 1 1 12 20 2 12 10 4 12 2 4
− − −
1 τ (9)
u˙ (t)= (αf(u (t)+y )+(1 α)f(u (t s)+y ))ds b u (t) b x ,
3 τ Z0 2 10 − 2 − 10 − 2 3 − 2 20
u˙ (t) = a y u (t)+u (t) (a +a y )u (t) a u (t)u (t).
4 12 20 2 3 2 12 10 4 12 2 4
− − −
4
System (9) has (0,0,0,0) as equilibrium point.
To investigate the local stability of the equilibrium state we linearize
system (9). We expand it in a Taylor series around the origin and neglect
the terms of higher order than the first order for the functions from the right
side of (9). We obtain:
1 τ
U˙(t) = AU(t)+ B U(t s)ds, (10)
τ Z0 −
where
b 0 0 0 0 0 0 0
1
−
1 (a +a y ) 0 a y 0 0 0 0
A= 00 − −1aα1ρ211y22020 −1b2 −(−a2+10a2121y010),B=00 (1−0α)ρ1 00 00 (11)
withρ = f′(y ), U(t) = (u (t),u (t),u (t),u (t))T, τ U(t s)ds = ( τ u (t
1 10 1 2 3 4 0 − 0 1 −
R R
s)ds, τ u (t s)ds, τ u (t s)ds, τ u (t s)ds)T.
0 2 − 0 3 − 0 4 −
R R R
The characteristic equation corresponding to system (10) is det(λI A
− −
1 τ
( e−λsds)B) = 0 which leads to:
τ Z
0
h τ
(λ+b )(λ3 +bλ2 +cλ+d+ e−λsds) = 0, (12)
1 τ Z
0
where
b = a +a +b +a (y +y ),
1 2 2 12 20 10
c = b (a +a )+b a (y +y )+a a +a (a y +a y ),
2 1 2 2 12 20 10 1 2 12 1 10 2 20 (13)
d = b a a +a b (y a +a y )+αa y ρ ,
2 1 2 12 2 20 2 1 10 12 10 1
h = (1 α)a y ρ .
12 10 1
−
The equilibrium point X∗ = (x ,y ,x ,y )T is locally asymptotically
10 10 20 20
stable if and only if all eigenvalues of (12) have negative real parts.
Because b > 0 we will analyze the function:
1
h τ
∆(λ,τ) = λ3 +bλ2 +cλ+d+ e−λsds, (14)
τ Z
0
5
with λ IR.
We ∈are going to show that the equilibrium point X∗ undergoes a Hopf
bifurcation. In this sense, we look for the existence of the purely imaginary
roots of ∆(λ,τ) = 0. First, we verify if X∗ is locally asymptotically stable
when τ = 0. In this case, the equation ∆(λ,τ) = 0 becomes:
λ3 +bλ2 +cλ+d+h = 0. (15)
Because the coefficients of equation (15) are positive then according to the
Routh-Hurwitz criterion we have:
Proposition 2. When there is no delay, the equilibrium point X∗ of
system (1) is locally asymptotically stable if and only if
cb > d+h,
where c,b,d,h are given by (13).
We are looking for the values τ so that the equilibrium point X∗ changes
0
from local asymptotic stability to instability or vice versa. This is specific for
the characteristic equation with pure imaginary solutions. Let λ = iω be
±
these solutions. We assume ω > 0. It is sufficient to look for λ = iω root of
∆(λ,τ) = 0. Separating real and imaginary parts of ∆(iω,τ) = 0 we obtain:
τω(bω2 d) τω2(c ω2)
sin(ωτ) = − , cos(ωτ) = 1 − . (16)
h − h
A solution of (16) is a pair (ω,τ) so that ωτ [0,2π], ω (0,√c] and
∈ ∈
τ = g (ω), where g : [0,√c] IR is given by:
1 1 +
−→
2h(c x2)
g (x) = − . (17)
1 (bx2 d)2 +x2(c x2)2
− −
2hc 2hc
From (17) it results that g (x) [0, ] and for all τ [0, ] there is
1 ∈ d2 ∈ d2
ω [0,√c] so that τ = g(ω).
∈
We have:
2hc
Proposition 3. For τ [0, ] there is ω [0,√c], so that (ω,τ) is
∈ d2 ∈
solution of (16).
According to the above proposition, system (16) has not unique solution.
6
In order to show that X∗ undergoes a Hopf bifurcation for τ = τ we
0
have to prove that iω are simple eigenvalues of ∆(λ,τ ) and satisfy the
0 0
±
dRe(λ)
transversality condition = 0.
dτ |τ=τ0 6
From (14), it results that:
∂∆ h
∆ (λ,τ) = (λ,τ) = 3λ2 +2bλ+c (1 (1+λτ)e−λτ). (18)
λ ∂λ − λ2τ −
From (18), for (ω,τ) a solution of (16) and τ = g(ω) it results that:
M = Re(∆ (iω,g (ω)) = 4ω2 +2c τ(bω2 d)
1 λ 1
− − −
3bω2 d h
M = Im(∆ (iω,g (ω)) = − − +τω(c ω2).
2 λ 1
ω −
Then, we have:
(3bω2 d h)2
M2+M2 = (2c 4ω2)2+τ2(bω2 d)2+ − − +τ2ω2(c ω2)2+
1 2 − − ω2 −
+ 2τ((3bω2 d h)(c ω2) (bω2 d)(2c 4ω2)).
− − − − − −
(19)
From (19), it results that:
Proposition 4. If bc > d + h, d > h, bc > 3d h, then λ = iω is a
0
−
simple root for the equation ∆(λ,τ ) = 0.
0
Now, we consider a branch of eigenvalues λ(τ) = ν(τ)+iω(τ) of (14) so
that ν(τ ) = 0 and ω(τ ) = ω , where (ω ,τ ) is a solution of (16). Differen-
0 0 0 0 0
tiating equation ∆(λ(τ),τ) = 0 with respect ro τ, we obtain:
dλ ∆ (iω ,τ )
′ τ 0 0
λ(τ) = = = M(ω ,τ )+iN(ω ,τ )
dτ|τ=τ0,ω=ω0 −∆ (iω ,τ ) 0 0 0 0
λ 0 0
where
M N +M N
1 1 2 2
M(ω ,τ ) = (20)
0 0 − M2 +M2
1 2
with
1
N = [h (bω2 d) τω2(c ω2)]
1
τ − − − −
1
N = [ω(c ω2) τω(bω2 d)]
2
τ − − −
and
M N M N
1 2 2 1
N = − . (21)
M2 +M2
1 2
7
By direct calculation M = 0.
6 2hc
We can conclude that when τ [0, ] the characteristic equation
0 ∈ d2
∆(λ,τ ) = 0 has a unique pair of purely imaginary simple eigenvalues satis-
0
dReλ(τ)
fied = 0 and the real roots are negative. Consequently, a Hopf
dτ |τ=τ0 6
bifurcation occurs at X∗ when τ = τ . Moreover, applying Rouche’s theorem
0
we can verify that every eigenvalue of ∆(λ,τ) = 0 with τ < τ has negative
0
real part. It follows that X∗ is locally asymptotically stable for 0 τ < τ .
0
≤
These results are summed up in the following theorem:
Theorem 1. Assume that cb > d + h. Then there exists values τ
0
2hc ∈
[0, ] of time delay so that the equilibrium point X∗ is locally asymptoti-
d2
cally stable when τ [0,τ ) and becomes unstable when τ = τ throughout a
0 0
∈
Hopf bifurcation. In particular, the periodic solutions appear for system (1)
when τ = τ .
0
3. Direction and stability of the Hopf bifurcation
In the previous section, we obtain some conditions with guarantee that
system (1) undergoes Hopf bifurcation at τ = τ .
0
In this section, we study the direction, the stability and the period of the
bifurcating periodic solutions. The used method is based on the normal form
theory and the center manifold theorem introduced by Hassard [6].
For an interval I IR, we define the space of continuous functions as
⊆
C(I,K) = f : I K,f continuous , where K = IR4 orCl4. When I =
{ → }
[ τ,0], τ = τ + µ, µ > 0 sufficiently small, we set C = C([ τ,0],K).
0 µ
− −
Expanding the functions from the right side of system (9) in Taylor series
around (0,0,0,0)T it results that:
1 τ τ
X˙(t)=AX(t)+ B X(t s)ds+F(X(t), (αX(t)+(1 α)X(t s))ds) (22)
τ Z0 − Z0 − −
where
X(t) = (u (t),u (t),u (t),u (t))T,
1 2 3 4
τ τ τ τ τ
X(t s)ds = ( u (t s)ds, u (t s)ds, u (t s)ds, u (t s)ds)T,
Z0 − Z0 1 − Z0 2 − Z0 3 − Z0 4 −
8
F(X(t), τ(αX(t)+(1 α)X(t s))ds) = (0,F2(u (t),u (t)),
F3( τ(αuR0(t)+(1 α)u−(t s))−ds),F4(u (t),u (t)))2T, 4 (23)
0 2 − 2 − 2 4
R
F2(u (t),u (t))= a u (t)u (t),
2 4 12 2 4
−
1 τ
F3( τ(αu (t)+(1 α)u (t s))ds)= ρ (αu (t)+(1 α)u (t s))2ds+
R0 2 − 2 − 2τ 2Z0 2 − 2 −
1 τ
+ ρ (αu (t)+(1 α)u (t s))3ds,
6τ 3Z0 2 − 2 −
F4(u (t),u (t)) = a u (t)u (t),
2 4 12 2 4
−
ρ = f′′(y ), ρ = f′′′(y ) and A,B are given by (11).
2 10 3 10
For Φ C with K =Cl4 we define a linear operator:
µ
∈
1 0
L (Φ) = AΦ(0) B Φ(s)ds
µ − τ0 Z−τ0
and a nonlinear operator:
0
F (Φ) = (0,F2(Φ (0),Φ (0)),F3( Φ (s)ds),F4(Φ (0),Φ (0)))T.
µ 2 4 Z−τ 2 2 4
For Φ C1([ τ ,0],Cl4) we define:
0
∈ −
dΦ(θ)
, θ [ τ ,0)
dθ ∈ − 0
A(µ)Φ(θ) = 1 0
AΦ(0) B Φ(s)ds, θ = 0,
− τ0 Z−τ0
0, θ [ τ ,0)
0
R(µ)Φ(θ) = ∈ −
F (Φ), θ = 0
µ
and for Ψ C1([0,τ ],Cl∗4), we define the adjoint operator ∗ of by:
0
∈ A A
dΨ(s)
, s (0,τ ]
∗ − ds ∈ 0
A Ψ(s) = 1 τ0
AΨ(0)+ ( Ψ(θ)dθ)B, s = 0.
τ Z
0 0
Then, we can rewrite (22) in the following vector form:
X˙ = A(µ)X +R(µ)X (24)
t t t
9
where X = X(t + θ) for θ [ τ ,0]. For Φ C([ τ ,0],Cl∗4) and Ψ
t 0 0
C([0,τ ],Cl∗4) we define the fo∈llow−ing bilinear form∈: − ∈
0
0 θ 1 ξ
< Ψ(s),Φ(θ) >= Ψ¯(0)Φ(0) Ψ¯(ξ θ)B( Φ(ξ′)dξ′)dξdθ,
−Z−τ0Zξ=0 − τ0 Z0
s [0,τ ], θ [ τ ,0].
0 0
∈ ∈ − ∗
Then, it can verified that and (0) are adjoint operators with respect
A A
to this bilinear form.
In the light of the obtained results in the last section, we assume that
iω are eigenvalues of (0). Thus, they are also eigenvalues of ∗. We can
0
± A A
easily obtain:
Φ(θ) = veλ1θ, θ [ τ ,0] (25)
0
∈ −
where v = (v ,v ,v ,v )T,
1 2 3 4
v = 0,v = a y ,v = a2 y y (λ +a +a y )(λ +a +a y ),
1 2 12 10 3 12 10 20 − 1 1 12 20 1 2 12 10
v = (λ +a +a y )
4 1 1 12 20
−
is the eigenvector of (0) corresponding to λ = iω and
1 0
A
Ψ(s) = weλ1s, s [0,τ ]
0
∈
where w = (w ,w ,w ,w ),
1 2 3 4
1 d d d
2 3 4
w = ,w = ,w = ,w = ,
1 2 3 4
η¯ η¯ η¯ η¯
a y (b +λ ) a y (b +λ )
12 10 1 1 12 10 1 1
d = b +λ ,d = ,d =
2 1 1 3 4
−(λ +a +a y )(b +λ ) −λ +a +a y
1 2 12 10 2 1 1 2 12 10
(1 α)ρ
η = v d¯ +v d¯ +v d¯ d¯v − 1(τ λ 2+2e−λ1τ0 +λ τ e−λ1τ0)
2 2 3 3 4 4 − 3 2 τ λ3 0 1 − 1 0
0 1
is the eigenvector of ∗ corresponding to λ = iω .
2 0
WecanverifythatA: < Ψ(s),Φ(s) >= 1,< Ψ(−s),Φ¯(s) >=< Ψ¯(s),Φ(s) >=
0, < Ψ¯(s),Φ¯(s) >= 1.
Using the approach in [1], we next compute the coordinates to describe
the center manifold Ω at µ = 0. Let X = X(t + θ),θ [ τ ,0], be the
0 t 0
∈ −
solution of equation (24) when µ = 0 and
z(t) =< Ψ,X >, w(t,θ) = X (θ) 2Re z(t)Φ(θ) .
t t
− { }
10
