 Research
 Open Access
 Published:
Dynamics of miRNA driven feedforward loop depends upon miRNA action mechanisms
BMC Genomics volume 15, Article number: S9 (2014)
Abstract
Background
We perform the theoretical analysis of a gene network subsystem, composed of a feedforward loop, in which the upstream transcription factor regulates the target gene via two parallel pathways: directly, and via interaction with miRNA.
Results
As the molecular mechanisms of miRNA action are not clear so far, we elaborate three mathematical models, in which miRNA either represses translation of its target or promotes target mRNA degradation, or is not reused, but degrades along with target mRNA. We examine the feedforward loop dynamics quantitatively at the whole time interval of cell cycle. We rigorously proof the uniqueness of solutions to the models and obtain the exact solutions in one of them analytically.
Conclusions
We have shown that different mechanisms of miRNA action lead to a variety of types of dynamical behavior of feedforward loops. In particular, we found that the ability of feedforward loop to dampen fluctuations introduced by transcription factor is the model and parameter dependent feature. We also discuss how our results could help a biologist to infer the mechanism of miRNA action.
Background
An important role in regulation of gene expression in higher eukaryotes, plants and animals belongs to miRNAs, the endogeneous small noncoding RNAs that bind to partially complementary sequences in target mRNAs. The miRNAs are involved in regulation of development, differentiation, apoptosis, cell proliferation [1–4], as well as in the progression of numerous human diseases, such as chronic lymphocytic leukemia, fragile × syndrome, and various tumor types [4–8]. Each miRNA molecule may target hundreds of mRNAs, and, vice versa, some targets are combinatorially affected by multiple miRNAs [9–11]. Comparative phylogenetic studies uncovered the conserved miRNAbinding sequences in more than one third of all genes, that lead to a suggestion that the miRNA regulation may be relevant to a large portion of cellular processes [12, 13].
However, it should be noted that the mechanisms of miRNA action are not unanimous. There are experimental evidences that miRNAs regulate gene expression through translational repression, mRNA deadenylation and decay, however the contribution and timing of these effects remain unclear [14, 15]. Although some studies show translational repression without mRNA decay [16], others point to decay as a primary effect [17, 18]. It was also demonstrated in several papers that the destabilization of target mRNA is accompanied by degradation of a miRNA molecule [18, 19].
Mathematical modeling is often used to discriminate between these tentative mechanisms. It was shown that variability in the miRNA action on protein production in different experiments could be due to differences in ratelimiting steps [20]. Morozova et al [21] developed a model, based on coexistence of all proposed mechanisms of microRNA action, and the apparent mechanism detected in an experiment is determined by the relative values of the intrinsic characteristics of the target mRNAs and associated biological processes. It was shown in the detailed review [22] that the analysis of the transient protein translation dynamics provides enough information to verify or reject a hypothesis about a particular molecular mechanism of microRNA action on protein translation. Besides both in vitro and in vivo assays with Caenorhabditis elegans showed that targets can also protect miRNAs from active degradation [18, 23]
The role of miRNA in gene regulatory networks becomes currently a subject of wide speculations. In general, it is expected that motifs with miRNA can buffer the consequences of noise action in gene expression in order to confer a robustness to environmental perturbation and genetic variation [12, 24]. The anticorrelative expression of miRNAs and their target mRNAs was documented in many cases [12, 25], where it points out that while transcription primarily controls target gene expression, the miRNAs lend further reinforcement to gene regulation by attenuation any unwanted transcripts. The second class of genetic buffering by miRNA emerges in cases, where both the miRNA and the target are coexpressed at intermediate levels [26–28]. It was proposed that intermediate miRNA quantities and low target avidity (the result of a single "seed" binding site) intentionally provide the target protein synthesis, but place a burden on this process, which, most likely, provides an approach to buffer stochastic fluctuations in the mRNA level [12]. However, it should be noted that up to date there is the only direct observation that miRNA buffers gene expression against perturbation [29].
The transcription networks contain several biochemical wiring patterns called network motifs, and one of the most significant recurring motifs is the feedforward loop (FFL) [30–32], in which the upstream transcription factor (TF) regulates the target gene via two parallel pathways: directly, and by interaction with a second molecule, which also regulates the target gene. An assigned value for a pathway is defined as positive, if the total number of negative interactions in the pathway is even, and negative otherwise. The FFL is named as "coherent" if the sign of direct regulation path coincides with the overall sign of the indirect regulation path, and "incoherent" otherwise.
Recent computational analysis demonstrated that FFL with TF and miRNAs are overrepresented in gene regulatory networks, assuming that they confer useful regulatory opportunities [33]. There are two possible structure configurations in each coherent and incoherent FFL, containing miRNA, as shown in Figure 1. We will specify these configurations as the type 1 or 2 coherent FFLs, and the type 1 or 2 incoherent FFLs, respectively.
A simple mathematical modelling was used recently to explore the capability of FFLs with miRNA to buffer fluctuations in gene expression. Osella et al. [34] introduced the model, describing a target gene regulation in the type 1 incoherent FFL, however, solved analytically the coupled algebraic equations only, obtained by trivial reduction of the coupled 1st order ordinary differential equations (O.D.E.) to a steady state. For stochastic statement of the problem in [34] O.D.E. were solved numerically. It was shown that with respect to the simple gene activation by TF, the introduction of the miRNAmediated repressing pathway can significantly dampen fluctuations in the target gene output for essentially all the choices of input parameters and initial conditions. Moreover, this noise buffering function was expected to be direct consequence of the peculiar topology of the FFL.
There is a critical imperfection in any mathematical model based on the steady state data analysis. Firstly, latter measurements [35, 36] of protein abundance and turnover by parallel metabolic pulse labelling for more than 5,000 genes in the NIH 3T3 mouse fibroblasts showed that the half live times of many proteins in these cells are longer than the cell cycle duration of 27.5 hours. For proteins, which half life time is comparable or longer than the cell cycle, the genuine steady state is not reached before the cell division. As a consequence, quantities per cell are not only defined by synthesis and degradation, but also by an initial number of protein molecules at the beginning of cell cycle. Secondly, an absence of time derivatives (as in the steady state approximation) is provided by simple replacement of the initial coupled kinetic nonlinear equations with corresponding algebraic equations, avoiding any information about the temporal behaviour of cell components. A variety of the cell process models usually demonstrate versatile dynamics at early stages even in case of coincidence of stationary values. For these reasons it is much more faithful and informative, however difficult, to study coupled ODEs at every time moment, aiming to decipher mechanisms underlying its dynamics. It is worth to mention the importance of analytical solutions, which provide a genuine 'check point' for any numerical simulation often obtained by means of different methods.
In general, it is not simple to derive (either partial or ordinary) differential equations, describing a biological system and obtain an analytic solution to them. Fortunately, many biological systems demonstrate the module structure, and therefore, at least, in theory, they can be decomposed in segregated subsystems, and, with luck, the equations can be solved in closed form.
Another challenge in mathematical modeling of biological systems consists in the nonuniqueness of solutions. This nonuniqueness is partly caused by imperfections of current experimental methods, which are often unable to provide the measurements of all the parameters required to describe the system dynamics. As a result the parameter values have usually to be found by numerical fitting of the solutions to data. This task, called "an inverse problem", is, in a way, tricky, because there are usually several parameter sets in good agreement with data. It means that within a given accuracy there are several models, which are, in a sense, equivalent in the data description. Nonuniqueness provided by numerical simulations represents a crucial problem for contemporary mathematical modelling in biology, and usually an appropriate solution may be found by invoking additional information from complicated experiments. On the other side, an intrinsic nonuniqueness may have transparent biological reasons, as biological subsystems are robust and can tolerate large parameter variation, as long as the core network topology is retained [37, 38]. Therefore, mathematical treatment and refinements are still required to better understand how the component interactions result in the system complex behaviour.
Here we consider a gene network subsystem composed of a FFL mediated by TF and miRNA. We perform the analysis of three mathematical models, which describe temporal dynamics of gene expression in FFL under an assumption of different mechanisms of miRNA action. In contrast to [34] we examine the FFL dynamics quantitatively at the whole time interval of cell cycle. We rigorously proof the uniqueness of solutions to these models and obtain the exact solutions in one of them analytically. We show how different mechanisms of miRNA action lead to distinctive dynamical behaviour of FFLs. In particular, we find that the ability of FFLs to dampen fluctuations introduced by TF is both the model and the parameter dependent feature. We also discuss how our results could help a biologist to infer the mechanism of miRNA action.
New computational experiments [21, 22] justified a hypothesis about coexistence of distinct miRNAmediated mechanisms of translation repression, however it is still open for modelling based on differential equation analysis. The mathematical approach developed here seems to be complimentary to analysis performed in [22] and useful in those problems, however, it will require future work.
Results and discussion
Mathematical models of FFLs
The biological system under study will be described with 5 variables, representing the number w of mRNA molecules transcribed from the TF gene, the number q of TF molecules, the number s of miRNA molecules, the number r of mRNA molecules transcribed from the target gene, and the number p of target protein molecules.
We deal with the models proposed in [34] and consider 4 conventional configurations of FFLs with miRNA, see Figure 1. For each gene in FFL we consider transcription, translation, degradation and interactions between genes in the regulatory network.
As molecular mechanisms of miRNA action are not clear so far we consider three different models:

the model, in which miRNA represses translation of its target  Stop model,

the model in which miRNA promotes target mRNA degradation  Target degradation model and

the model in which miRNA is not reused but degrades along with target mRNA  Dual degradation model.
and the 'working titles' for the models introduced above will be used for brevity.
Let us write five coupled differential equations proposed in [34] for a feedforward loop (FFL)
Here k_{ w } and k_{ q } are rates of TF mRNA and TF synthesis, k_{ s }(q) and k_{ r } (q) are rates of transcription of the regulated gene, k_{ p }(s) is the rate of target protein synthesis; g_{ w }, g_{ q }, g_{ s }, g_{ s } and g_{ p } represent the degradation rates of the corresponding species.
Before we shall specify the types of production functions in equations (14) let us remind that various dynamic processes in a complex system can be described as a progression from an initial quantity that accelerates (or decelerates in case of repression) and approaches a plateau over time. When a detailed description is lacking, a sigmoid function is used, that is based on an idea by A.V. Hill (1910) to describe the sigmoidal oxygen binding curve of haemoglobin. The general form of the Hill function is written as:
where both α, β are arbitrary given parameters of the process. The rate c represents either production (c ≥ 2) of a quantity x or repression of it (c ≤ −2) in time. Consequently, the Hill function will be represented by either uprising or falling down graph, and its slope depends on the value of c. The easiest way to get a nontrivial regulation type is to prescribe the rate values c = 2 and c = −2 for activation and repression, respectively, however, even in this case the mathematical problem becomes the nonlinear one.
In the problem considered the production functions k_{ s }(q) and k_{ r }(q) are assumed to be the classical Hill functions in the form k(q) = (k_{ max }q^{c})/(h^{c} + q^{c}). The parameters h_{ s } and h_{ r } specify the amount of TF, at which the transcription rate of the miRNA or target gene reaches one half of its maximal value (k_{ s } or k_{ r }), and c is the Hill coefficient, representing the steepness of the regulation curve, see (2).
Therefore for the type 1 coherent FFL (see Figure 1) the Hill functions will have the following form:
while for type 1 incoherent FFL the production function for target mRNA will be different:
The difference between two types of each of coherent or incoherent loops lies in form of the production functions for both target mRNA and miRNA.
We shall consider three different mechanisms of the miRNA action, following [34], and to model the effect of direct translational repression of target mRNA we consider the translation rate of the target k_{ p }(s) to be a repressive Hill function of the number of miRNA molecules: ${k}_{p}\left(s\right)=\left({k}_{p}{s}^{c}\right)/\left({h}^{c}+{s}^{c}\right)$. The parameter h specifies the quantity of miRNAs, at which the translation rate reaches one half of its maximal value k_{ p }, and c = −2 is again the Hill coefficient.
To model miRNA action in the destabilization of target mRNA we add to the basal rate of mRNA degradation g_{ r } (in absence of miRNAs), a term, which represents an increasing Hill function of a copy number of miRNAs, where g_{ max } is the maximal value of the degradation rate in case of high miRNA concentration, h_{ deg } is the dissociation constant of miRNAmRNA interaction, and c = 2 is the Hill coefficient.
As in [34] we discuss also a tentative destabilization of target mRNA, accompanied by the miRNA degradation. The miRNA forms a complex with its target, which degrades with it, instead of being reused. This complex degradation constant will be denoted as k_{ rs }, and the resulting nonlinear (due to the multiplicative term rs) equations will have a form:
Details of mathematical analysis of the coupled ODE are given in Additional file 1. For the Stop model we obtained the exact solutions to the nonlinear ODE, which coefficients explicitly depend on initial values. Main advantage of exact solutions to corresponding nonlinear differential equations is that they do not require any parameter fitting and provide also a reliable basis for verification of numerical results. For the Target and Dual degradation models the numerical solutions to the problem are obtained. For each model considered we proved also the uniqueness of solutions, i.e., the onetoone correspondence between given parameters and solutions.
Comparative analysis of FFL temporal behaviour under different models
The mechanisms underlying miRNAmediated repression are not clear so far, and for this reason we consider three models of the miRNA action described above.
Four different topologies of FFLs mediated by TF and miRNA are possible in theory, as in Figure 1. We shall analyse the dynamical behaviour of all these network topologies in frame of the models described above, that leads, in total, to consideration of 12 different variants of regulation in FFLs.
We shall use for brevity the following abbreviations for the FFL identification: 1C will mean the type 1 coherent loop, 1In  type 1 incoherent loop, 2C  type 2 coherent loop and 2In  type 2 incoherent loop, resp. We assume that the initial number of molecules in a loop is equal to one half of that obtained just before cell division, which often (however, not necessarily) corresponds to the steady state level. The results for the Stop model are based on exact and explicit solutions obtained to the nonlinear O.D.E. 1 (see Additional file 1), and, therefore are free of any data fitting procedure and numerical approximations. For this reason exact solutions may provide also a reliable base and "check points" for numerical solutions in those (Target and Dual degradation) models, where an exact solution is hardly obtainable. We begin with brief description of the temporal variation of the molecule quantity of each player and in each of three models. To simplify comparison we shall use one and the same parameter set for all FFL and models.
Stop model
In the Stop model the behaviour of target mRNA and protein is different in all FFLs because miRNA stops translation of the former and does not promote its degradation. Both 1C and 2In loops form the identical bellshaped target mRNA profiles due to repression by TF, in both 1In and 2C loop the target mRNA profiles are also identical and increase in time to a steady state value (Figure 2). In both 1In and 1C loop, in which TF activates miRNA gene, the target protein shows pulselike behaviour due to repression mediated by miRNA (Figure 2A,B,E,F). However at steady state in the 1C loop the number of target protein molecules is much lower than in the 1In loop.
In 2C loop the number of target protein molecules firstly slightly rises, than falls down due to miRNA action), but with continuous repression of the miRNA synthesis, the target protein begins to rise up to the steady state level (Figure 2H,D).
In 2In loop the target protein profile also exhibits the wavelike behaviour (Figure 2C,G). Firstly, it rises to almost stationary level, after that it slightly decreases due to repression of mRNA synthesis by TF and repression of mRNA translation by miRNA, and eventually grows again to the stationary level due to decrease of the miRNA molecules number.
Target degradation model
Contrary to the Stop model in the Target degradation model the temporal dynamics of both target mRNA and protein is similar. In this model the target protein production is a linear function of the target mRNA molecules number, and miRNA promotes the degradation of target mRNA. In 1C loop two mechanisms are active, namely, mRNA degradation under miRNA action, and the repression of target mRNA synthesis. In 2In loop TF represses miRNA production. This explains the difference in dynamics of target RNA and protein in these loops (Figure 3B,F,C,G). In both loops the dynamics shows pulselike behaviour, however in the 1C loop the numbers of target mRNA and protein molecules are smaller (about 1300 vs. 1500 molecules) at steady state, and the steady state level is reached later (about 4000 sec vs. 3000 sec) in comparison with 2In loop.
The dynamics of target mRNA and protein in both the 1In and the 2C loops has a form of increasing function, tending to a constant value (Figure 3A,D,E,H). The loops behaviour differs in time, when the steady state levels of the target mRNA and protein are reached, as well as in numbers of target mRNA and protein molecules at steady state. In 1In loop these numbers are smaller (about 3000 molecules vs. 3500 molecules) than in 2C loop, that can be explained by promotion of the target mRNA degradation by miRNA in the 1In loop.
Dual degradation model
In this model both miRNA and target mRNA degrade due to the duplex formation, and the miRNA molecules are not reused. As in the Target degradation model both target mRNA and protein show similar dynamical behaviour. The dynamics of target mRNA and target protein both in 1C and 2In loops shows pulselike behaviour (shown in Figure 1 of Additional file 2), similar to that observed in these loops in the Target degradation model (as shown in Figure 3B,C,F,G). In this model in 21n loop the steady state level is approached earlier (3000 sec vs 4000 sec), than in 1C loop. The target protein molecules number at steady state in 21n loop is also higher than in 1C loop (1800 molecules vs 1400 molecules). Moreover, all these numbers in the Dual degradation model are higher than in Target degradation model (Figure 1 of Additional file 2).
In both 1In and 2C loops the dynamics of target mRNA and protein takes a form of increasing function, tending to constant value (Figure 1 of Additional file 2). As in Target degradation model in 1In loop the steady state is reached earlier (2500 sec vs. 4000 sec) and the level of target protein at steady state is lower (3000 molecules vs. 4000 molecules) than in 2C loop. Again, in the Dual degradation model these numbers are bigger than in the Target degradation model.
Dynamical behaviour of FFLs
Next step is to describe the behaviour of each type of FFL in the framework of each model under consideration. We should demonstrate that FFLs mediated by miRNA and TF may have many possible outcomes, depending on the nature of the relationships between the loop elements.
The study of dynamical behaviour of FFLs is based on the variation of initial conditions, and model coefficients with subsequent analysis, how these changes affect the target protein molecules number. To simplify comparison in each numerical experiment we used one and the same parameter set, described in the Section above, except of the coefficient value, which effect on the loop behaviour will be analyzed. The dynamical analysis is important not only for the noise buffering ability, but also for demonstration of how sensitive are the control pathways to parameters and initial state.
The efficiency in control of the target protein synthesis in all three models depends upon the quantity of TFs (which, in turn, is a function of k_{ w } and k_{ q }), the number of miRNA copies (i.e., the function of k_{ s } and h_{ s } defining the affinity of TF to the promoter of miRNA gene), as well as the strength of the miRNA action. In the Stop model the strength of repression of mRNA translation by miRNA is defined as 1/h_{ p }. In the Target degradation model the degradation of the target mRNA is described by a term, which represents an increasing Hill function of a copy number of miRNAs, while 1/h_{ g } represents the strength of miRNA action. In the Dual degradation model the degradation constant of the mRNAmiRNA complex k_{ rs } is introduced. Therefore, in short, the idea of an analysis given below is to fix a type of FFL and investigate how the target protein molecule number p will vary in all models considered.
However, the detailed description of each FFL behaviour under variations of every parameter is quite lengthy, and for lack of space we shall do it for the type 1 incoherent (1In) loop only. For other FFL subjected to parameter variation we have to refer to corresponding Figures and captions; the descriptive analysis will be similar to the one made for 1In loop.
Type 1 incoherent (1In) loop
This FFL is characterized by direct activation and indirect repression pathways, in which TF acts on target protein production.
Variation of synthesis and degradation parameters
Mathematical analysis showed that the increase of h_{ p } in the Stop model and h_{ g } coefficient in the Target degradation model results in the increase of the target protein quantity (Figure 2 of Additional file 2). The difference in target protein production is the largest for big values of these coefficients. The increase of the k_{ rs } coefficient in the Dual degradation model results in the fall of the target protein quantity. Noteworthy, at early times (up to about 1000 seconds) the variation of h_{ g } or k_{ rs } coefficients has little influence on target protein production. In all models the difference in target protein production at different values of coefficients increases with time, as shown in Figure 2 of Additional file 2.
In all models the increase of h_{ s } leads to increase in the target protein quantity (Figure 4). In both Target and Dual degradation models at early times (up to ca. 1000 seconds) the target protein production does not depend on h_{ s } variation; at later times the difference between two bounding h_{ s } values in the Target degradation model is much smaller than in other models. In all models the largest difference in target protein production is observed for large h_{ s } values (Figure 4).
In all models the increase of kr results in the increase of the target protein production, the effect of the k_{ r } variation being larger at a later time (Figure 2 of Additional file 2). In all models the profiles nearly coincide up to a moment ca. 250 seconds and diverge afterwards. Later in the Stop model the profiles form a peak, which amplitude rises as k_{ r } increases. Both the Target and Dual degradation models exhibit similar behaviour: the profiles firstly tend to a minimum and increase to a constant value afterwards. In both models the steady state is achieved faster (around 2000 seconds) than in the Stop model (around 5000 seconds) (Figure 2 of Additional file 2).
In both the Target degradation and Dual degradation models applied to the 1In loop the increase of the TF translation rate k_{ q } leads to increase of target protein quantity, the difference in target protein being the largest for small values of the coefficient (Figure 4). In Dual degradation model the steady state is reached earlier that in Target degradation model. In Stop model for this loop the largest number of target protein molecules is observed at intermediate values of the kq coefficient (Figure 4).
Variations in initial data
In all models the variation of the initial number of TF and miRNA molecules does not change the number of target protein molecules at steady state. In both Stop and Dual degradation model the form of target protein profile changes from the bellshaped to the Ushaped one as the initial number of miRNA molecules rises (Figure 4G, I). In Target degradation model all profiles have a form of increasing curve tending to a steady state and show moderate dependence on the change of the initial number of miRNA molecules (Figure 4H). The increase of the initial number of TF molecules results in the change of the target protein profile from the Ushaped form to the bellshaped one in both the Target degradation and the Dual degradation models (Figure 2H, I of Additional file 2). In the Stop model all the target protein profiles have the bellshaped forms, and their amplitudes decrease with the TF molecules number growth (Figure 2G of Additional file 2).
Remarks concerning the analysis for Type 1 coherent (1C), Type 2 incoherent (2In) and Type 2 coherent (2C) loops
1C FFL is characterized by a synergetic action via both the direct and indirect pathways;
2C FFL is characterized by coherent activation of a target via direct and indirect pathways;
2In FFL describes an indirect pathway of activation and direct pathway of the target repression; all of them are schematically presented in Figure 1.
We omitted for brevity a detailed analysis of these FFLs response on coefficient variation and refer to the Figures, containing graphs for every FFL. There are the Figure 5 and Fig. A3 in Additional file 2 for 1C loop; the Figure 6 and Fig. A4 in Additional file 2 for 2C loop: the Figs. A5 and A6 in Additional file 2 for 2In loop, respectively.
In short, our results show that FFLs mediated by miRNA may have many possible outcomes, depending on interaction between the loop elements. The target protein profiles can take different forms, which are unambiguously defined by initial conditions and model coefficients. In most cases the variation of model coefficients leads to results, which could be intuitively explained by consideration of the miRNA action in a model and the topology of a loop, however, in several cases the response of the system is hardly predictable. This especially concerns the variation of k_{ q }, the parameter, which defines the quantity of TF. In the Stop model in 1C loop the TF represses and in 2C loop it activates the target protein synthesis in parallel via direct or indirect pathways, respectively. Therefore in 2C loop the rise of k_{ q } leads to an increase in target protein quantity, while in 1C loop the relation is opposite (Figure 5 and Figure 4 of Additional file 2). However, in the Stop model applied to both 1In and 2In loops the maximal number of the target protein molecules is observed at intermediate k_{ q } values, moreover, in the last one the effect becomes visible only after 750 seconds (Figure 4 andFigure 6 of Additional File 2). In general, any noticeable variation in molecule quantity provided by an intermediate value of a parameter can be explained by simple mathematical analysis of the Hill function sigmoid, see Additional file 1.
Noise buffering by miRNA
The comparative analysis performed above shows significantly diverse reaction of each FFL to variation of coefficients in each of the models considered.
In a wet lab it seems to be easier to identify a type of the FFL rather than to reveal the regulation details, which will be helpful in selection of the miRNA action model. The analysis of the temporal behaviour of the FFL together with the ability to find out a unique solution for every set of parameters may provide the way to select the most probable mechanism of miRNA action if the type of FFL is known. However, a noise in data can corrupt an ideal behaviour of FFL in absence of any perturbation.
It is widely believed that miRNA can buffer the consequences of noise in a cell. A simple mathematical model based on assumption that miRNA represses translation of its target was introduced to explore the ability of the 1In FFL to buffer fluctuations in upstream TF at steady state [34]. We demonstrated already that the behaviour of a FFL is modeldependent, and it is reasonable to study the ability of FFLs to buffer fluctuations in upstream regulator quantity on the assumption of the miRNA action, different from the translational repression. Besides, the results of our analysis also demonstrate an imperfection of approaches based on the analysis of model behaviour at steady state only: as it is shown in Figures 5, 6 and Figures 2 and 3 of Additional file 2 the dynamical behaviour of FFLs can be multivariant at early times even when the quantity of target protein at steady state is the same.
Consequently, we decide to investigate the ability of all FFLs with miRNA to buffer a noise caused by TF at all time moments and in all models. At each time moment we introduce the random fluctuations in TF quantity and measure how these fluctuations affect the target protein molecule number. We also consider how the variation of model coefficients influences the ability of the loops in noise damping. The efficiency of the FFLs in controlling the fluctuations of the target protein in response to noise introduced by TF depends on the number of TF molecules (which is a function of both k_{ w } and k_{ q }), the number of miRNA copies (depending on k_{ s } and h_{ s }) and the strength of miRNA action on target mRNA (defined by h_{ p }, h_{ g } and k_{ rs } coefficients in the models). We studied the ability of FFLs to buffer noise as a function of each of these three quantities, changing a corresponding coefficient and keeping fixed all others. For each combination of a model, loop and coefficient value we performed 100 numerical simulation runs to estimate the average value of the parameter e and the number of experiments with positive value of this parameter as described in section. Positive values of ε mean that a loop cannot dampen noise introduced by TF, and vice versa negative values of this coefficient testify the ability of the loop to reduce TF noise.
The 1In FFL shows the best ability to buffer noise introduced by TF: the noise is strongly decreased at the level of target protein production in all models and within a wide range of parameter variation Figure 7, Table 1. In all models the noise buffering increases with k_{ q } parameter for TF translation rate. At small values of k_{ q } the Stop model is more effective in noise reduction than two other models, while at large values of this parameter all models show similar ability to reduce noise. The variation of parameters that define the miRNA level in a loop, namely h_{ s } and k_{ s }, has small effect on ability of this loop to buffer noise in both the Target degradation and Dual degradation models. In the Stop model the maximal effect is achieved at intermediate values of these coefficients, the rate of noise reduction being the highest among all the models. Both in Target and Dual degradation models the ability of 1In FFL to buffer noise does not significantly change when a parameter, which defines the strength of miRNA action (h_{ g } or k_{ rs }), has large variation. In the Stop model the maximal effect is achieved at intermediate values of the h_{ p } coefficient.
The ability of the 2In loop to decrease noise introduced by TF essentially depends on the k_{ q } coefficient values: all the models are not able to efficiently buffer noise for large values of this coefficient (at k_{ q } = 0.16 the ε coefficient was positive in 9  14% of experiments), however, in the Stop and Target degradation models this effect becomes evident at larger k_{ q } values that in the other model (Table 2). In the Target degradation and Dual degradation models the variation of other coefficients has small effect on the ability of the loop to reduce noise. It is worth to note, that for larger values of both k_{ s } and h_{ s } the Stop model is more effective in noise buffering, than two other models. In frame of the Stop model the higher are h_{ s } values, the higher is the ability of 2In loop to reduce TF noise. On the contrary, the strongest noise reduction in this loop is achieved at intermediate values of the k_{ s } and h_{ p } coefficients.
Both 1C and 2C loops are bad buffers in frame of the Stop model, but not in both the Target and Dual degradation models (Figure 7), for which the variation of all coefficients (except of k_{ q } in both loops) exposes to the ability of these loops to buffer noise weaker than in the Stop model (Tables 3 and 4). In Stop model the behaviour of the 2C loop and, to smaller extent, of the 1C loop shows a strong dependence on parameters. Below we address the ability of coherent loops to buffer noise in detail.
In the Stop model the 1C loop is inefficient in reduction of TF noise at any value of the k_{ q } coefficient, showing the worst reduction at intermediate values (k_{ q } = 0.04) (Table 3). For large k_{ q } values (k_{ q } = 0.08 and higher) and in both the Target and Dual degradation models this loop loses the ability to efficiently reduce noise and starts to deal with the noise as in the Stop model.
The 1C loop is able to effectively buffer noise only for very small values of h_{ s } and k_{ s } coefficients in the Stop model: the larger are these coefficients, the less efficient is the noise reduction (Table 3). The loop shows nonability to reduce noise in small number of experiments at h_{ s } = 400 in both Target and Dual degradation models, as well as at intermediate k_{ s } = 0.25 value in the Target degradation model and at k_{ s } = 0.75 in the Dual degradation model.
In the Stop model the h_{ p } increase improves the ability of the loop to buffer noise as both the average ε value and the number of experiments, in which TF noise was not dampen (positive ε value) decrease (Table 3). In the Dual degradation model the 1C loop shows nonability to buffer TF noise in 25% of experiments.
The ability of 2C loop to reduce TF noise increases with the value of k_{ q } in the Stop model (Table 4). For large values of it (k_{ q } = 0.08 and higher) the TF noise was reduced in all numerical simulations. In two other models the TF noise is always dampen (with one exception of k_{ q } = 0.02 in the Dual degradation model, see (Table 4), however the efficiency of noise buffering increases as the value of the k_{ q } rises.
The dependence of the efficiency of noise damping on h_{ s } variation in the 2C loop and in frame of Stop model is complex: noise is efficiently reduced for very small h_{ s } values, while the efficiency of noise reduction decreases as the h_{ s } value grows, however the worst noise reduction happens at intermediate h^{s} values (h_{ s } = 200). In the Stop model the ability of 2C loop to buffer noise decreases with the growth of k_{ s }: at small k_{ s } ≤ 0.25 values the noise damping was observed in all simulations, while for higher k_{ s } values the fluctuations of the target protein molecules number are not reduced in many simulation runs (Table 4).
The 2C loop is unable to effectively buffer TF noise in all experiments when the values of the h_{ p } coefficient are small (below or equal to 60), at larger coefficient values this loop starts to efficiently buffer noise (Table 4). The noise is effectively reduced in the Dual degradation model of 2C loop for all values of the k_{ rs } coefficient, however the ability of dampening decreases as the value of this coefficient rises.
Conclusions
We performed here the theoretical analysis of a gene network subsystem, containing a FFL mediated by TF and miRNA. We have shown that different mechanisms of miRNA action lead to a variety of types of dynamical behaviour of FFLs during cell cycle and govern their ability to dampen noise caused by TF fluctuations.
The molecular mechanisms of miRNA action are not clear so far, and we elaborate three mathematical models introduced in [34], that describe the gene expression in miRNA mediated FFL under the assumption of different mechanisms of miRNA action. In the Stop model miRNA represses translation of its target mRNA, in the Target degradation model miRNA promotes the target mRNA degradation, and in the Dual degradation model miRNA is not reused, but degrades along with target mRNA.
Due to an intrinsic complexity and nonlinearity of biological systems it is a hard task to obtain any analytic solution to differential equations, which describe the regulation in FFL with miRNA under the models considered. These equations are nonlinear, fortunately, we were able to obtain the exact solutions to some of them, namely, to those, describing target mRNA and miRNA production in the Stop model and for several biologically relevant values, i.e.,

for slow degradation of miRNA (when miRNA is degraded two times slower than TF or its mRNA),

fast degradation of miRNA (when miRNA is degraded two times faster than TF or its mRNA) and

very fast degradation of miRNA (when miRNA is degraded three times faster than TF) (see Additional file 1 for details).
Despite of the fact that miRNAs are generally stable molecules, it was shown recently that individual miRNAs may be exposed to an accelerated decay [18].
We used the exact solutions to check the results of numerical simulations obtained for all other coefficient sets and initial conditions, and proven the validity of these results. In general, the exact solutions obtained can be used as a genuine check point in numerical simulations of larger networks with miRNA embedded into a FFL motif.
We have rigorously proven the uniqueness of solutions to all the models under consideration, i.e., in the models considered there is the onetoone correspondence between the given parameter set and the solution, describing the dynamics of target protein production in FFL.
Study of cell components behaviour at steady state is conventional, however it is worth to note that the steady states are not completely informative even in nonbiotic systems, consisting of almost identical clusters of several atoms/molecules. In biology the steady states are not unique: very different pathways may lead to the same stationary position. From the general viewpoint of dynamic control theory an early time seems to be the most promising one for tentative control/influence onto the loop dynamics, either by noise or by any external factor. That is why we first examined the FFL dynamics quantitatively, and at the whole time interval of cell cycle, alternatively to recent qualitative consideration [34].
We have shown that FFLs mediated by miRNA and TF may have many possible outcomes, depending on interaction between the loop elements. The target protein profiles can take different forms, which are unambiguously defined by initial conditions and model coefficients. This can be illustrated by analysis of the behaviour of one and the same FFL under different models, and when both initial conditions and the level of a target protein at steady state are the same. Due to the difference in mechanisms of miRNA action the behaviour of a FFL in time will be quite different in frame of different models. This situation is reconstructed for the 1C loop in Figure 8. It is evident that the behaviour of the models at early times is different: the maximum of target protein production is the smallest in the Dual degradation model, while the time, at which this protein reaches the steady state is the largest one for the Stop model. It is noteworthy that these graphs correspond to models with different degradation coefficients of the target mRNA and protein, i.e. to different biological situations.
Contemporary experimental set up seems hardly be able to capture many facets of miRNA function. Indeed, in spite of a bunch of publications, describing the crucial role of miRNA in control of many biological processes and in progression of various diseases, the molecular mechanisms of miRNA action are still not evident [14, 15]. In a wet lab it is easier to identify a type of FFL, rather than to reveal the details of regulation, which will be helpful in selection of the miRNA action model. Then the analysis of temporal behaviour of the FFL together with the ability to find out a unique solution for every set of parameters may help to select the most feasible mechanism of miRNA action for the type of FFL given.
The results obtained allow us to propose the following strategy for an ideal FFL: let us consider how a miRNA acts, having an information about the FFL topology and a given set of quantitative measurements of each player in the loop. If the initial quantities of molecules for each player in FFL are not known, we should firstly determine them in experiments. Using our analytical results, we shall be able to calculate in detail the temporal dependence for each player in the FFL and in each possible model for various values of coefficients. Next, we identify that graph, which passes through the given set of the (molecule quantities) points among all others. By virtue of the uniqueness theorem this graph will correspond to most feasible regulation type in the FFL among all models considered.
However FFLs are subjected to noise influence, and for this reason we studied the noise buffering in both coherent and incoherent FFLs, that led to conclusion that incoherent FFLs are better noise buffers than coherent ones (Figure 7). Moreover, even the parameter variation does not seriously affect the noise buffering ability of incoherent FFLs. Therefore the selection strategy proposed seems to be suitable to predict a model for incoherent FFLs.
The FFL dynamic behaviour analysis performed for different models and parameter sets shows that an extremal value of target protein quantity (max or min, depending on a loop) may be accomplished for intermediate values of coefficients. This is valid also for the noise buffering problem, that required a simple analysis of the Hill sigmoid behaviour, see Additional file 1.
In general, an action of any disturbance, e.g., a noise, having sufficiently small amplitude, (ca. 10  15%; otherwise it would prevail over a signal in FFL) depends on the Hill function type and can be described as follows. For the Hill sigmoid functions, governing either an activation or a repression, the noise will be recognizable for intermediate molecule quantity. By virtue of analysis in the whole time interval we found that when the noise amplitude is high, but the total amount (the regular one plus the noise associated fluctuations) of miRNA is intermediate (the Hill function is far from any of almost constant limits), then the noise is translated directly into the protein production process and will not be suppressed by a loop.
When the noise amplitude is still high, however, the total amount of miRNA is sufficiently high, too (the Hill function is close to any of limits), then the noise level is relatively small to disturb the protein production and will be suppressed by a loop.
Methods
Many details of mathematical analysis are given in Section and in the Supporting Information files; here we briefly describe the numerical simulations. The coefficients used in mathematical modeling were mostly taken from [34], and for the Stop model they are, as follows:
The other coefficients used in both the Degradation and the Dual degradation models are:
Numerical simulation of dynamics in each loop was made for the time period of t = 5000 seconds. The initial numbers of components in each loop were equal to one half of their steady state values. The dynamics of solutions at early stages depends not only on coefficients, but on initial values, too. To analyze how the quantity of target protein molecules depends on the numbers of TF and miRNA molecules at initial time moment these numbers were changed from 0 to one quarter and one half, as well as to the values of two and four times higher, than at steady state (in normalized quantities).
For the noise simulations we introduced a vector z with the correlation function $K\left(\tau \right)=\mathsf{\text{exp}}\left(\left\tau \right\right)$, that was obtained in a form:
where ξ is another vector, obtained by the MATLAB procedure, which contains the Gaussian white noise. The vector z represents the model of a stationary Markov process, and we may add its values at the moments t_{ i } to the number of molecules of TF, which is used in equations for the quantities of miRNA and target mRNA, and solve these equations numerically. In fact, the procedure mentioned is based on the UhlenbeckOrnstein process (1930), the only one, which is stationary random, Gaussian and Marcovian simultaneously. The process is widely used now in mathematical modelling of noise instead of any numerical version of the so called "white noise". The particle velocity in this process is finite, and taking into account the difference between a random force and the idealised white noise, one may provide a finite acceleration of the particle, too.
We calculated the numbers q of TF molecules, of miRNA (s), of the target mRNA (r) and of target protein (p) with the noise component and without it, and afterwards we found the values ${N}_{A}=ma{x}_{t}\mathsf{\text{}}q\left(t\right){q}_{n}\left(t\right)\mathsf{\text{}}{N}_{B}=ma{x}_{t}\mathsf{\text{}}p\left(t\right){p}_{n}\left(t\right)$ to estimate the relative power of noise. Here q(t),p(t) are the numbers of molecules of transcription factor and target protein without the noise component, respectively, while q_{ n }(t), p_{ n }(t) are corresponding quantities with noise. Therefore we may find the value of the estimation parameter ε:
to conclude whether the relative noise level in target protein quantity is higher or lower, than in TF. Consequently, ε < 0 means that the noise level in target protein is lower than in TF, and the noise is buffered in the loop.
We performed calculations with 100 different noise vectors z for each loop and each model, estimated ε in every calculation, studied and calculated an average value of ε, as well as the number of calculations with positive ε value. We used the Wilcoxon test for these averaged values of ε to compare results for the same loops and same models under variation of the model coefficients.
References
 1.
He L, Hannon GJ: MicroRNAs: small RNAs with a big role in gene regulation. Nat Rev Genet. 2004, 5: 522531. 10.1038/nrg1379.
 2.
Harfe BD: MicroRNAs in vertebrate development. Curr Opin Genet Dev. 2005, 15 (4): 410415. 10.1016/j.gde.2005.06.012.
 3.
Bushati N, Cohen SM: microRNA functions. Annu Rev Cell Dev Biol. 2007, 23: 175205. 10.1146/annurev.cellbio.23.090506.123406.
 4.
Avraham R, Yarden Y: Regulation of signalling by microRNAs. Biochem Soc Trans. 2012, 40 (1): 2630. 10.1042/BST20110623.
 5.
Calin GA, Ferracin M, Cimmino A, Di Leva G, Shimizu M, Wojcik SE, lorio MV, Visone R, Sever NI, Fabbri M, Luliano R, Palumbo T, Pichiorri F, Roldo C, Garzon R, Sevignani C, Rassenti L, Alder H, Volinia S, Liu CG, Kipps TJ, Negrini M, Croce CM: A microRNA signature associated with prognosis and progression in chronic lymphocytic leukemia. N Engl J Med. 2005, 353: 17931801. 10.1056/NEJMoa050995.
 6.
AlvarezGarcia I, Miska EA: MicroRNA functions in animal development and human disease. Development. 2005, 132: 46534662. 10.1242/dev.02073.
 7.
Beezhold KJ, Castranova V, Chen F: Microprocessor of microRNAs: regulation and potential for therapeutic intervention. Mol Cancer. 2010, 9: 134
 8.
Flynt AS, Lai EC: Biological principles of microRNAmediated regulation: shared themes amid diversity. Nature Reviews Genetics. 2008, 9 (11): 83142. 10.1038/nrg2455.
 9.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS: MicroRNA targets in Drosophila. Genome Biol. 2003, 5 (1): 1
 10.
Brennecke J, Stark A, Russell RB, Cohen SM: Principles of microRNAtarget recognition. PLoS Biol. 2005, 1 (1): 1310.1371/journal.pcbi.0010013.
 11.
Grun D, Wang YL, Langenberger D, Gunsalu KC, Rajewsky N: microRNA target predictions across seven drosophila species and comparison to mammalian targets. PLoS Comput Biol. 2005, 1 (1): 1310.1371/journal.pcbi.0010013.
 12.
Hornstein E, Shomron N: Canalization of development by microRNAs. Nat Genet. 2006, 38: 2024. 10.1038/ng1803.
 13.
Friedman RC, Farh KKH, Burge CB, Bartel DP: Most mammalian mRNAs are conserved targets of microRNAs. Genome Research. 2008, 19 (1): 92105. 10.1101/gr.082701.108.
 14.
Baek D, Villen J, Shin C, Camargo FD, Gygi SP, Bartel DP: The impact of microRNAs on protein output. Nature. 2008, 455 (7209): 6471. 10.1038/nature07242.
 15.
Selbach M, Schwanhausser B, Thierfelder N, Fang Z, Khanin R, Rajewsky N: Widespread changes in protein synthesis induced by microRNAs. Nature. 2008, 455 (7209): 5863. 10.1038/nature07228.
 16.
Bazzini AA, Lee MT, Giraldez AJ: Ribosome profiling shows that miR430 reduces translation before causing mRNA decay in zebrafish. Science. 2012, 336 (6078): 2337. 10.1126/science.1215704.
 17.
Guo H, Ingolia NT, Weissman JS, Bartel DP: Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature. 2010, 466: 835841. 10.1038/nature09267.
 18.
Ruegger S, Grosshans H: Micro rna turnover: when, how, and why. Trends Biochem Sci. 2012, 37 (10): 436446. 10.1016/j.tibs.2012.07.002.
 19.
Baccarini A, Chauhan H, Gardner TJ, Jayaprakash AD, Sachidanandam R, Brown BD: Kinetic analysis reveals the fate of a microRNA following target regulation in mammalian cells. Current Biology. 2011, 21 (5): 369376. 10.1016/j.cub.2011.01.067.
 20.
Nissan T, Parker R: Computational analysis of mirnamediated repression of translation: implications for models of translation initiation inhibition. RNA. 2008, 14 (8): 148091. 10.1261/rna.1072808.
 21.
Morozova N, Zinovyev A, Nonne N, Pritchard LL, Gorban AN, HarelBellan A: Kinetic signatures of microrna modes of action. RNA. 2012, 18 (9): 163555. 10.1261/rna.032284.112.
 22.
Zinovyev A, Morozova N, Gorban AN, HarelBelan A: Mathematical modeling of micrornamediated mechanisms of translation repression. Adv Exp Med Biol. 2013, 774: 189224. 10.1007/9789400755901_11.
 23.
Zhang Z, Qin YW, Brewer G, Jing Q: Microrna degradation and turnover: regulating the regulators wires. RNA. 2012, 3: 593600.
 24.
Herranz H, Cohen SM: MicroRNAs and gene regulatory networks: managing the impact of noise in biological systems. Genes & Development. 2010, 24 (13): 133944. 10.1101/gad.1937010.
 25.
Stark A, Brennecke J, Bushati N, Russell RB, Cohen SM: Animal MicroRNAs confer robustness to gene expression and have a significant impact on 3'UTR evolution. Cell. 2005, 123 (6): 11331146. 10.1016/j.cell.2005.11.023.
 26.
Farh KK, Grimson A, Jan C, Lewis BP, Johnston WK, Lim LP, Burge CB, Bartel DP: The widespread impact of mammalian microRNAs on mRNA repression and evolution. Science. 2005, 310: 18171821. 10.1126/science.1121158.
 27.
Sood P, Krek A, Zavolan M, Macino G, Rajewsky N: Celltypespecific signatures of microRNAs on target mRNA expression. Proc Natl Acad Sci USA. 2006, 103: 27462751. 10.1073/pnas.0511045103.
 28.
Giraldez AJ, Mishima Y, Rihel J, Grocock RJ, Van Dongen S, Inoue K, Enright AJ, Schier AF: Zebrafish mir430 promotes deadenylation and clearance of maternal mRNAs. Science. 2006, 312: 7579. 10.1126/science.1122689.
 29.
Li X, Cassidy JJ, Reinke CA, Fischboeck S, Carthew RW: A microRNA imparts robustness against environmental fluctuation during development. Cell. 2010, 137 (2): 273282.
 30.
ShenOrr SS, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of escherichia coli. Sciencet. 2002, 31 (1): 6468.
 31.
Milo R, ShenOrr SS, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: simple building blocks of complex networks. Nat Genet. 2002, 298: 824827.
 32.
Alon U: An introduction to systems biology. Design principles of biological circuits. 2006, Chapman & Hall/CRC, 301
 33.
Yu X, Lin J, Zack DJ, Mendell JT, Qian J: Analysis of regulatory network topology reveals functionally distinct classes of microRNAs. Nucleic Acids Research. 2008, 36: 64946503. 10.1093/nar/gkn712.
 34.
Osella M, Bosia C, Cora D, Caselle M: The role of incoherent microRNAmediated feedforward loops in noise buffering. PLoS Computational Biology. 2011, 7 (3): 100110110.1371/journal.pcbi.1001101.
 35.
Schwanhausser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, Chen W, Selbach M: Global quantification of mammalian gene expression control. Nature. 2011, 473 (7347): 337342. 10.1038/nature10098.
 36.
Li JJ, Bickel PJ, Biggin MD: System wide analyses have underestimated protein abundances and the importance of transcription in mammals. PeerJ. 2014, 2: 270
 37.
von Dassow G, Meir E, Munro EM, Odell GM: The segment polarity network is a robust developmental module. Nature. 2000, 406: 188192. 10.1038/35018085.
 38.
Kitano H: Biological robustness. Nat Rev Genet. 2004, 5 (11): 826837. 10.1038/nrg1471.
Acknowledgements
We are thankful to Sergey Berezin for valuable advises. This work is partly supported by RFBR grants 140100334, 140401522, EC collaborative project HealthF52010260429 SysPatho and the Russian Federation Ministry of Education and Science program 51002020.
Declarations
Publication of this article has been funded by the BGRS\SB2014 Organizing Committee and the Russian Science Foundation (RSF).
This article has been published as part of BMC Genomics Volume 15 Supplement 12, 2014: Selected articles from the IX International Conference on the Bioinformatics of Genome Regulation and Structure\Systems Biology (BGRS\SB2014): Genomics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/15/S12.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
AS and MS proposed and designed the study, MD and AS obtained all the mathematical solutions to the problem, MD and MS analyzed the modeling results. AS and MS wrote and typeset the manuscript with considerable input from MD.
Rights and permissions
About this article
Published
DOI
Keywords
 dynamic behaviour
 feedforward loop
 miRNA
 noise buffering
 nonlinear equations