# Bifurcations and Dynamics of the Rb-E2F Pathway Involving miR449

Authors: Lingling Li; Jianwei Shen
Journal: Complexity
 (2017)

---

## Abstract

We focused on the gene regulative network involving Rb-E2F pathway and microRNAs (miR449) and studied the influence of time delay on the dynamical behaviors of Rb-E2F pathway by using Hopf bifurcation theory. It is shown that under certain assumptions the steady state of the delay model is asymptotically stable for all delay values; there is a critical value under another set of conditions; the steady state is stable when the time delay is less than the critical value, while the steady state is changed to be unstable when the time delay is greater than the critical value. Thus, Hopf bifurcation appears at the steady state when the delay passes through the critical value. Numerical simulations were presented to illustrate the theoretical results.


---

## Body

## 1. Introduction

In the last few years, an increasing number of noncoding RNA (or microRNA) have been discovered to play central regulatory roles in gene regulation processes of prokaryotes and eukaryotes [1]. It has become evident that miRNAs regulate a variety of biological processes, and their expression is often deregulated in human malignancy. On one hand, miRNAs play roles in tumorigenesis by modulating oncogenic and tumor suppressor pathways. On the other hand, the expressions of miRNAs can be regulated by several oncogenic or tumor suppressor transcription factors [2]. In this paper, we focus on miR449, which can induce cell senescence and apoptosis and act as a tumor suppressor through regulating Rb/E2F activity [2, 3]. In recent years, large numbers of researches have focused on the mechanisms controlling cellular proliferation associated with human cancer regulated by Rb-E2F pathway experimentally [4–7]. Rb and E2F proteins play important roles in the regulation of cell division, cell growth, and programmed cell death by controlling the expressions of genes involved in these processes, which are best known for their regulation of the cell cycle at the G1/S transition [8]. As the first identified tumor suppressor gene [9], Rb is recognized to play a fundamental role in a signaling pathway that controls cell proliferation [10]. Rb regulates the transcription of genes that are essential for DNA replication and cell cycle progression by binding and inhibiting E2F transcription factors [11]. In the Rb-E2F pathway including negative feedback loops involving miR449, miR449 provides a twofold safety mechanism to avoid excessive E2F-induced proliferation by cell cycle arrest and apoptosis [12]. Mathematical models have been established to explain the nonlinear dynamical behaviors of the Rb-E2F pathway [12–14], which mainly concentrate on the stability and bifurcation of the deterministic systems, but not taking into account the effects of time delay.Time delay is one of the most important characteristics of gene regulation. In many cases, a gene regulates the expression of another gene by its products (RNAs or proteins). Since it takes time to generate those products and different processes need different amounts of time, time delayed regulation is ubiquitous in cellular processes [15]. Time-delay system is also called system with aftereffect or hereditary system [16, 17], and it is common in mathematical biology, such as population dynamics, the chemostat, neural network, blood cell maturation, transcriptional regulator dynamics, virus dynamics, and genetic network [18–23]. Recent papers have demonstrated that complex dynamical behaviors can arise as a consequence of time delays in biological systems. For example, in gene regulatory networks, time delays may lead to oscillations in protein levels and existing oscillations may become more robust [24, 25]. Oscillatory cellular dynamics, in particular periodic oscillation, plays an important role in maintaining homeostasis of living organisms. In addition, taking into account time delays in models of gene regulatory networks is often essential to capture the whole range of dynamic behaviors. For example, in experiments, a single self-repressed gene has been observed to display oscillatory behavior, but which cannot be deduced by models that ignore time delays. However, this oscillatory behavior is reproduced by a mathematical model including time delays. In addition, theoretical analysis that ignored time delays led to the erroneous conclusion that oscillations were not possible for this single gene [18, 26]. Many kinds of methods have been introduced to infer time delayed gene regulatory network [27–31]. Here, we mainly use the local linearization approach to analyze the nonlinear system.In this paper, firstly, we will introduce Rb-E2F pathway network involving miR449 modeled by time delayed differential equations. Secondly, we will study the dynamical behaviors of the model and derive sufficient conditions of the oscillation by using Hopf bifurcation theory. Particularly, we will also prove that there are periodic solutions under certain conditions. Finally, numerical simulations will be showed to illustrate the theoretical results.


## 2. A Simple Gene Regulatory Network Mediated by miR449

### 2.1. A Simple Mathematical Model of Gene Regulation with a Delayed Negative Feedback Loop

A simplified model (Figure1) of the Rb-E2F pathway begins with growth signals activating CycD. Initially, E2F is bound to and repressed by Rb, a tumor suppressor protein that is dysfunctional in several major cancers. CycD represses the repressor Rb and allows E2F to be turned on. Myc also induces E2F transcription. Subsequently, E2F activates the transcription of CycE, which forms a complex with Cdk2 to further remove Rb repression, establishing a positive feedback loop. E2F also activates its own transcription, constituting another positive feedback loop. An interesting addition to the regulatory mechanism of Rb-E2F network is the recent discovery that miR449 modulates E2F activity. It has been demonstrated that E2F strongly upregulates the expression of miR449. In turn, E2F is inhibited by miR449 through regulating different transcripts. On one hand, miR449 directly affects level of its target transcript Myc and therefore lowers E2F concentration. On the other hand, miR449 directly affects E2F inducer Cdk6 and CycE, thus forming negative feedback loops [12, 14]. Yao et al. [14] provided a mathematical model in the absence of miR449 and indicated that the Rb-E2F pathway acts as a bistable switch to convert signal inputs into all-or-none E2F responses. Yan et al. [12] gave another mathematical model and further investigated the stabilities and bifurcations of E2F, CycE, and miR449 in the participation of miR449. Our main work in this paper is considering the effects of time delays on the dynamic behaviors of the model including miR449. When time delays are taken into account, the time delayed differential equation model of the network including miR449 is described by the following system:(1)dx1dt=k11x3k12+x3x1k13+x1+k14x3k15+x3+k16x4x8k17+x8+k18x5x8k19+x8-k110x1-k111x6x1,dx2dt=k21x1k22+x1-k23x2-k24x2x3k25+x3-k26x2x4k27+x4-k28x2x5k29+x5,dx3dt=k31Sk32+S-k33x3-k34x2t-τ1x3k35+x3,dx4dt=k41Sk42+S+k43x3k44+x3-k45x4-k46x2t-τ2x4k47+x4,dx5dt=k51x1k52+x1-k53x5-k54x2t-τ3x5k55+x5,dx6dt=k61+k62x7k63+x7-k64x6x1-k65x4x8k78+x8-k67x5x6k68+x6-k69x6,dx7dt=k71x4x6k72+x6+k73x5x6k74+x6+k75x4x8k76+x8+k77x5x8k78+x8-k79x7k710+x7-k711x7,dx8dt=k81x6x1-k82x4x8k83+x8-k84x5x8k85+x8-k86x8.x1,x2,x3,x4,x5,x6,x7,x8 represent the concentrations of E2F, miR449, Myc, CycD-Cdk4/6 complex, CycE-Cdk2 complex, Rb, phosphorylated Rb, and Rb/E2F complex, respectively. And S is intensity of growth factor. In the following simulations, all the values of parameters are shown in Table 1 unless specified elsewhere.Table 1
Parameters for the model.
Rate constant
Value
k11
0.4
k12
0.15
k13
0.15
k14
0.003
k15
0.15
k16
18
k17
0.92
k18
18
k19
0.92
k110
0.25
k111
180
k21
1.4
k22
0.15
k23
0.02
k24
0.6
k25
0.15
k26
1
k27
0.92
k28
0.7
k29
0.92
k31
1
k32
0.5
k33
0.7
k34
0.6
k35
0.15
k41
0.45
k42
0.5
k43
0.03
k44
0.15
k45
1.5
k46
1
k47
0.92
k51
0.35
k52
0.15
k53
1.5
k54
0.7
k55
0.92
k61
0.18
k62
3.6
k63
0.01
k64
180
k65
18
k66
0.92
k67
18
k68
0.92
k69
0.06
k71
18
k72
0.92
k73
18
k74
0.92
k75
18
k76
0.92
k77
18
k78
0.92
k79
3.6
k710
0.01
k711
0.06
k81
180
k82
18
k83
0.92
k84
18
k85
0.92
k86
0.03
Note. The descriptions of the parameters in Table 1 are shown in [13, 15].Figure 1
Rb-E2F pathway mediated by miR449.


### 2.2. Oscillation Induced by Time Delay

In this subsection, we consider system (1) with (x1,x2,x3,x4,x5,x6,x7,x8) as state variables. The linearized system of (1) at equilibrium point (x10,x20,x30,x40,x50,x60) is as follows:(2)dxdt=A0x+B0xt-τ,x=x1,x2,x3,x4,x5,x6,x7,x8T.(3)A0=a110a13a14a15a160a18a21a22a23a24a2500000a330000000a43a440000a51000a55000a6100a64a65a66a670000a74a75a76a77a78a8100a84a85a860a88,B0=00000000000000000b320000000b420000000b52000000000000000000000000000000,where(4)a11=k11x30k12+x30k13+x10-k11x30x10k12+x30k13+x102-k110-k111x60,a13=k11x10k12+x30k13+x10-k11x30x10k13+x10k12+x302k14k15+x30-k14x30k15+x302,a14=k16x80k17+x80,a15=k18x80k19+x80,a16=-k111x10,a18=k16x40k17+x80-k16x40x80k17+x802+k18k50k19+x80-k18x50x80k19+x802,a21=k21k22+x10-k21x10k22+x102,a22=-k23-k24x30k25+x30-k26x40k27+x40-k28x50k29+x50,a23=-k24x20k25+x30+k24x20x30k25+x302,a24=-k26x20k27+x40+k26x20x40k27+x402,a25=-k28x20k29+x50+k28x20x50k29+x502,a33=-k33-k34x20k35+x30+k34x20x30k35+x302,a43=k43k44+x30-k43x30k44+x302,a44=-k45-k46x20k47+x40+k46x20x40k47+x402,a51=k51k52+x10-k51x10k52+x102,a53=-k53-k54x20k55+x50+k54x20x50k55+x502,a61=-k64x60,a64=-k65x60k66+x60,a65=-k67x60k68+x60,a66=-k64x10-k65x40k66+x60+k65x40x60k66+x602-k67x50k68+x60+k67x50x60k68+x602-k69,a67=k62k63+x70-k62x70k63+x702,a74=k71x60k72+x60+k75x80k76+x80,a75=k73x60k74+x60+k77x80k78+x80,a76=k71x40k72+x60-k71x40x60k72+x602+k73x50k74+x60-k73x50x60k74+x602,a77=-k79k710+x70+k79x70k710+x702-k711,a78=k75x40k76+x80-k75x40x80k76+x802+k77x50k78+x80-k77x50x80k78+x802,a81=k81x60,a84=-k82x80k83+x80,a85=-k84x80k85+x80,a86=k81x10,a88=-k82x40k83+x80+k82x40x80k83+x802-k84x50k85+x80+k84x50x80k85+x802-k86,b32=-k34x30k35+x30,b42=-k46x40k47+x40,b52=-k54x50k55+x50.Then we can obtain the characteristic equation of (2) at the equilibrium (x10,x20,x30,x40,x50,x60,x70,x80) as follows:(5)λI-A0-B0e-λτ=0,where I is the 8∗8 identity matrix, and the characteristic equation (5) has the following form:(6)λ8+A1λ7+A2λ6+A3λ5+A4λ4+A5λ3+A6λ2+A7λ+A8+B1λ6+B2λ5+B3λ4+B4λ3+B5λ2+B6λ+B7e-λτ1+C1λ6+C2λ5+C3λ4+C4λ3+C5λ2+C6λ+C7e-λτ2+E1λ6+E2λ5+E3λ4+E4λ3+E5λ2+E6λ+E7e-λτ3=0,where the values of A1–8,B1–7,C1–7,E1–7 are showed in the Appendix.If we assume thatτ1=τ2=τ3=τ, we will have (7)λ8+A1λ7+A2λ6+A3λ5+A4λ4+A5λ3+A6λ2+A7λ+A8+D1λ6+D2λ5+D3λ4+D4λ3+D5λ2+D6λ+D7e-λτ=0,where Dj=Bj+Cj+Ej,j=1,2,3,4,5,6,7.(1) Ifτ=0, (7) becomes (8)λ8+G1λ7+G2λ6+G3λ5+G4λ4+G5λ3+G6λ2+G7λ+G8=0,where(9)G1=A1,G2=A2+D1,G3=A3+D2,G4=A4+D3,G5=A5+D4,G6=A6+D5,G7=A7+D6,G8=A8+D7.According to the Routh-Hurwitz criterion, all roots of (8) have negative real parts if and only if all the subdeterminants in the diagonal are positive; that is, (10)H1=G1>0,H2=G11G3G2>0,H3=G110G3G2G1G5G4G3>0,H4=G1100G3G2G11G5G4G3G2G7G6G5G4>0,H5=G11000G3G2G110G5G4G3G2G1G7G6G5G4G30G8G7G6G5>0,H6=G110000G3G2G1100G5G4G3G2G11G7G6G5G4G3G20G8G7G6G5G4000G8G7G6>0,H7=G1100000G3G2G11000G5G4G3G2G110G7G6G5G4G3G2G10G8G7G6G5G4G3000G8G7G6G500000G8G7>0,H8=G8H7>0.(2) Ifτ>0, considering the transcendental equation (7), clearly iω(ω>0) is a root of (7) if and only if (11)-ω8+A2ω6-A4ω4+A6ω2-A8+iA1ω7-A3ω5+A5ω3-A7ω=-D1ω6+D3ω4-D5ω2+D7+iD2ω5-D4ω3+D6ωcos⁡ωτ-isin⁡ωτ.Separating the real and imaginary parts of (11), we have (12)-ω8+A2ω6-A4ω4+A6ω2-A8=-D1ω6+D3ω4-D5ω2+D7cos⁡ωτ+D2ω5-D4ω3+D6ωsin⁡ωτ,A1ω7-A3ω5+A5ω3-A7ω=D1ω6-D3ω4+D5ω2-D7sin⁡ωτ+D2ω5-D4ω3+D6ωcos⁡ωτ.Adding up the squares of both equations of (12), we have(13)ω16+F1ω14+F2ω12+F3ω10+F4ω8+F5ω6+F6ω4+F7ω2+F8=0,where(14)F1=-2A2+A12,F2=-D12-2A1A3+A22+2A4,F3=2A1A5-2A2A4-D22+A32-2A6+2D1D3,F4=-2A3A5-2A1A7-D32+2D2D4+2A2A6+2A8+A42-2D1D5,F5=-D42-2A4A6-2A2A8-2D2D6+2D3D5+2D1D7+2A3A7+A52,F6=-2D3D7+2D4D6+2A4A8-2A5A7+A62-D52,F7=-D62-2A6A8+2D5D7+A72,F8=-D72+A82.Let z=ω2; (13) becomes (15)z8+F1z7+F2z6+F3z5+F4z4+F5z3+F6z2+F7z+F8=0.Denote (16)hz=z8+F1z7+F2z6+F3z5+F4z4+F5z3+F6z2+F7z+F8.Lemma 1.
IfF8<0, (15) has at least one positive root.Proof.
Clearly,h(0)=F8<0, and limz→∞h(z)=∞. Hence, there exists z0∈(0,∞), so that h(z0)=0. This completes the proof.Lemma 2.
IfF8>0, the sufficient condition for (15) has positive roots being z′(7)>0 and h(z′(7))<0.Proof.
From (16), we have h′(z)=8z7+7F1z6+6F2z5+5F4z+4F4z3+3F5z2+2F6z+F7; suppose the equation h′(z)=0 has seven real roots and satisfies z′(1)<z′(2)<z′(3)<z′(4)<z′(5)<z′(6)<z′(7), and z′(7) is the local minimum value, if z′(7)>0 and h(z′(7))<0; there exists z∈(z′(7),∞), so that h(z)=0; this completes the proof.Lemma 3.
IfF8>0, according to the Routh-Hurwitz criterion, all the roots of (15) have negative real parts if and only if the subdeterminants in the diagonal are positive, that is, Δj>0, j=1,2,…,8.Suppose that (15) has positive roots; without loss of generality, we assume that it has eight positive roots, denoted by z1,z2,z3,z4,z5,z6,z7,z8, respectively. Hence, (13) has eight positive roots, say ω1=z1,ω2=z2,ω3=z3,ω4=z4,ω5=z5,ω6=z6,ω7=z7,ω8=z8.From (12), we can get (17)τkj=arccosR1R2+R3R4R22+R42,k=1,2,3,…,8,j=0,1,2,3,…,where(18)R1=-ωk8+A2ωk6-A4ωk4+A6ωk2-A8,R2=-D1ωk6+D3ωk4-D5ωk2+D7,R3=A1ωk7-A3ωk5+A5ωk3-A7ωk,R4=D2ωk5-D4ωk3+D6ωk.Define(19)τ0=τk0j0=min1≤k≤8j≥0τkj.Let λ(τ)=η(τ)+iω(τ) be the root of (7) satisfying η(τ0)=0,ω(τ0)=ω0.Lemma 4.
Consider the exponential polynomialP(λ,e-λτ1,…,e-λτm)=p1(λ)+p2(λ)e-λτ1+⋯+pm(λ)e-λτm, where τi>0(i=1,2,…,m) and pi(λ) is polynomial about λ. As (τ1,τ2,…,τm) vary, the sum of the orders of the zeros of P(λ,e-λτ1,…,e-λτm) on the open right half plane can change only if a zero appears on or cross the imaginary axis.
Then, we have the following theoretical results.Theorem 5.
Suppose that conditions (10) are satisfied.(i)
IfΔj>0(j=1,2,…,8), then all roots of (7) have negative real parts for all τ≥0; thus the steady state (x10,x20,x30,x40,x50,x60,x70,x80) of system (2) is absolutely stable.(ii)
IfF8<0 or F8≥0,z′(7)>0, and h(z′(7))<0, then all the roots of (7) have negative real parts when τ∈[0,τ0); thus the steady state of system (2) is asymptotically stable.(iii)
If the condition of (ii) is satisfied,τ=τ0, and h′(z0)≠0,z0=ω02, then ±iω0 is a pair of simple purely imaginary roots of (15) and all other roots have negative real parts. Moreover, dReλ(τ0)/dττ=τ0>0. Thus, system (2) exhibits the Hopf bifurcation at (x10,x20,x30,x40,x50,x60,x70,x80).


## 2.1. A Simple Mathematical Model of Gene Regulation with a Delayed Negative Feedback Loop

A simplified model (Figure1) of the Rb-E2F pathway begins with growth signals activating CycD. Initially, E2F is bound to and repressed by Rb, a tumor suppressor protein that is dysfunctional in several major cancers. CycD represses the repressor Rb and allows E2F to be turned on. Myc also induces E2F transcription. Subsequently, E2F activates the transcription of CycE, which forms a complex with Cdk2 to further remove Rb repression, establishing a positive feedback loop. E2F also activates its own transcription, constituting another positive feedback loop. An interesting addition to the regulatory mechanism of Rb-E2F network is the recent discovery that miR449 modulates E2F activity. It has been demonstrated that E2F strongly upregulates the expression of miR449. In turn, E2F is inhibited by miR449 through regulating different transcripts. On one hand, miR449 directly affects level of its target transcript Myc and therefore lowers E2F concentration. On the other hand, miR449 directly affects E2F inducer Cdk6 and CycE, thus forming negative feedback loops [12, 14]. Yao et al. [14] provided a mathematical model in the absence of miR449 and indicated that the Rb-E2F pathway acts as a bistable switch to convert signal inputs into all-or-none E2F responses. Yan et al. [12] gave another mathematical model and further investigated the stabilities and bifurcations of E2F, CycE, and miR449 in the participation of miR449. Our main work in this paper is considering the effects of time delays on the dynamic behaviors of the model including miR449. When time delays are taken into account, the time delayed differential equation model of the network including miR449 is described by the following system:(1)dx1dt=k11x3k12+x3x1k13+x1+k14x3k15+x3+k16x4x8k17+x8+k18x5x8k19+x8-k110x1-k111x6x1,dx2dt=k21x1k22+x1-k23x2-k24x2x3k25+x3-k26x2x4k27+x4-k28x2x5k29+x5,dx3dt=k31Sk32+S-k33x3-k34x2t-τ1x3k35+x3,dx4dt=k41Sk42+S+k43x3k44+x3-k45x4-k46x2t-τ2x4k47+x4,dx5dt=k51x1k52+x1-k53x5-k54x2t-τ3x5k55+x5,dx6dt=k61+k62x7k63+x7-k64x6x1-k65x4x8k78+x8-k67x5x6k68+x6-k69x6,dx7dt=k71x4x6k72+x6+k73x5x6k74+x6+k75x4x8k76+x8+k77x5x8k78+x8-k79x7k710+x7-k711x7,dx8dt=k81x6x1-k82x4x8k83+x8-k84x5x8k85+x8-k86x8.x1,x2,x3,x4,x5,x6,x7,x8 represent the concentrations of E2F, miR449, Myc, CycD-Cdk4/6 complex, CycE-Cdk2 complex, Rb, phosphorylated Rb, and Rb/E2F complex, respectively. And S is intensity of growth factor. In the following simulations, all the values of parameters are shown in Table 1 unless specified elsewhere.Table 1
Parameters for the model.
Rate constant
Value
k11
0.4
k12
0.15
k13
0.15
k14
0.003
k15
0.15
k16
18
k17
0.92
k18
18
k19
0.92
k110
0.25
k111
180
k21
1.4
k22
0.15
k23
0.02
k24
0.6
k25
0.15
k26
1
k27
0.92
k28
0.7
k29
0.92
k31
1
k32
0.5
k33
0.7
k34
0.6
k35
0.15
k41
0.45
k42
0.5
k43
0.03
k44
0.15
k45
1.5
k46
1
k47
0.92
k51
0.35
k52
0.15
k53
1.5
k54
0.7
k55
0.92
k61
0.18
k62
3.6
k63
0.01
k64
180
k65
18
k66
0.92
k67
18
k68
0.92
k69
0.06
k71
18
k72
0.92
k73
18
k74
0.92
k75
18
k76
0.92
k77
18
k78
0.92
k79
3.6
k710
0.01
k711
0.06
k81
180
k82
18
k83
0.92
k84
18
k85
0.92
k86
0.03
Note. The descriptions of the parameters in Table 1 are shown in [13, 15].Figure 1
Rb-E2F pathway mediated by miR449.


## 2.2. Oscillation Induced by Time Delay

In this subsection, we consider system (1) with (x1,x2,x3,x4,x5,x6,x7,x8) as state variables. The linearized system of (1) at equilibrium point (x10,x20,x30,x40,x50,x60) is as follows:(2)dxdt=A0x+B0xt-τ,x=x1,x2,x3,x4,x5,x6,x7,x8T.(3)A0=a110a13a14a15a160a18a21a22a23a24a2500000a330000000a43a440000a51000a55000a6100a64a65a66a670000a74a75a76a77a78a8100a84a85a860a88,B0=00000000000000000b320000000b420000000b52000000000000000000000000000000,where(4)a11=k11x30k12+x30k13+x10-k11x30x10k12+x30k13+x102-k110-k111x60,a13=k11x10k12+x30k13+x10-k11x30x10k13+x10k12+x302k14k15+x30-k14x30k15+x302,a14=k16x80k17+x80,a15=k18x80k19+x80,a16=-k111x10,a18=k16x40k17+x80-k16x40x80k17+x802+k18k50k19+x80-k18x50x80k19+x802,a21=k21k22+x10-k21x10k22+x102,a22=-k23-k24x30k25+x30-k26x40k27+x40-k28x50k29+x50,a23=-k24x20k25+x30+k24x20x30k25+x302,a24=-k26x20k27+x40+k26x20x40k27+x402,a25=-k28x20k29+x50+k28x20x50k29+x502,a33=-k33-k34x20k35+x30+k34x20x30k35+x302,a43=k43k44+x30-k43x30k44+x302,a44=-k45-k46x20k47+x40+k46x20x40k47+x402,a51=k51k52+x10-k51x10k52+x102,a53=-k53-k54x20k55+x50+k54x20x50k55+x502,a61=-k64x60,a64=-k65x60k66+x60,a65=-k67x60k68+x60,a66=-k64x10-k65x40k66+x60+k65x40x60k66+x602-k67x50k68+x60+k67x50x60k68+x602-k69,a67=k62k63+x70-k62x70k63+x702,a74=k71x60k72+x60+k75x80k76+x80,a75=k73x60k74+x60+k77x80k78+x80,a76=k71x40k72+x60-k71x40x60k72+x602+k73x50k74+x60-k73x50x60k74+x602,a77=-k79k710+x70+k79x70k710+x702-k711,a78=k75x40k76+x80-k75x40x80k76+x802+k77x50k78+x80-k77x50x80k78+x802,a81=k81x60,a84=-k82x80k83+x80,a85=-k84x80k85+x80,a86=k81x10,a88=-k82x40k83+x80+k82x40x80k83+x802-k84x50k85+x80+k84x50x80k85+x802-k86,b32=-k34x30k35+x30,b42=-k46x40k47+x40,b52=-k54x50k55+x50.Then we can obtain the characteristic equation of (2) at the equilibrium (x10,x20,x30,x40,x50,x60,x70,x80) as follows:(5)λI-A0-B0e-λτ=0,where I is the 8∗8 identity matrix, and the characteristic equation (5) has the following form:(6)λ8+A1λ7+A2λ6+A3λ5+A4λ4+A5λ3+A6λ2+A7λ+A8+B1λ6+B2λ5+B3λ4+B4λ3+B5λ2+B6λ+B7e-λτ1+C1λ6+C2λ5+C3λ4+C4λ3+C5λ2+C6λ+C7e-λτ2+E1λ6+E2λ5+E3λ4+E4λ3+E5λ2+E6λ+E7e-λτ3=0,where the values of A1–8,B1–7,C1–7,E1–7 are showed in the Appendix.If we assume thatτ1=τ2=τ3=τ, we will have (7)λ8+A1λ7+A2λ6+A3λ5+A4λ4+A5λ3+A6λ2+A7λ+A8+D1λ6+D2λ5+D3λ4+D4λ3+D5λ2+D6λ+D7e-λτ=0,where Dj=Bj+Cj+Ej,j=1,2,3,4,5,6,7.(1) Ifτ=0, (7) becomes (8)λ8+G1λ7+G2λ6+G3λ5+G4λ4+G5λ3+G6λ2+G7λ+G8=0,where(9)G1=A1,G2=A2+D1,G3=A3+D2,G4=A4+D3,G5=A5+D4,G6=A6+D5,G7=A7+D6,G8=A8+D7.According to the Routh-Hurwitz criterion, all roots of (8) have negative real parts if and only if all the subdeterminants in the diagonal are positive; that is, (10)H1=G1>0,H2=G11G3G2>0,H3=G110G3G2G1G5G4G3>0,H4=G1100G3G2G11G5G4G3G2G7G6G5G4>0,H5=G11000G3G2G110G5G4G3G2G1G7G6G5G4G30G8G7G6G5>0,H6=G110000G3G2G1100G5G4G3G2G11G7G6G5G4G3G20G8G7G6G5G4000G8G7G6>0,H7=G1100000G3G2G11000G5G4G3G2G110G7G6G5G4G3G2G10G8G7G6G5G4G3000G8G7G6G500000G8G7>0,H8=G8H7>0.(2) Ifτ>0, considering the transcendental equation (7), clearly iω(ω>0) is a root of (7) if and only if (11)-ω8+A2ω6-A4ω4+A6ω2-A8+iA1ω7-A3ω5+A5ω3-A7ω=-D1ω6+D3ω4-D5ω2+D7+iD2ω5-D4ω3+D6ωcos⁡ωτ-isin⁡ωτ.Separating the real and imaginary parts of (11), we have (12)-ω8+A2ω6-A4ω4+A6ω2-A8=-D1ω6+D3ω4-D5ω2+D7cos⁡ωτ+D2ω5-D4ω3+D6ωsin⁡ωτ,A1ω7-A3ω5+A5ω3-A7ω=D1ω6-D3ω4+D5ω2-D7sin⁡ωτ+D2ω5-D4ω3+D6ωcos⁡ωτ.Adding up the squares of both equations of (12), we have(13)ω16+F1ω14+F2ω12+F3ω10+F4ω8+F5ω6+F6ω4+F7ω2+F8=0,where(14)F1=-2A2+A12,F2=-D12-2A1A3+A22+2A4,F3=2A1A5-2A2A4-D22+A32-2A6+2D1D3,F4=-2A3A5-2A1A7-D32+2D2D4+2A2A6+2A8+A42-2D1D5,F5=-D42-2A4A6-2A2A8-2D2D6+2D3D5+2D1D7+2A3A7+A52,F6=-2D3D7+2D4D6+2A4A8-2A5A7+A62-D52,F7=-D62-2A6A8+2D5D7+A72,F8=-D72+A82.Let z=ω2; (13) becomes (15)z8+F1z7+F2z6+F3z5+F4z4+F5z3+F6z2+F7z+F8=0.Denote (16)hz=z8+F1z7+F2z6+F3z5+F4z4+F5z3+F6z2+F7z+F8.Lemma 1.
IfF8<0, (15) has at least one positive root.Proof.
Clearly,h(0)=F8<0, and limz→∞h(z)=∞. Hence, there exists z0∈(0,∞), so that h(z0)=0. This completes the proof.Lemma 2.
IfF8>0, the sufficient condition for (15) has positive roots being z′(7)>0 and h(z′(7))<0.Proof.
From (16), we have h′(z)=8z7+7F1z6+6F2z5+5F4z+4F4z3+3F5z2+2F6z+F7; suppose the equation h′(z)=0 has seven real roots and satisfies z′(1)<z′(2)<z′(3)<z′(4)<z′(5)<z′(6)<z′(7), and z′(7) is the local minimum value, if z′(7)>0 and h(z′(7))<0; there exists z∈(z′(7),∞), so that h(z)=0; this completes the proof.Lemma 3.
IfF8>0, according to the Routh-Hurwitz criterion, all the roots of (15) have negative real parts if and only if the subdeterminants in the diagonal are positive, that is, Δj>0, j=1,2,…,8.Suppose that (15) has positive roots; without loss of generality, we assume that it has eight positive roots, denoted by z1,z2,z3,z4,z5,z6,z7,z8, respectively. Hence, (13) has eight positive roots, say ω1=z1,ω2=z2,ω3=z3,ω4=z4,ω5=z5,ω6=z6,ω7=z7,ω8=z8.From (12), we can get (17)τkj=arccosR1R2+R3R4R22+R42,k=1,2,3,…,8,j=0,1,2,3,…,where(18)R1=-ωk8+A2ωk6-A4ωk4+A6ωk2-A8,R2=-D1ωk6+D3ωk4-D5ωk2+D7,R3=A1ωk7-A3ωk5+A5ωk3-A7ωk,R4=D2ωk5-D4ωk3+D6ωk.Define(19)τ0=τk0j0=min1≤k≤8j≥0τkj.Let λ(τ)=η(τ)+iω(τ) be the root of (7) satisfying η(τ0)=0,ω(τ0)=ω0.Lemma 4.
Consider the exponential polynomialP(λ,e-λτ1,…,e-λτm)=p1(λ)+p2(λ)e-λτ1+⋯+pm(λ)e-λτm, where τi>0(i=1,2,…,m) and pi(λ) is polynomial about λ. As (τ1,τ2,…,τm) vary, the sum of the orders of the zeros of P(λ,e-λτ1,…,e-λτm) on the open right half plane can change only if a zero appears on or cross the imaginary axis.
Then, we have the following theoretical results.Theorem 5.
Suppose that conditions (10) are satisfied.(i)
IfΔj>0(j=1,2,…,8), then all roots of (7) have negative real parts for all τ≥0; thus the steady state (x10,x20,x30,x40,x50,x60,x70,x80) of system (2) is absolutely stable.(ii)
IfF8<0 or F8≥0,z′(7)>0, and h(z′(7))<0, then all the roots of (7) have negative real parts when τ∈[0,τ0); thus the steady state of system (2) is asymptotically stable.(iii)
If the condition of (ii) is satisfied,τ=τ0, and h′(z0)≠0,z0=ω02, then ±iω0 is a pair of simple purely imaginary roots of (15) and all other roots have negative real parts. Moreover, dReλ(τ0)/dττ=τ0>0. Thus, system (2) exhibits the Hopf bifurcation at (x10,x20,x30,x40,x50,x60,x70,x80).


## 3. Numerical Analysis

In this section, we demonstrate the above theoretical results by numerical method. When we takeS=5,k21=1.4, and the other parameters are shown in Table 1, system (1) becomes(20)dx1dt=0.4x3x10.15+x30.15+x1+0.003x30.15+x3+18x4x80.92+x8+18x5x80.92+x8-0.25x1-180x6x1,dx2dt=1.4x10.15+x1-0.02x2-0.6x2x30.15+x3-x2x40.92+x4-0.7x2x50.92+x5,dx3dt=0.9090909090-0.7x3-0.6x2t-τ1x30.15+x3,dx4dt=0.4090909090+0.03x30.15+x3-1.5x4-x2t-τ2x40.92+x4,dx5dt=0.35x10.15+x1-1.5x5-0.7x2t-τ3x50.92+x5,dx6dt=0.18+3.6x70.01+x7-180x6x1-18x4x60.92+x6-18x5x60.92+x6-0.06x6,dx7dt=18x4x60.92+x6+18x5x60.92+x6+18x4x80.92+x8+18x5x80.92+x8-3.6x70.01+x7-0.06x7,dx8dt=180x6x1-18x4x80.92+x8-18x5x80.92+x8-0.03x8.The system has a positive equilibrium point Z=(0.4102,1.6027,0.3430,0.1429,0.0985,0.0478,1.2301, 3.4433). Using Theorem 5, there is a critical value of the time delay τ0=0.122632125. The equilibrium point is stable when τ<τ0 (see Figures 2 and 3); the equilibrium point becomes unstable and a Hopf bifurcation occurs when τ passes through the critical value τ0 (see Figure 4). The bifurcation diagrams of system (20) are shown in Figure 5, where the control parameter is the time delay τ.Figure 2
An asymptotically stable equilibrium forτ=τ1=τ2=τ3=0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=5,k21=1.4.Figure 3
An asymptotically stable equilibrium forτ=τ1=τ2=τ3<τ0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=5,k21=1.4.Figure 4
A periodic solution bifurcated from equilibrium forτ=τ1=τ2=τ3>τ0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=5,k21=1.4.Figure 5
Bifurcation induced by time delay withS=5,k21=1.4. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=5,k21=1.4.When we takeS=3,k21=2.5, and the other parameters are shown in Table 1, system (1) becomes(21)dx1dt=0.4x3x10.15+x30.15+x1+0.003x30.15+x3+18x4x80.92+x8+18x5x80.92+x8-0.25x1-180x6x1,dx2dt=1.4x10.15+x1-0.02x2-0.6x2x30.15+x3-x2x40.92+x4-0.7x2x50.92+x5,dx3dt=0.9090909090-0.7x3-0.6x2t-τ1x30.15+x3,dx4dt=0.4090909090+0.03x30.15+x3-1.5x4-x2t-τ2x40.92+x4,dx5dt=0.35x10.15+x1-1.5x5-0.7x2t-τ3x50.92+x5,dx6dt=0.18+3.6x70.01+x7-180x6x1-18x4x60.92+x6-18x5x60.92+x6-0.06x6,dx7dt=18x4x60.92+x6+18x5x60.92+x6+18x4x80.92+x8+18x5x80.92+x8-3.6x70.01+x7-0.06x7,dx8dt=180x6x1-18x4x80.92+x8-18x5x80.92+x8-0.03x8.The system has a positive equilibrium point Z′=(0.0561,1.0057,0.5478,0.1689,0.0427,0.2816,1.4875, 2.4618). Using Theorem 5, there is a critical value of the time delay τ0′=0.1334173317. The equilibrium point is stable when τ<τ0′ (see Figures 6 and 7); the equilibrium point becomes unstable and a Hopf bifurcation occurs when τ passes through the critical value τ0′ (see Figure 8). The bifurcation diagrams of system (21) are shown in Figure 9, where the control parameter is the time delay τ.Figure 6
An asymptotically stable equilibrium forτ=τ1=τ2=τ3=0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=3,k21=2.5.Figure 7
An asymptotically stable equilibrium forτ=τ1=τ2=τ3<τ0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=3,k21=2.5.Figure 8
A periodic solution bifurcated from equilibrium forτ=τ1=τ2=τ3>τ0. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=3,k21=2.5.Figure 9
Bifurcation induced by time delay withS=3,k21=2.5. Assume initial conditions are x1=1.2,x2=0,x3=0,x4=0,x5=0,x6=0.55,x7=0,x8=0, and S=3,k21=2.5.


## 4. Discussion

In this paper, we studied the dynamical behaviors of the Rb-E2F pathway including negative feedback loops involving miR449 by using Hopf bifurcation theory. On one hand, we gave the detailed theoretical analysis of system (1) in which there are some dynamical behaviors, and it is shown that, under certain conditions, the equilibrium point is asymptotically stable for all the delay τ≥0; however, if these conditions are not met, there is a critical delay τ0; when τ<τ0, the equilibrium point is asymptotically stable; when τ>τ0, periodic oscillations occur. Thus, Hopf bifurcation appears at the steady state when the delay passes through the critical value τ0. On the other hand, through numerical simulations, we calculated the critical value of τ0, giving the time courses of E2F, CycE, and miR449 both τ<τ0 and τ>τ0, drawing the bifurcation diagrams of E2F, CycE, and miR449, respectively, and testifying the correctness of the theory.



---

*Source: 1409865-2017-10-19.xml*