Dynamical Behaviors of Rb-E2F Pathway Including Negative Feedback Loops Involving miR449

MiRNAs, which are a family of small non-coding RNAs, regulate a broad array of physiological and developmental processes. However, their regulatory roles have remained largely mysterious. E2F is a positive regulator of cell cycle progression and also a potent inducer of apoptosis. Positive feedback loops in the regulation of Rb-E2F pathway are predicted and shown experimentally. Recently, it has been discovered that E2F induce a cluster of miRNAs called miR449. In turn, E2F is inhibited by miR449 through regulating different transcripts, thus forming negative feedback loops in the interaction network. Here, based on the integration of experimental evidence and quantitative data, we studied Rb-E2F pathway coupling the positive feedback loops and negative feedback loops mediated by miR449. Therefore, a mathematical model is constructed based in part on the model proposed in Yao-Lee et al. (2008) and nonlinear dynamical behaviors including the stability and bifurcations of the model are discussed. A comparison is given to reveal the implication of the fundamental differences of Rb-E2F pathway between regulation and deregulation of miR449. Coherent with the experiments it predicts that miR449 plays a critical role in regulating the cell cycle progression and provides a twofold safety mechanism to avoid excessive E2F-induced proliferation by cell cycle arrest and apoptosis. Moreover, numerical simulation and bifurcation analysis shows that the mechanisms of the negative regulation of miR449 to three different transcripts are quite distinctive which needs to be verified experimentally. This study may help us to analyze the whole cell cycle process mediated by other miRNAs more easily. A better knowledge of the dynamical behaviors of miRNAs mediated networks is also of interest for bio-engineering and artificial control.


Introduction
MicroRNAs (miRNAs) have been demonstrated to play crucial roles both in prokaryotes and eukaryotes [1,2]. Their importance is suggested in [1][2][3] by (i) the predictions that each miRNA targets hundreds of genes and that the majority of protein-coding genes are miRNA targets [4][5][6][7][8], (ii) their abundance, with some miRNAs expressed as highly as 50,000 copies per cell [9], and (iii) their sequence conservation, with some miRNAs conserved from sea urchins to humans [10]. MiRNAs can regulate a large variety of cellular processes, from differentiation and proliferation to apoptosis [11][12][13][14][15], by determining how and when genes turn on and off. Thus, from a biological point of view, miRNAs are challenging objects to study as they regulate cohorts of target genes, which are not readily identified. From a therapeutical point of view, miRNAs are highly interesting as several studies have demonstrated the power of miRNAs as biomarkers and initial preclinical studies have established that miRNAs may be therapeutically targeted in vivo [16].
MiRNAs regulate protein synthesis in the cell cytoplasm by promoting target mRNAs' degradation and/or inhibiting their translation. It is thought that the primary role of miRNAs is to modulate or fine-tune the dynamics of regulatory networks [17][18][19][20]. The significance of this role is now increasingly recognized as there are now many reported cases in which miRNA deregulation in a broad spectrum of diseases including all major cancers [21]. The regulatory roles of miRNAs have been a subject of research for the last several years, both experimentally and theoretically [6][7][8][9][10]. Although some of the miRNAs have been well studied, the information about possible functions and biological significance of miRNAs still need to be fully understood due to the diversity of mechanisms by which miRNAs may regulate biological processes. Here, we focus on miR449, which has been identified significantly down-regulated or lost in gastric cells, testicular cancer cells as well as in lung adenocarcinoma cell line [21][22][23][24]. Its reintroduction into these cancer cell lines leads to inhibition of cell proliferation and induces senescence and apoptosis by targeting different cell cycle regulators. Hence, these studies further underly that miR449 may act as a tumor suppressor [21,22].
The E2F family is best known for its critical role in cell proliferation. It is intimately connected to cell proliferation by coordinating a large group of genes involving G1 to S transition. In normal cells the E2F activity is tightly controlled at multiple levels. Deregulated E2F activity occurs in the vast majority of human tumors through several different mechanisms. These changes often lead to the inappropriate activation of E2F and its transcriptional target genes that further result in improper cell cycle control and unrestricted proliferation. Thus, repressing E2F activity is the key mechanism to exert the tumor-suppresive function [25]. Although the detailed mechanism of posttranscriptional regulation by miRNAs is not fully understood, evidence for the functional roles of miRNAs is accumulating. Recent study demonstrates that miRNAs mainly control transcripts coding for proteins involved in cell damage responses, cell cycle control, inflammation and cancer pathways [22]. Interestingly, it has been discovered that E2F strongly upregulate miR449 [21,22,26]. In turn, E2F is inhibited by miR449 through regulating different transcripts, thus forming negative feedback loops in the interaction network. This leads to following natural questions: what will occur by adding these negative feedback loops into Rb-E2F pathway and how does miR449 fine-tune the dynamics in this cancer network.
To address these questions in this paper, a mathematical model is constructed based in part on the model proposed in Yao-Lee et al. (2008) and quantitative comparison of dynamical features associated with or without miR449 regulation is given. They largely differ in the onset of bistability and oscillations and further affects E2F activity and G1/S transition. It is expected that the difference will generate a detailed and precise insight of mi449 mediated regulation. In agreement with experimental observations, the model verified that miR449 plays a critical role in regulating the cell cycle progression and provides a twofold safety mechanism to avoid excessive E2F-induced proliferation by cell cycle arrest and apoptosis. Moreover, bifurcation analysis of the model predicts that different mechanisms are employed by miR449 to negative regulate three different transcripts, which needs to be verified experimentally. These results are important for understanding the dynamics and functions of miRNAs and provide clues for therapeutic manipulation of Rb-E2F pathway in the treatment of cancer.

Results
In this research, particular attentions are payed to the bistability and oscillations. Bistability, the capacity to achieve two alternative internal states in response to different stimuli, exists ubiquitously in synthetic and natural biomolecular networks and has fundamental biological significance [27][28][29][30]. Cells can switch between two internal states to accommodate environmental and intercellular conditions owing to regulatory interactions among cellular components. Oscillatory behaviors are generated by complex interactions among genes, proteins and metabolites. They are used to control every aspect of cell physiology, from signalling, motility and development to growth, division and death [30][31][32]. Hopf bifurcations, including supercritical and subcritical Hopf bifucations, are the key principle of designing biochemical oscillators. In the following, numerical simulation and bifurcation analysis in our paper were carried out with XPPAUT. Moreover, XPPAUT codes are available in the Supporting Information File S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11, S12, S13, 14.

Modeling Rb-E2F pathway mediated by miR449
We focus on Rb-E2F pathway involved multiple positive feedback loops and negative feedback loops mediated by miR449. The network that we considered in this paper is shown schematically in Figure 1. E2F is a transcriptional factor that can activate genes encoding proteins involved in DNA replication and cell cycle progression [25,[33][34][35], and is essential in the regulation of cell proliferation [25,36,37]. In quiescent cells, E2F is bound to and repressed by retinoblastoma protein (Rb), a tumor suppressor protein that is dysfunctional in several major cancers. With sufficient growth stimulation, phosphorylation by Myc-induced cyclin D (CycD)-Cdk4,6 removes Rb repression; Myc also induces E2F transcription. Subsequently, E2F activates the transcription of CycE, which forms a complex with Cdk2 to further remove Rb repression by phosphorylation, establishing a positive-feedback loop. E2F also activates its own transcription, constituting another positive-feedback loop [25]. An intriguing addition to the already complex regulatory mechanism of Rb-E2F network is the recent discovery that miR449 modulate E2F activity [21,22,26,38]. It has been demonstrated that E2F strongly upregulate the expression of miR449. In turn, E2F is inhibited by miR449 through regulating different transcripts, thus forming negative feedback loops. One of the miR449 regulation is that it directly affects level of its target transcript Myc and therefore lower E2F concentration. The other one is directly affects E2F inducer Cdk6 and CycE. Besides these two scenarios, miR449 also target Cdc25A, a oncoprotein that positively regulate E2F and is required for progression from G1 to the S phase of the cell cycle. It is competent to activate the G1/S cyclin-dependent kinases Cdk4 and Cdk2 by removing inhibitory phosphate groups from adjacent tyrosine and threonine residues. To some extent because the negative effect of miR449 on Cdc25A can be considered by choosing different inhibitory parameters of miR449 on CycD-Cdk4/6 and CycE-Cdk2 complex, our network does not explicitly include the negative effect of miR449 on Cdc25A.
Mathematical models of cell cycle regulation including Rb-E2F pathway have been previously proposed [25,[39][40][41][42], however, the effects of miRNAs on this network have not been investigated. Meanwhile, miRNAs are presented as a family of critical regulators of almost all cellular processes. So investigation on the regulatory mechanisms of miRNAs is urgent and meaningful especially for therapeutic manipulation in the treatment of cancer. In the present paper, we dedicate to explore the regulatory effects of miR449 on Rb-E2F pathway. According to the scheme of Figure 1, the dynamical relations of the network are characterized by following nonlinear ordinary differential equations (ODEs):    [25], and on the biochemical constraints [21,22,26,42] to give simulations and bifurcation diagrams that are consistent with known dynamical behaviors of miR449 and Rb-E2F pathway. In the following simulations, all the values of employed parameters are shown in Table 1 unless specified elsewhere.

The effects of miR449 on Rb-E2F pathway
When analyzing qualitative information about Rb-E2F pathway, we choose the growth factor S as a governing parameter because it is the key parameter responsible for the G1/S transition in mammalian and plant cell cycle. Moreover, E2F plays a major role during the G1/S transition and CycE-Cdk2 as an output of Rb-E2F pathway is an important mark of cell cycle transition from G 1 to S phase. Thus, in order to investigate the dynamical potential of miR449 in regulation of cell cycle, we will next examine and compare the dynamical behaviors of E2F and CycE-Cdk2 complex without regulation of miR449 or under this regulation.
In the following, if the concentration of E2F lingers on a high stable level state, then the cell cycle is at a state of uncontrolled cell proliferation [43][44][45][46]. When the concentration of E2F is at a state of sustained oscillations with high amplitudes, then the cell cycle is at a state of normal cell cycle process [46]. If the concentration of E2F is at a state of sustained oscillations with low amplitudes or a low stable level state, then the cell cycle is arrested in G1 phase [46]. While the concentration of E2F is at a state of relatively high level (lower than the high stable state and higher than the low stable state), then the cell undergoes apoptosis [46,47]. Note that, in the following bifurcation diagrams, the red and black lines respectively represent stable and unstable steady states, green dots are the maxima and minima of the stable limit cycles, while blue open circles denote the maxima and minima of the unstable limit cycles, and X h is the horizontal coordinate of point X .
Without the effects of miR449, bifurcation diagrams of E2F and CycE-Cdk2 complex concentrations versus the growth factor S are obtained in Figure 2 and 3, where A is a subcritical Hopf bifurcation point and B is a saddle-node bifurcation point. Contrasting Figure 2 with Figure 3, it can be seen that E2F and CycE-Cdk2 are so sensitive to growth factor S that S~B h &1:042 is sufficient to activate them. And the dynamical behaviors of both E2F and CycE-Cdk2 are almost synchronous, that is because CycE is produced by E2F. From Figure 2 and 3, there exists an unique stable steady state corresponding to the quiescent state when S is small enough. As S increases gradually, bistability occurs in the range starting from A h to B h . Particularly, a subcritical Hopf bifurcation appears at A h &0:1468, which leads to a series of sustained oscillations of both E2F and CycE-Cdk2 corresponding to the normal cell cycle. When the growth factor S keeps on going out of the bistable range, these oscillations vanish and then E2F and CycE-Cdk2 step into their high stable steady states. Meanwhile, the value of the high stable state in Figure 2 is almost close to the maximum of the oscillation. It is more obvious especially for CycE-Cdk2 as shown in Figure 3. This means that if S surpasses the bistable range, cell undergoes excessive proliferation which will leads to uncontrolled cell proliferation.
Under the regulation of miR449, bifurcation diagrams for the concentrations of E2F and CycE-Cdk2 complex corresponding to the growth factor S are shown in Figure 4 and 5. To highlight distinct features associated with the one deregulation of miR449, we consider similar parameters to above except k smiR~1 :4 instead of k smiR~0 . Obviously, both A and B are saddle-node bifurcation points; C and D are supercritical Hopf bifurcation points. From Figure 4 and 5, it can be seen that there exists an unique stable steady state corresponding to the quiescent state if S is small enough. As S increases gradually, bistability occurs in the range starting from A h &0:9312 to B h &1:117. A lower and a higher stable steady state coexist. However, when S keeps on going out of the range, the higher stable steady state loses its stability and a supercritical Hopf bifurcation appears at C h &1:119, which leads to a stable branch of limit cycle corresponding to a series of sustained oscillations. When S is increased to D h &3:854, these oscillations vanished. In fact, the amplitudes of these oscillations near C h and D h are low, which indicates the G1/S transition does not happen and corresponds the cell cycle arrest. While at an intermediate k smiR value between C h and D h , one can notice that the amplitudes of these oscillations are high, which may corresponds to the normal cell cycle process. When SwD h , E2F steps into the highest stable state which monotonically increases at the beginning and then tends almost to the horizontal line. While the concentration of CycE-Cdk2 reaches its saturation state more quickly. But it is much lower than the maximum of the oscillations corresponding to the normal cell cycle. This means that the cell undergoes apoptosis instead of proliferation and cell cycle arrest. For thorough testify the dynamical potential of miR449 in regulation of cell cycle, we also compared the time courses of E2F and CycE-Cdk2 complex without regulation of miR449 or under this regulation as show in Figure 6 and 7. Without the regulation of miR449, dynamical processes of E2F and CycE-Cdk2 complex with S~3 and S~5 are show in Figure 6(a) and (b), respectively. It can be seen that the concentration of E2F and CycE-Cdk2 complex keep in a low level at the beginning and suddenly jump to a high level. The only difference is that the jump is occurs at t&85 in Figure 6(a) and t&69 in Figure 6(b). That is, when S~3, the cell transits from a quiescent state to an excessive E2F-induced proliferation state at t&85. And when S~5, this transition happens at t&69. Under the regulation of miR449, the time courses of E2F and CycE-Cdk2 complex with S~3 and S~5 are show in Figure 7(a) and (b), respectively. The initial state of E2F is an inappropriate activated state with [E2F] = 1.2mM. It was easy to see that there are essential differences between the two diagrams. The concentration of E2F is driven from a high level state to a sustained oscillatory state as shown in Figure 7(a). Specially, the concentration of E2F quickly drops at the beginning due to strong up-regulation of E2F on miR449 and heavy downregulation of miR449 on E2F. And then the concentrations of E2F, CycE-Cdk2 and miR449 undergo sustained oscillations with low amplitudes in the primary stage and high amplitudes later due to the effects of the negative feedback loops between E2F and miR449. This means that the cell cycle is quickly driven from an excessive E2F-induced proliferation state to a cell cycle arrest state, which rendered the cell to repair DNA damage. After all the damage is repaired, the concentrations of E2F, CycE-Cdk2 and miR449 undergo sustained oscillations with high amplitudes corresponding to a normal proliferation state. In Figure 7(b), the concentration of E2F are driven from a high level state to a relatively high level state compared with the low level in Figure 6. This means that miR449 has the ability to drive cell cycle from excessive E2F-induced proliferation state to a state of apoptosis. Therefore, miR449 plays a dual role in regulation of cell cycle. Specially, miR449 drives the cell cycle from an excessive E2Finduced proliferation state returning to a normal proliferation state in a moderate stimulation of S by promoting cell cycle arrest, but leads to apoptosis for a high stimulation of S.

The effects of miR449 on E2F activity
In order to further explain the effectiveness of miR449 in inhibiting E2F activity, we consider the effects of variations in the parameter k smiR (a rate constant of miR449 production) and S on the concentration of E2F. The results in Figure 8 show that the concentration of E2F is sensitive to the variation of k smiR and not particularly sensitive to the variation of S. Therefore, E2F activity can be significantly affected by the regulation of miR449.
As shown in Figure 8, the concentration of E2F decreases with k smiR because a larger rate of miR449 production means lower E2F activity. The stability of the equilibrium can be changed with the variation in k smiR . In absence of miR449 regulation, i.e., k smiR~0 , the unique equilibrium is stable. E2F settles on a high stable steady state corresponding to an excessive E2F-induced proliferation state. When miR449 is introduced, i.e., k smiR w0, E2F drops rapidly with increasing k smiR . It means that the cell cycle is quickly driven from an excessive E2F-induced proliferation state to a cell cycle arrest state (because the concentration of E2F is   lower than the maximum of the oscillations corresponding to the normal cell cycle in the left-neighborhood of point A). As k smiR increases gradually, the stability of the equilibrium transforms from stable to unstable, meanwhile sustained oscillations appear. In fact, one can notice that the amplitudes of these oscillations near A and B are low, which indicates that the G1/S transition does not happen and corresponds to the cell cycle arrest. However, the amplitudes of these oscillations near a certain value of k smiR are high (specially, k smiR &1:053 in Figure 8(a), k smiR &1:402 in Figure 8(b), k smiR &1:533 in Figure 8(c)), which may corresponds to the normal cell cycle process. When the parameter k smiR keeps on going out of the oscillatory range, oscillations vanish and the unique equilibrium regains its stability. E2F steps into a low stable steady state corresponding to the state of cell cycle arrest. If the value of k smiR is large enough, the activity of E2F is almost completely inhibited and cell growth is permanently arrested in the G1 phase. However, it is worth noting that this inhibitory process of miR449 on E2F has a little difference for different level of S in the middle range of k smiR . For example, in Figure 8(a) with S~1:2, there exist a subcritical Hopf bifurcation A at k smiR &1:017 and a supercritical Hopf bifurcation B at k smiR &1:464. Stable oscillation and equilibrium may coexist for values of k smiR in the interval enclosed by A and a fold limit cycle bifurcation (k smiR &1:015) due to the occurrence of the subcritical Hopf bifurcation. In Figure 8(b) with S~3, there exist two supercritical Hopf bifurcations A at k smiR &1:339 and B at k smiR &2:089 as k smiR increases. The first one is supercritical, resulting in a stable branch of limit cycles. While the second one is also supercritical because the equilibrium branch loses the stability going left and the periodic orbit branch goes left too. Similar to Figure 8  Comparing the three diagrams, one can notice that the larger S is, the larger k smiR and the longer oscillation range will be needed to completely inhibit E2F activity. In other words, the more growth factor, the more miR449 is needed to inhibit the activity of E2F.
Therefore, no matter how large S is, a value of k smiR , which always exists, is sufficient to inhibit the concentration of E2F from a high stable steady state to a state of sustained oscillations whose amplitudes undergo a change from low to high and then again to low. In other words, no matter how large S is, there always exists a value of k smiR that is sufficient to drive the cell cycle from an excessive E2F-induced proliferation state (high stable steady state with high concentration) to a normal proliferation state (sustained oscillations with high amplitudes) by promoting the cell cycle arrest (oscillations with low amplitudes or high stable steady state with low concentration in the left-neighborhood of point A). For example, assuming that S~3 and the initial state of E2F is an inappropriate activated state with [E2F] = 1.2mM, there exists k smiR~1 :5 that is sufficient to drive E2F from a high stable steady state to a sustained oscillatory state as shown in Figure 9(a). It can be seen that the concentration of E2F quickly drops at the beginning due to strong positive regulation of E2F on miR449 and strong inhibitory regulation of miR449 on E2F. And then the concentrations of E2F, CycE-Cdk2 and miR449 undergo sustained oscillations with low amplitudes in the primary stage and high amplitudes later due to the effects of the negative feedback loops between E2F and miR449. This means that the cell cycle is quickly driven from an excessive E2F-induced proliferation state to a cell cycle arrest state, which rendered the cell to repairing DNA damage. After all the damage is repaired, the concentrations of E2F, CycE-Cdk2 and miR449 undergo sustained oscillations with high amplitudes corresponding to a normal proliferation state. Moreover, no matter how large S is, a value of k smiR , which always exists, is sufficient to inhibit the concentration of E2F from a high stable steady state to a low stable steady state even to zero. In other words, no matter how large S is, there always exists a value of k smiR that is sufficient to drive the cell cycle from excessive E2F-induced proliferation state to a state of cell cycle permanent arrested in the G1 phase (senescence). For example, there exists k smiR~2 :5 that is sufficient to drive cell cycle from excessive E2Finduced proliferation state to a state of cell cycle permanent arrested in the G1 phase as shown in Figure 9(b). Therefore, E2F is so sensitive to the rate constant of miR449 production and miR449 remarkably inhibits E2F activity. These numerical results and analysis are coherent with the experiments. Thus, miR449 acts as an effective tumor suppressor.

Different repressions of miR449 on E2F via three different targets
In order to further compare the repressions of miR449 on E2F by targeting three different transcripts named Myc, CycD-Cdk4/6 and CycE-Cdk2, the dynamical behaviors of E2F versus the inhibition rate of miR449 on these three transcripts are investigated numerically in Figure 10. Particularly, during studying the negative regulatory effect of miR449 on E2F by targeting Myc, we assume that k dCD2 and k dCE2 , the inhibition rates of miR449 on CycD-Cdk4/6 and CycE-Cdk2, are zero. Then a bifurcation diagram in Figure 10(a) shows the function of E2F concentration corresponding to k dM2 , the inhibition rate of miR449 on Myc. Using the same idea, other two bifurcation diagrams respectively describing the dynamical behaviors of E2F level as a function on k dCD2 and k dCE2 are shown in Figure 10 Interestingly, these three bifurcation diagrams are quite distinct. Firstly, Figure 10(a) shows that E2F always settles on a high stable steady state without the regulation of miR449 (k dM2~0 ). When the inhibitory effect of miR449 on Myc is introduced (k dM2 w0), the high stable steady state of E2F decays rapidly with increasing k dM2 . If the value of k dM2 is limited between two supercritical Hopf bifurcation values A h and B h , the high stable steady state of E2F losses its stability and sustained oscillations of E2F appears.  When k dM2 keeps on going out of the range, these oscillations vanish and then E2F steps into its low stable steady states. In Figure 10(c), under the regulation of miR449 on E2F by targeting CycE-Cdk2, the concentration of E2F always stays at a stable state, which descends quickly and then tends to a horizontal level. However, in Figure 10(b), the high stable level of E2F decreases gradually at the beginning of introducing repression of miR449 on E2F by targeting CycD-Cdk4/6. As the k dCD2 increases, E2F steps into a bistable range between two subcritical Hopf bifurcation values A h and B h . Different from the case of monostability in Figure 10(a) and (c), where a lower inhibition rate corresponds to a higher E2F concentration and a higher inhibition rate corresponds to a lower protein concentration due to the negative regulation mediated by miR449, when the inhibition rate lies between the two subcritical bifurcation points A and B, despite the negative regulation mediated by the miR449, a higher inhibition rate may correspond to a higher E2F concentration or a lower inhibition rate may also correspond to a lower E2F concentration, depending on the initial conditions. Thus, the numerical results and analysis suggest that the dynamical behaviors of the negative effects of miR449 on E2F by targeting three different transcript largely differ in the onset of oscillations, bistability and monostability, which needs to be verified experimentally. It is expected that the difference will generate a detailed and precise insight of miR449 mediated Rb-E2F pathway.

Discussion
MiRNAs have been shown to function as integral components of a wide range of cellular processes including cellular growth, differentiation, and disease [11][12][13][14][15]. As such, miRNA dysregulation can have a profound effect on cancer development [21][22][23][24]. Previous studies have shown that miR449 is down-regulated or lost in gastric cells, testicular cancer cells as well as in lung adenocarcinoma cell line and possesses potential tumor suppressor function [21][22][23][24]. In this paper, to investigate the dynamical potential of miR449 in regulation of cell cycle, a mathematical model for Rb-E2F pathway mediated by miR449 are constructed   based in part on the model proposed in Yao-Lee et al. (2008) and numerical simulation and dynamical analysis are performed. By comparing the dynamical behaviors of E2F pathway deregulation and regulation of miR449, we theoretically verified that the miR449 plays a critical role in regulating the cell cycle progression and provides a twofold safety mechanism to avoid excessive E2Finduced unrestricted proliferation by cell cycle arrest and apoptosis. Moreover, numerical simulation suggests that the mechanisms of the negative regulation of miR449 on three different transcripts are quite distinctive which needs to be verified experimentally.
The results in this paper turned out that the model without miR449 has several disadvantages than the one with miR449: N Model without miR449 is risky. As S increases, E2F concentration quickly transit to a high level that will induce unrestricted proliferation. However, the model with miR449 is safer. It is the effective inhibition of E2F by miR449 that can prevent this abnormal proliferation by inducing cell cycle arrest or irreversible apoptosis.
N Model without miR449 has little practical significance. The periodic orbit corresponding to the normal cell proliferation is unstable in Figure 2 and 3, that is, the orbits can scarcely be reached in real environment. But the model with miR449 is proper in Figure 4 and 5 because the limit cycle corresponds to the cell proliferation are stable.
N Model without miR449 is less robust. It can be seen that the oscillatory range of the model without miR449 is smaller, indicating that the model under regulation of miR449 is more robust. Thus, miR449 has the ability to enhance the fidelity, robustness and flexibility in temporal regulation.
The study of Rb-E2F pathway mediated by miR449 may help us to analyze the whole cell cycle process mediated by other miRNAs more easily and may even provide an new insight for therapeutic manipulation in the treatment of cancer. A better knowledge of the miRNAs mediated gene regulation is also of interest for the bio-engineering or artificial control of specified components, interactions, and even network functions [1,2,[21][22][23][24].
Certainly, besides Rb-E2F pathway mediated by miR449, more complex models will be needed soon to accommodate the complexity of miRNAs in gene regulatory networks. For example, further improvements can be done to the present model in future. Firstly, for simplicity, we only considered the cell cycle of G 1 to S phase in this study. Further exploration on the whole cell cycle process mediated by miRNAs still needed be researched intensively, which will provide a more comprehensive view on how cell cycle are regulated by miRNAs. Secondly, how do miR449 and miR34 cooperatively fine tune the E2F and p53 pathways in the case of DNA damage is another interesting and challenging problem [21]. After DNA is damaged, p53 and E2F coordinate cell cycle progression and cell fate decision including growth arrest, DNA repair and apoptosis [47]. Moreover, p53 induces miR34 and E2F increases the levels of miR449 [21,22,26]. The most exciting aspect is miR449 and miR34 each repress E2F, but promote p53 activity, allowing efficient cross-talk between two major DNA damage-responsive gene regulators, which may thus represent a determinant for therapeutic manipulation in the treatment of cancer by irradiation or chemotherapy [21]. Therefore, to obtain a comprehensive and coherent picture of the regulatory roles of miRNAs are urgent, and such important aspects still need to be considered.

Materials and Methods
Numerical simulation and bifurcation analysis of the ODEs in our paper were carried out with XPPAUT, a software program freely downloadable from http://www.math.pitt.edu/*bard/ xpp/xpp.html. XPPAUT can automatic detect bifurcations of fixed points and limit cycles. It is base upon the following strategy: N Use Sing pts Go command to compute stable steady point for the system. XPPAUT will used this point as initial conditions for drawing bifurcation diagrams.
N Use the Axes command to tell XPPAUT what parameters are varied and what is to be plotted, as well as the range of the graphs.
N Use the Numerics command to define all the AUTO numerical parameters such as the direction, step size, and so on.
N Use the Run command to run the bifurcation. Grab the Hopf bifurcation and track the periodic orbits.
Values for the rate constants in our model are shown in Table 1 unless specified elsewhere. They were chosen based on the literature whenever possible [25], and on the biochemical constraints [21,22,26,42] to give simulations and bifurcation diagrams that are consistent with known behaviors of miR449 and Rb-E2F pathway. XPPAUT codes including all the parameter sets are available in the Supporting Information Section.

Supporting Information
File S1 In order to calculate the bifurcation diagram of [E2F] and [CycE-Cdk2] with S as a control parameter at k smiR~0 , we provide XPPAUT code of Figure 2  File S8 In order to obtain the bifurcation diagram of [E2F] with k smiR as a control parameter at S~3, in this Supporting Information File S8 we provide XPPAUT code of Figure 8 File S12 In order to obtain the bifurcation diagram of [E2F] with k dM2 as a control parameter at k dCD2~0 and k dCE2~0 , we provide XPPAUT code of Figure 10(a) here.

(ODE)
File S13 To plot the bifurcation diagram of [E2F] with k dCD2 as a control parameter at k dM2~0 and k dCE2~0 , XPPAUT code of Figure 10(b) is provided in this Supporting Information File S13.

(ODE)
File S14 To calculate the bifurcation diagram of [E2F] with k dCE2 as a control parameter at k dM2~0 and k dCD2~0 , in this Supporting Information File S14 we include XPPAUT code of Figure 10