- Email: [email protected]
Hopf bifurcation in a DDE model of gene expression
Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242 www.elsevier.com/locate/cnsns
Short communication
Hopf bifurcation in a DDE model of gene expression Anael Verdugo a, Richard Rand
b,*
a
b
Center for Applied Mathematics, Cornell University, Ithaca, NY 14850, United States Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14850, United States Received 1 May 2006; accepted 2 May 2006 Available online 22 June 2006
Abstract We analyze a model of gene transcription and protein synthesis which has been previously presented in the biological literature. The model takes the form of an ODE (ordinary differential equation) coupled to a DDE (delay differential equation), the state variables being concentrations of messenger RNA and protein. Linear analysis gives a critical time delay beyond which a periodic motion is born in a Hopf bifurcation. Lindstedt’s method is applied to the nonlinear system, resulting in closed form approximate expressions for the amplitude and frequency of oscillation. A parameter study shows that the Hopf bifurcation may not occur if the rates of degradation are too large. Ó 2006 Elsevier B.V. All rights reserved. PACS: 02.30.Ks; 02.30.Oz; 82.39.Rt Keywords: Delay; DDE; Hopf bifurcation; Gene transcription
1. Introduction This work deals with a mathematical model of gene expression [4]. The biology of the problem may be described as follows: A gene, i.e., a section of the DNA molecule, is copied (transcribed) onto messenger RNA (mRNA), which diffuses out of the nucleus of the cell into the cytoplasm, where it enters a subcellular structure called a ribosome. In the ribosome the genetic code on the mRNA produces a protein (a process called translation). The protein then diffuses back into the nucleus where it represses the transcription of its own gene. Dynamically speaking, this process may result in a steady state equilibrium, in which case the concentrations of mRNA and protein are constant, or it may result in an oscillation. In this paper we analyze a simple model previously proposed in the biological literature [4], and we show that the transition between equilibrium and oscillation is a Hopf bifurcation. The model takes the form of two equations, one an ordinary differential equation (ODE) and the other a delayed differential equation (DDE). The delay is due to an observed time lag in the transcription process.
*
Corresponding author. Tel.: +1 6072557145; fax: +1 6072552011. E-mail address: [email protected] (R. Rand).
1007-5704/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2006.05.001
236
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Oscillations in biological systems with delay have been dealt with previously in [1–3]. Mahaffy [1] studied a system in which concentrations of mRNA and cell repressor are analyzed by varying several parameters, such as diffusivity and cell radii. Delay is introduced into the system and the model is linearized to find stability changes and associated critical delays which give rise to Hopf bifurcations. In a later study, Mahaffy et al. [2] investigated a transport mechanism in cells to obtain nutrients. Their model examined how the change in diffusivities and cell radii caused biochemical oscillatory responses in the concentrations of the nutrients. Their model was reduced to a system of DDEs and stability analysis was used to show that the system can undergo Hopf bifurcations for certain parameter values. In a more recent study, Mocek et al. [3] studied biochemical systems with delay. They approximated the DDE system with an ODE system by means of characterizing critical delays. In all of these works, the presence of Hopf bifurcations was indicated by the existence of a periodic solution in the linearized equations. In the present work we go beyond the linearized equations, and by considering nonlinear effects we are able to predict the amplitude and frequency of the resulting limit cycle, and its stability. The model equations investigated here involve the variables M(t), the concentration of mRNA, and P(t), the concentration of the associated protein [4] 0 1 _ ¼ am B M @
1þ
1 n C A lm M Pd P0
P_ ¼ ap M lp P
ð1Þ ð2Þ
where dots represent differentiation with respect to time t, and where we use the subscript d to denote a variable which is delayed by time T. Thus Pd = P(t T). The model constants are as given in [4]: am is the rate at which mRNA is transcribed in the absence of the associated protein, ap is the rate at which the protein is produced from mRNA in the ribosome, lm and lp are the rates of degradation of mRNA and of protein, respectively, P0 is a reference concentration of protein, and n is a parameter. We assume lm = lp = l. 2. Stability of equilibrium We begin by rescaling Eqs. (1) and (2). We set m ¼ aMm , p ¼ amPap , and p0 ¼ amP 0ap , giving m_ ¼ 1þ
1 n lm pd p0
p_ ¼ m lp
ð3Þ ð4Þ
Equilibrium points, (m*, p*), for (3) and (4) are found by setting m_ ¼ 0 and p_ ¼ 0 lm ¼ 1þ
1 n
ð5Þ
p p0
m ¼ lp
ð6Þ
Eliminating m* from Eqs. (5) and (6), we obtain an equation on p* ðp Þnþ1 þ pn0 p
pn0 ¼0 l2
ð7Þ
Next we define n and g to be deviations from equilibrium: n = n(t) = m(t) m*, g = g(t) = p(t) p*, and gd = g(t T). This results in the nonlinear system n_ ¼ 1þ
1
gd þp p0
g_ ¼ n lg
n lðm þ nÞ
ð8Þ ð9Þ
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
237
Expanding for small values of gd, Eq. (8) becomes n_ ¼ ln Kgd þ H 2 g2d þ H 3 g3d þ
ð10Þ
where K, H2 and H3 depend on p*, p0, and n as follows: n nb p ; where b ¼ K¼ 2 p0 p ð1 þ bÞ H2 ¼
bnðbn n þ b þ 1Þ
H3 ¼
ð11Þ ð12Þ
3
2ðb þ 1Þ p2 bnðb2 n2 4bn2 þ n2 þ 3b2 n 3n þ 2b2 þ 4b þ 2Þ 6ðb þ 1Þ4 p3
ð13Þ
Next we analyze the linearized system coming from Eqs. (10) and (9) n_ ¼ ln Kgd g_ ¼ n lg
ð14Þ ð15Þ
Stability analysis of Eqs. (14) and (15) shows that for T = 0 (no delay), the equilibrium point (m*, p*) is a stable spiral. Increasing the delay, T, in the linear system (14) and (15), will yield a critical delay, Tcr, such that for T > Tcr, (m*, p*) will be unstable, giving rise to a Hopf bifurcation. For T = Tcr the system (14) and (15) will exhibit a pair of pure imaginary eigenvalues ±xi corresponding to the solution nðtÞ ¼ B cosðxt þ /Þ gðtÞ ¼ A cos xt
ð16Þ ð17Þ
where A and B are the amplitudes of the g(t) and n(t) oscillations, and where / is a phase angle. Note that we have chosen the phase of g(t) to be zero without loss of generality. Then for values of delay T close to Tcr, T ¼ T cr þ D
ð18Þ
the nonlinear system (3) and (4) is expected to exhibit a periodic solution (a limit cycle) which can be written in the approximate form of Eqs. (16) and (17). Substituting Eqs. (16) and (17) into Eqs. (14) and (15) and solving for x and Tcr we obtain pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x ¼ K l2 ð19Þ pffiffiffiffiffiffiffiffi 2l Kl2 arctan K2l2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð20Þ T cr ¼ K l2 3. Lindstedt’s method We use Lindstedt’s Method [5,6] on Eqs. (10) and (9). We begin by changing the first order system into a second order DDE. This results in the following form: €g þ 2lg_ þ l2 g ¼ Kgd þ H 2 g2d þ H 3 g3d þ
ð21Þ
where K, H2 and H3 are defined by Eqs. (11)–(13). We introduce a small parameter via the scaling g ¼ u
ð22Þ 2
2
The detuning D of Eq. (18) is scaled like , D = d T ¼ T cr þ D ¼ T cr þ 2 d
ð23Þ
238
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Next we stretch time by replacing the independent variable t by s, where s ¼ Xt
ð24Þ
This results in the following form of Eq. (21): X2
d2 u du þ 2lX þ l2 u ¼ Kud þ H 2 u2d þ 2 H 3 u3d 2 ds ds
ð25Þ
where ud = u(s XT). We expand X in a power series in , omitting the O() term for convenience, since it turns out to be zero X ¼ x þ 2 k 2 þ
ð26Þ
Next we expand the delay term ud ud ¼ uðs XT Þ ¼ uðs ðx þ 2 k 2 þ ÞðT cr þ 2 dÞÞ
ð27Þ
2
¼ uðs xT cr ðk 2 T cr þ xdÞ þ Þ 2
0
ð28Þ 3
¼ uðs xT cr Þ ðk 2 T cr þ xdÞu ðs xT cr Þ þ Oð Þ
ð29Þ
Now we expand u(s) in a power series in uðsÞ ¼ u0 ðsÞ þ u1 ðsÞ þ 2 u2 ðsÞ þ
ð30Þ
Substituting and collecting terms, we find d2 u0 du0 þ Ku0 ðs xT cr Þ þ l2 u0 ¼ 0 þ 2lx ds2 ds d2 u1 du1 x2 2 þ 2lx þ Ku1 ðs xT cr Þ þ l2 u1 ¼ H 2 u20 ðs xT cr Þ ds ds d2 u2 du2 x2 2 þ 2lx þ Ku2 ðs xT cr Þ þ l2 u2 ¼ . . . ds ds x2
ð31Þ ð32Þ ð33Þ
where . . . stands for terms in u0 and u1, omitted here for brevity. We take the solution of the u0 equation as ^ cos s u0 ðsÞ ¼ A
ð34Þ
^ Next we substitute (34) into (32) and obtain the following where from Eqs. (17) and (22) we know A ¼ A. expression for u1: u1 ðsÞ ¼ m1 sin 2s þ m2 cos 2s þ m3 where m1 is given by the equation pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ 2 H 2 l K l2 ðl2 KÞð2l2 3KÞ 2A m1 ¼ Kð16l6 39Kl4 þ 18K 2 l2 þ 9K 3 Þ
ð35Þ
ð36Þ
and where m2 and m3 are given by similar equations, omitted here for brevity. We substitute Eqs. (34) and (35) into (33), and, after trigonometric simplifications have been performed, we equate to zero the coefficients of the resonant terms sin s and cos s. This yields the amplitude, A, of the limit cycle that was born in the Hopf bifurcation A2 ¼
P D Q
ð37Þ
where P ¼ 8K 2 ðl2 KÞðl2 þ KÞð16l6 39Kl4 þ 18K 2 l2 þ 9K 3 Þ Q ¼ Q0 T cr þ Q1
ð38Þ ð39Þ
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
239
and Q0 ¼ 48H 3 K 2 l8 þ 16H 22 Kl8 69H 3 K 3 l6 þ 32H 22 K 2 l6 63H 3 K 4 l4 162H 22 K 3 l4 þ 81H 3 K 5 l2 þ 108H 22 K 4 l2 þ 27H 3 K 6 þ 30H 22 K 5 Q1 ¼ 96H 3 Kl þ 64H 22 l9 138H 3 K 2 l7 16H 22 Kl7 126H 3 K 3 l5 308H 22 K 2 l5 þ 162H 3 K 4 l3 þ 296H 22 K 3 l3 þ 54H 3 K 5 l þ 12H 22 K 4 l
ð40Þ
9
ð41Þ
Eq. (39) depends on l, K, H2, H3, and Tcr. By using Eq. (20) we may express Eq. (39) as a function of l, K, H2, and H3 only. Removal of secular terms also yields a value for the frequency shift k2, where, from Eq. (26), we have X = x + 2k2 k2 ¼
R d Q
ð42Þ
where Q is given by (39) and pffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ K l2 Q 0
ð43Þ
An expression for the amplitude B of the periodic solution for n(t) (see Eq. (16)) may be obtained directly from Eq. (9) by writing n ¼ g_ þ lg, where g A cos xt. We find pffiffiffiffi B ¼ KA ð44Þ where K and A are given as in (11) and (37), respectively. 4. Numerical example Using the same parameter values as in [4] l ¼ 0:03=min; p0 ¼ 100;
n¼5
ð45Þ
we obtain p ¼ 145:9158;
m ¼ 4:3774 H 2 ¼ 6:2778 105 ; H 3 ¼ 6:4101 107 2p ¼ 114:5432 X ¼ 5:4854 102 ; X
ð46Þ
K ¼ 3:9089 103 ;
ð47Þ
T cr ¼ 18:2470;
ð48Þ
Fig. 1. Period of oscillation, 2 Xp , plotted as a function of delay T, where X is given by Eq. (52). The initiation of oscillation at T = Tcr = 18.2470 is due to a supercritical Hopf bifurcation, and is marked in the figure with a dot.
240
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Here the delay Tcr and the response period 2p/x are given in minutes. Substituting (46)–(48) into (37)–(44) yields the following equations: pffiffiffiffi A ¼ 27:0215 D ð49Þ k 2 ¼ 2:4512 103 d pffiffiffiffi B ¼ 1:6894 D
ð50Þ ð51Þ
Note that since Eq. (49) requires D > 0 for the limit cycle to exist, and since we saw in Eqs. (14) and (15) that the origin is unstable for T > Tcr, i.e., for D > 0, we may conclude that the Hopf bifurcation is supercritical, i.e., the limit cycle is stable. Multiplying (50) by 2 and substituting into (26) we obtain X ¼ 5:4854 102 2:4512 103 D
ð52Þ
where D = T Tcr = T 18.2470. Plotting the period, 2p , against the delay, T, yields the graph shown in X Fig. 1. These results are in agreement with those obtained by numerical integration of the original Eqs. (1) and (2) and with those presented in [4]. 5. Effect of changing parameters An advantage of the closed form approximate solution presented in this paper is that the effect of changes in parameters may be easily studied. In this section we present a few results obtained from our solution. The equilibrium concentration p* is determined by solving Eq. (7) for given values of l, p0, and n. Fig. 2 shows p* displayed as a function of l for p0 = 10, 50, 100, and 200. Here and in the following plots we follow [4] and take n = 5. We note from Eq. (20) that the quantity K 2l2 must be non-negative in order that Tcr > 0, that is in order that parameters corresponding to the Hopf bifurcation occur in a delay equation. (If Tcr < 0 then we would have a future equation, which is physically unreasonable.) Fig. 3 shows that this condition restricts the values of degradation rate l for a given value of reference concentration p0. For values of l which are greater than this critical value lcritical, the system will not exhibit a Hopf bifurcation and no oscillation result. qffiffiffi will pffiffiffi ffi P and D . Fig. 4 disFrom Eq. (37) we see that the amplitude A of protein oscillation is the product of qffiffiffi Q plays QP as a function of l for p0 = 10, 50, 100, and 200 and for n = 5. Note that the maximum permissible value of l depends on p0 as shown in Fig. 3. Eqs. (26), (42), (43), (19) and (23) give that X ¼ xð1 QQ0 DÞ, where X is the frequency of oscillation for delay T = Tcr + D, and x is the frequency of oscillation for delay T = Tcr. Fig. 5 displays QQ0 as a function
Fig. 2. The equilibrium concentration p* displayed as a function of l for p0 = 10, 50, 100 and 200 and for n = 5.
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
241
Fig. 3. Values of degradation rate l which are greater than lcritical correspond to negative values of Tcr and will prevent the system from oscillating. Here lcritical is shown to depend on the reference concentration p0.
Fig. 4. Eq. (37) shows that the amplitude A of protein oscillation is the product of for p0 = 10, 50, 100 and 200 and for n = 5.
qffiffiffi P Q
and
qffiffiffi pffiffiffiffi D. Here QP is displayed as a function of l
Fig. 5. Our solution gives that X ¼ xð1 QQ0 DÞ where X is the frequency of oscillation for delay T = Tcr + D and x is the frequency of oscillation for delay T = Tcr. Here QQ0 is displayed as a function of l for p0 = 10, 50, 100 and 200 and for n = 5.
242
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
of l for p0 = 10, 50, 100, and 200 and for n = 5. Note again that the maximum permissible value of l depends on p0 as shown in Fig. 3. 6. Conclusions The biological literature [4] shows that long time behavior of gene expression dynamics can consist of both stable equilibrium as well as periodic behavior. By analyzing a DDE model originally proposed in [4], we have shown that the transition between these states is due to a Hopf bifurcation. Our nonlinear analysis provides approximate expressions for the amplitude and frequency of the resulting limit cycle as a function of the model parameters. Fig. 3 shows that the Hopf bifurcation may not occur if the rates of degradation l are too large. Inspection of Fig. 4 shows that for a given detuning D off of the Hopf bifurcation, the amplitude of the oscillation depends on both p0 and l. We see that increasing p0 for a fixed value of l causes an increase in amplitude. However, for a fixed value of p0, the limit cycle amplitude is largest for a certain optimal value of l. Fig. 5 shows a similar behavior regarding the period of the limit cycle oscillation. Here again, for a fixed value of p0 we see that the quantity QQ0 achieves a maximum for a certain optimal value of l. In this case the peak values of QQ0 correspond to minimal values of frequency X of the limit cycle, and thus to maximal values for the period of the limit cycle. References [1] Mahaffy JM. Genetic control models with diffusion and delays. Math Biosci 1988;90:519–33. [2] Mahaffy JM, Jorgensen DA, Vanderheyden RL. Oscillations in a model of repression with external control. J Math Biol 1992;30:669–91. [3] Mocek WT, Rudbicki R, Voit EO. Approximation of delays in biochemical systems. Math Biosci 2005;198:190–216. [4] Monk NAM. Oscillatory expression of Hes1, p53, and NF-jB driven by transcriptional time delays. Curr Biol 2003;13:1409–13. [5] Rand RH. Lecture notes on nonlinear vibrations (version 52). Available from: http://www.tam.cornell.edu/randdocs/nlvibe52.pdf (2005). [6] Rand R, Verdugo A. Hopf bifurcation formula for first order differential-delay equations. Commun Nonlinear Sci Numer Simul, in press, doi:10.1016/j.cnsns.2005.08.005.
Short communication
Hopf bifurcation in a DDE model of gene expression Anael Verdugo a, Richard Rand
b,*
a
b
Center for Applied Mathematics, Cornell University, Ithaca, NY 14850, United States Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14850, United States Received 1 May 2006; accepted 2 May 2006 Available online 22 June 2006
Abstract We analyze a model of gene transcription and protein synthesis which has been previously presented in the biological literature. The model takes the form of an ODE (ordinary differential equation) coupled to a DDE (delay differential equation), the state variables being concentrations of messenger RNA and protein. Linear analysis gives a critical time delay beyond which a periodic motion is born in a Hopf bifurcation. Lindstedt’s method is applied to the nonlinear system, resulting in closed form approximate expressions for the amplitude and frequency of oscillation. A parameter study shows that the Hopf bifurcation may not occur if the rates of degradation are too large. Ó 2006 Elsevier B.V. All rights reserved. PACS: 02.30.Ks; 02.30.Oz; 82.39.Rt Keywords: Delay; DDE; Hopf bifurcation; Gene transcription
1. Introduction This work deals with a mathematical model of gene expression [4]. The biology of the problem may be described as follows: A gene, i.e., a section of the DNA molecule, is copied (transcribed) onto messenger RNA (mRNA), which diffuses out of the nucleus of the cell into the cytoplasm, where it enters a subcellular structure called a ribosome. In the ribosome the genetic code on the mRNA produces a protein (a process called translation). The protein then diffuses back into the nucleus where it represses the transcription of its own gene. Dynamically speaking, this process may result in a steady state equilibrium, in which case the concentrations of mRNA and protein are constant, or it may result in an oscillation. In this paper we analyze a simple model previously proposed in the biological literature [4], and we show that the transition between equilibrium and oscillation is a Hopf bifurcation. The model takes the form of two equations, one an ordinary differential equation (ODE) and the other a delayed differential equation (DDE). The delay is due to an observed time lag in the transcription process.
*
Corresponding author. Tel.: +1 6072557145; fax: +1 6072552011. E-mail address: [email protected] (R. Rand).
1007-5704/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2006.05.001
236
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Oscillations in biological systems with delay have been dealt with previously in [1–3]. Mahaffy [1] studied a system in which concentrations of mRNA and cell repressor are analyzed by varying several parameters, such as diffusivity and cell radii. Delay is introduced into the system and the model is linearized to find stability changes and associated critical delays which give rise to Hopf bifurcations. In a later study, Mahaffy et al. [2] investigated a transport mechanism in cells to obtain nutrients. Their model examined how the change in diffusivities and cell radii caused biochemical oscillatory responses in the concentrations of the nutrients. Their model was reduced to a system of DDEs and stability analysis was used to show that the system can undergo Hopf bifurcations for certain parameter values. In a more recent study, Mocek et al. [3] studied biochemical systems with delay. They approximated the DDE system with an ODE system by means of characterizing critical delays. In all of these works, the presence of Hopf bifurcations was indicated by the existence of a periodic solution in the linearized equations. In the present work we go beyond the linearized equations, and by considering nonlinear effects we are able to predict the amplitude and frequency of the resulting limit cycle, and its stability. The model equations investigated here involve the variables M(t), the concentration of mRNA, and P(t), the concentration of the associated protein [4] 0 1 _ ¼ am B M @
1þ
1 n C A lm M Pd P0
P_ ¼ ap M lp P
ð1Þ ð2Þ
where dots represent differentiation with respect to time t, and where we use the subscript d to denote a variable which is delayed by time T. Thus Pd = P(t T). The model constants are as given in [4]: am is the rate at which mRNA is transcribed in the absence of the associated protein, ap is the rate at which the protein is produced from mRNA in the ribosome, lm and lp are the rates of degradation of mRNA and of protein, respectively, P0 is a reference concentration of protein, and n is a parameter. We assume lm = lp = l. 2. Stability of equilibrium We begin by rescaling Eqs. (1) and (2). We set m ¼ aMm , p ¼ amPap , and p0 ¼ amP 0ap , giving m_ ¼ 1þ
1 n lm pd p0
p_ ¼ m lp
ð3Þ ð4Þ
Equilibrium points, (m*, p*), for (3) and (4) are found by setting m_ ¼ 0 and p_ ¼ 0 lm ¼ 1þ
1 n
ð5Þ
p p0
m ¼ lp
ð6Þ
Eliminating m* from Eqs. (5) and (6), we obtain an equation on p* ðp Þnþ1 þ pn0 p
pn0 ¼0 l2
ð7Þ
Next we define n and g to be deviations from equilibrium: n = n(t) = m(t) m*, g = g(t) = p(t) p*, and gd = g(t T). This results in the nonlinear system n_ ¼ 1þ
1
gd þp p0
g_ ¼ n lg
n lðm þ nÞ
ð8Þ ð9Þ
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
237
Expanding for small values of gd, Eq. (8) becomes n_ ¼ ln Kgd þ H 2 g2d þ H 3 g3d þ
ð10Þ
where K, H2 and H3 depend on p*, p0, and n as follows: n nb p ; where b ¼ K¼ 2 p0 p ð1 þ bÞ H2 ¼
bnðbn n þ b þ 1Þ
H3 ¼
ð11Þ ð12Þ
3
2ðb þ 1Þ p2 bnðb2 n2 4bn2 þ n2 þ 3b2 n 3n þ 2b2 þ 4b þ 2Þ 6ðb þ 1Þ4 p3
ð13Þ
Next we analyze the linearized system coming from Eqs. (10) and (9) n_ ¼ ln Kgd g_ ¼ n lg
ð14Þ ð15Þ
Stability analysis of Eqs. (14) and (15) shows that for T = 0 (no delay), the equilibrium point (m*, p*) is a stable spiral. Increasing the delay, T, in the linear system (14) and (15), will yield a critical delay, Tcr, such that for T > Tcr, (m*, p*) will be unstable, giving rise to a Hopf bifurcation. For T = Tcr the system (14) and (15) will exhibit a pair of pure imaginary eigenvalues ±xi corresponding to the solution nðtÞ ¼ B cosðxt þ /Þ gðtÞ ¼ A cos xt
ð16Þ ð17Þ
where A and B are the amplitudes of the g(t) and n(t) oscillations, and where / is a phase angle. Note that we have chosen the phase of g(t) to be zero without loss of generality. Then for values of delay T close to Tcr, T ¼ T cr þ D
ð18Þ
the nonlinear system (3) and (4) is expected to exhibit a periodic solution (a limit cycle) which can be written in the approximate form of Eqs. (16) and (17). Substituting Eqs. (16) and (17) into Eqs. (14) and (15) and solving for x and Tcr we obtain pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x ¼ K l2 ð19Þ pffiffiffiffiffiffiffiffi 2l Kl2 arctan K2l2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð20Þ T cr ¼ K l2 3. Lindstedt’s method We use Lindstedt’s Method [5,6] on Eqs. (10) and (9). We begin by changing the first order system into a second order DDE. This results in the following form: €g þ 2lg_ þ l2 g ¼ Kgd þ H 2 g2d þ H 3 g3d þ
ð21Þ
where K, H2 and H3 are defined by Eqs. (11)–(13). We introduce a small parameter via the scaling g ¼ u
ð22Þ 2
2
The detuning D of Eq. (18) is scaled like , D = d T ¼ T cr þ D ¼ T cr þ 2 d
ð23Þ
238
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Next we stretch time by replacing the independent variable t by s, where s ¼ Xt
ð24Þ
This results in the following form of Eq. (21): X2
d2 u du þ 2lX þ l2 u ¼ Kud þ H 2 u2d þ 2 H 3 u3d 2 ds ds
ð25Þ
where ud = u(s XT). We expand X in a power series in , omitting the O() term for convenience, since it turns out to be zero X ¼ x þ 2 k 2 þ
ð26Þ
Next we expand the delay term ud ud ¼ uðs XT Þ ¼ uðs ðx þ 2 k 2 þ ÞðT cr þ 2 dÞÞ
ð27Þ
2
¼ uðs xT cr ðk 2 T cr þ xdÞ þ Þ 2
0
ð28Þ 3
¼ uðs xT cr Þ ðk 2 T cr þ xdÞu ðs xT cr Þ þ Oð Þ
ð29Þ
Now we expand u(s) in a power series in uðsÞ ¼ u0 ðsÞ þ u1 ðsÞ þ 2 u2 ðsÞ þ
ð30Þ
Substituting and collecting terms, we find d2 u0 du0 þ Ku0 ðs xT cr Þ þ l2 u0 ¼ 0 þ 2lx ds2 ds d2 u1 du1 x2 2 þ 2lx þ Ku1 ðs xT cr Þ þ l2 u1 ¼ H 2 u20 ðs xT cr Þ ds ds d2 u2 du2 x2 2 þ 2lx þ Ku2 ðs xT cr Þ þ l2 u2 ¼ . . . ds ds x2
ð31Þ ð32Þ ð33Þ
where . . . stands for terms in u0 and u1, omitted here for brevity. We take the solution of the u0 equation as ^ cos s u0 ðsÞ ¼ A
ð34Þ
^ Next we substitute (34) into (32) and obtain the following where from Eqs. (17) and (22) we know A ¼ A. expression for u1: u1 ðsÞ ¼ m1 sin 2s þ m2 cos 2s þ m3 where m1 is given by the equation pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ 2 H 2 l K l2 ðl2 KÞð2l2 3KÞ 2A m1 ¼ Kð16l6 39Kl4 þ 18K 2 l2 þ 9K 3 Þ
ð35Þ
ð36Þ
and where m2 and m3 are given by similar equations, omitted here for brevity. We substitute Eqs. (34) and (35) into (33), and, after trigonometric simplifications have been performed, we equate to zero the coefficients of the resonant terms sin s and cos s. This yields the amplitude, A, of the limit cycle that was born in the Hopf bifurcation A2 ¼
P D Q
ð37Þ
where P ¼ 8K 2 ðl2 KÞðl2 þ KÞð16l6 39Kl4 þ 18K 2 l2 þ 9K 3 Þ Q ¼ Q0 T cr þ Q1
ð38Þ ð39Þ
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
239
and Q0 ¼ 48H 3 K 2 l8 þ 16H 22 Kl8 69H 3 K 3 l6 þ 32H 22 K 2 l6 63H 3 K 4 l4 162H 22 K 3 l4 þ 81H 3 K 5 l2 þ 108H 22 K 4 l2 þ 27H 3 K 6 þ 30H 22 K 5 Q1 ¼ 96H 3 Kl þ 64H 22 l9 138H 3 K 2 l7 16H 22 Kl7 126H 3 K 3 l5 308H 22 K 2 l5 þ 162H 3 K 4 l3 þ 296H 22 K 3 l3 þ 54H 3 K 5 l þ 12H 22 K 4 l
ð40Þ
9
ð41Þ
Eq. (39) depends on l, K, H2, H3, and Tcr. By using Eq. (20) we may express Eq. (39) as a function of l, K, H2, and H3 only. Removal of secular terms also yields a value for the frequency shift k2, where, from Eq. (26), we have X = x + 2k2 k2 ¼
R d Q
ð42Þ
where Q is given by (39) and pffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ K l2 Q 0
ð43Þ
An expression for the amplitude B of the periodic solution for n(t) (see Eq. (16)) may be obtained directly from Eq. (9) by writing n ¼ g_ þ lg, where g A cos xt. We find pffiffiffiffi B ¼ KA ð44Þ where K and A are given as in (11) and (37), respectively. 4. Numerical example Using the same parameter values as in [4] l ¼ 0:03=min; p0 ¼ 100;
n¼5
ð45Þ
we obtain p ¼ 145:9158;
m ¼ 4:3774 H 2 ¼ 6:2778 105 ; H 3 ¼ 6:4101 107 2p ¼ 114:5432 X ¼ 5:4854 102 ; X
ð46Þ
K ¼ 3:9089 103 ;
ð47Þ
T cr ¼ 18:2470;
ð48Þ
Fig. 1. Period of oscillation, 2 Xp , plotted as a function of delay T, where X is given by Eq. (52). The initiation of oscillation at T = Tcr = 18.2470 is due to a supercritical Hopf bifurcation, and is marked in the figure with a dot.
240
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
Here the delay Tcr and the response period 2p/x are given in minutes. Substituting (46)–(48) into (37)–(44) yields the following equations: pffiffiffiffi A ¼ 27:0215 D ð49Þ k 2 ¼ 2:4512 103 d pffiffiffiffi B ¼ 1:6894 D
ð50Þ ð51Þ
Note that since Eq. (49) requires D > 0 for the limit cycle to exist, and since we saw in Eqs. (14) and (15) that the origin is unstable for T > Tcr, i.e., for D > 0, we may conclude that the Hopf bifurcation is supercritical, i.e., the limit cycle is stable. Multiplying (50) by 2 and substituting into (26) we obtain X ¼ 5:4854 102 2:4512 103 D
ð52Þ
where D = T Tcr = T 18.2470. Plotting the period, 2p , against the delay, T, yields the graph shown in X Fig. 1. These results are in agreement with those obtained by numerical integration of the original Eqs. (1) and (2) and with those presented in [4]. 5. Effect of changing parameters An advantage of the closed form approximate solution presented in this paper is that the effect of changes in parameters may be easily studied. In this section we present a few results obtained from our solution. The equilibrium concentration p* is determined by solving Eq. (7) for given values of l, p0, and n. Fig. 2 shows p* displayed as a function of l for p0 = 10, 50, 100, and 200. Here and in the following plots we follow [4] and take n = 5. We note from Eq. (20) that the quantity K 2l2 must be non-negative in order that Tcr > 0, that is in order that parameters corresponding to the Hopf bifurcation occur in a delay equation. (If Tcr < 0 then we would have a future equation, which is physically unreasonable.) Fig. 3 shows that this condition restricts the values of degradation rate l for a given value of reference concentration p0. For values of l which are greater than this critical value lcritical, the system will not exhibit a Hopf bifurcation and no oscillation result. qffiffiffi will pffiffiffi ffi P and D . Fig. 4 disFrom Eq. (37) we see that the amplitude A of protein oscillation is the product of qffiffiffi Q plays QP as a function of l for p0 = 10, 50, 100, and 200 and for n = 5. Note that the maximum permissible value of l depends on p0 as shown in Fig. 3. Eqs. (26), (42), (43), (19) and (23) give that X ¼ xð1 QQ0 DÞ, where X is the frequency of oscillation for delay T = Tcr + D, and x is the frequency of oscillation for delay T = Tcr. Fig. 5 displays QQ0 as a function
Fig. 2. The equilibrium concentration p* displayed as a function of l for p0 = 10, 50, 100 and 200 and for n = 5.
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
241
Fig. 3. Values of degradation rate l which are greater than lcritical correspond to negative values of Tcr and will prevent the system from oscillating. Here lcritical is shown to depend on the reference concentration p0.
Fig. 4. Eq. (37) shows that the amplitude A of protein oscillation is the product of for p0 = 10, 50, 100 and 200 and for n = 5.
qffiffiffi P Q
and
qffiffiffi pffiffiffiffi D. Here QP is displayed as a function of l
Fig. 5. Our solution gives that X ¼ xð1 QQ0 DÞ where X is the frequency of oscillation for delay T = Tcr + D and x is the frequency of oscillation for delay T = Tcr. Here QQ0 is displayed as a function of l for p0 = 10, 50, 100 and 200 and for n = 5.
242
A. Verdugo, R. Rand / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 235–242
of l for p0 = 10, 50, 100, and 200 and for n = 5. Note again that the maximum permissible value of l depends on p0 as shown in Fig. 3. 6. Conclusions The biological literature [4] shows that long time behavior of gene expression dynamics can consist of both stable equilibrium as well as periodic behavior. By analyzing a DDE model originally proposed in [4], we have shown that the transition between these states is due to a Hopf bifurcation. Our nonlinear analysis provides approximate expressions for the amplitude and frequency of the resulting limit cycle as a function of the model parameters. Fig. 3 shows that the Hopf bifurcation may not occur if the rates of degradation l are too large. Inspection of Fig. 4 shows that for a given detuning D off of the Hopf bifurcation, the amplitude of the oscillation depends on both p0 and l. We see that increasing p0 for a fixed value of l causes an increase in amplitude. However, for a fixed value of p0, the limit cycle amplitude is largest for a certain optimal value of l. Fig. 5 shows a similar behavior regarding the period of the limit cycle oscillation. Here again, for a fixed value of p0 we see that the quantity QQ0 achieves a maximum for a certain optimal value of l. In this case the peak values of QQ0 correspond to minimal values of frequency X of the limit cycle, and thus to maximal values for the period of the limit cycle. References [1] Mahaffy JM. Genetic control models with diffusion and delays. Math Biosci 1988;90:519–33. [2] Mahaffy JM, Jorgensen DA, Vanderheyden RL. Oscillations in a model of repression with external control. J Math Biol 1992;30:669–91. [3] Mocek WT, Rudbicki R, Voit EO. Approximation of delays in biochemical systems. Math Biosci 2005;198:190–216. [4] Monk NAM. Oscillatory expression of Hes1, p53, and NF-jB driven by transcriptional time delays. Curr Biol 2003;13:1409–13. [5] Rand RH. Lecture notes on nonlinear vibrations (version 52). Available from: http://www.tam.cornell.edu/randdocs/nlvibe52.pdf (2005). [6] Rand R, Verdugo A. Hopf bifurcation formula for first order differential-delay equations. Commun Nonlinear Sci Numer Simul, in press, doi:10.1016/j.cnsns.2005.08.005.