Dynamics of miRNA driven feed-forward loop depends upon miRNA action mechanisms

Abstract

Background

We perform the theoretical analysis of a gene network sub-system, composed of a feed-forward 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 re-used, but degrades along with target mRNA. We examine the feed-forward 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 feed-forward loops. In particular, we found that the ability of feed-forward 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 non-coding 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 miRNA-binding 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 desta-bilization 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 rate-limiting 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 anti-correlative 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 feed-forward 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.

Figure 1

The incoherent and coherent feed-forward loops. Arrows mean activation, the turned over T-bars indicate repression. TF-transcription factor, miR-miRNA, Target - target protein. A: 1In - type 1 incoherent FFL, TF activates both target mRNA and miRNA synthesis. B: 1C - type 1 coherent FFL, TF represses traget mRNA and activates miRNA synthesis. C: 2In - type 2 incoherent FFL, TF represses both target mRNA and miRNA synthesis. D: 2C - type 2 coherent FFL, TF activates target mRNA and represses miRNA synthesis.

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 miRNA-mediated 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 non-linear 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 non-uniqueness of solutions. This non-uniqueness 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. Non-uniqueness 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 non-uniqueness 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 sub-system 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 co-existence of distinct miRNA-mediated 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 re-used 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 feed-forward loop (FFL)

$$\left\{\begin{array}{cc}\frac{dw}{dt}=k_w-g_{w}w,\hfill & \left(mRNAofTF\right)\hfill \\ \frac{dq}{dt}=k_{q}w-g_{q}q,\hfill & \left(TF\right)\hfill \\ \frac{ds}{dt}=k_s\left(q\right)-g_{s}s,\hfill & \left(miRNA\right)\hfill & \cdot \hfill \\ \frac{dr}{dt}=k_r\left(q\right)-g_{r}r,\hfill & \left(mRNA\right)\hfill \\ \frac{dp}{d}=k_p\left(s\right)r-g_{p}p.\hfill & \left(targetprotein\right)\hfill \end{array}\right.$$

(1)

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 (1-4) 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:

$$y\left(x\right)=\frac{\alpha x^c}{{\beta }^c+x^c},$$

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 non-linear 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:

$$k_s\left(q\right)=\frac{k_s{q}^2}{h_s^2+q^2};k_r\left(q\right)=\frac{k_r{h}_r^2}{q^2+h_r^2};$$

(2)

while for type 1 incoherent FFL the production function for target mRNA will be different:

$$k_r\left(q\right)=\frac{k_r{q}^2}{h_r^2+q^2}.$$

(3)

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 miRNA-mRNA 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 re-used. This complex degradation constant will be denoted as k _rs , and the resulting non-linear (due to the multiplicative term rs) equations will have a form:

$$\begin{array}{c}\frac{dw}{dt}=k_w-g_{w}w;\frac{dq}{dt}=k_{q}w-g_{q}q;\frac{ds}{dt}=k_s\left(q\right)-g_{s}s-k_{rs}rs;\hfill \\ \frac{dr}{dt}=k_r\left(q\right)-g_{r}r-k_{rs}rs;\frac{dp}{dt}=k_{p}r-g_{s}p.\hfill \end{array}$$

(4)

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 non-linear ODE, which coefficients explicitly depend on initial values. Main advantage of exact solutions to corresponding non-linear 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 one-to-one correspondence between given parameters and solutions.

Comparative analysis of FFL temporal behaviour under different models

The mechanisms underlying miRNA-mediated 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 non-linear 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 bell-shaped 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.

Figure 2

The solutions to the Stop model in various FFLs. 1In - type 1 incoherent FFL, 1C-type 1 coherent FFL, 2In - type 2 incoherent FFL, 2C - type 2 coherent FFL; w, q, s, r and p denote graphs of solutions for TF mRNA, TF, miRNA, target mRNA and target protein, resp. A - D: temporal dynamics of absolute number of each molecule species is presented. E - H: the molecules numbers for each species are normalized on steady state values to better visualize the behaviour of RNA species.

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 wave-like 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 pulse-like 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.

Figure 3

The solutions to the Target degradation model in various FFLs. 1In - type 1 incoherent FFL, 1C-type 1 coherent FFL, 2In - type 2 incoherent FFL, 2C - type 2 coherent FFL; w, q, s, r and p denote graphs of solutions for TF mRNA, TF, miRNA, target mRNA and target protein correspondingly. A - D: the temporal dynamics of absolute number of each molecule species is presented, E - H: molecules numbers for each species are normalized on steady state values to better visualize the behaviour of RNA species.

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 re-used. 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 pulse-like 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 mRNA-miRNA 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).

Figure 4

The variation of target protein quantities in different models of 1 In loop in response to variation of h _s , k _q and initial miRNA molecules numbers. The values of h _s coefficient defining the amount of TFs, at which the transcription rate of miRNA gene is half of its maximum value, were within 0 - 400 mol. interval, the translation rates k _q for TF were taken from 0 - 0.16sec⁻¹ interval, initial quantities of miRNA were changed as described in section. Left column - Stop model, central column - Target degradation model, right column - Dual degradation model. A - C: In all models the quantity of target protein increases as h _s rise. D: In the Stop model the largest number of target protein molecules is observed at intermediate values of the kq coefficient. E-F: In the Target and Dual degradation models kq increase leads to increase of target protein quantity. G and I: In both Stop and Dual degradation models the form of target protein profile changes from the bell-shaped to the U-shaped one as the initial number of miRNA molecules rises. H: In the Target degradation model all profiles show moderate dependence on the change of the initial number of miRNA molecules.

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 bell-shaped to the U-shaped 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 U-shaped form to the bell-shaped 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 bell-shaped 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.

Figure 5

The variation of target protein quantities in different models of 1 C loop in response to variation of h _s , k _q and initial miRNA quantities. The parameter values are the same as of Fig. Left column - Stop model, central column - Target degradation model, right column - Dual degradation model. A-C: The dependence of target protein profiles on h _s variation. D - F: The dependence of target protein profiles on k _q variation. G - I: The dependence of target protein profiles on initial miRNA quantities. Initial quantities of miRNA were changed as described in section.

Figure 6

The variation of target protein quantities in different models of 2 C loop in response to variation of hp, h _g , k _rs and h _s coefficients, and initial miRNA and TF quantities. Parameters h _p , h _g , k _rs define the action of miRNA on target mRNA, h _s is the dissociation coefficient. The miRNA and TF initial quantities were changed as described in section. Left column - Stop model, central column - Target degradation model, right column - Dual degradation model. A: In the Stop model the quantity of target protein increases as the value of h _p increases from 0 to 240 molecules. B: In the Target degradation model the quantity of target protein increases as the value of h _g increases from 0 to 240 molecules. C: The increase of the k _rs coefficient value from 0 to 8 × 10⁻⁵ mol⁻¹sec⁻¹ results in the fall of the target protein quantity in the Dual degradation model. D - F: In all models the increase of h _s from 0 to 400 leads to decrease in target protein quantities. textbfG - I: The dependence of target protein quantities on initial miRNA molecules number in the models. J - L: The patterns of dependence of target protein quantities on initial TF quantities in the models.

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 model-dependent, 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.

Figure 7

The ability of FFLs to buffer noise introduced by TF in frame of different models. In each panel the values of parameter ε calculated in 100 experiments are shown as narrow vertical lines. ε < 0 means that the noise is buffered in a loop, ε > 0 means non-ability of the loop to dampen noise. The coefficients and initial conditions for each experiment are given in section. A-C: Type 1 incoherent loop is able to buffer noise introduced by TF in all models. D - F: Type 2 incoherent loop buffers TF noise in all models. G: Type 1 coherent loop is a bad buffer in frame of the Stop model. I, H: Type 1 coherent loop is able to dampen TF noise in frame of the Target and Dual degradation models. J: Type 2 coherent loop is a bad buffer in frame of the Stop model. K, L: Type 2 coherent loop is able to buffer TF noise under the Target and Dual degradation models.

Table 1 The ability of 1In loop to buffer TF noise under variation of parameters in all three models.

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.

Table 2 The ability of 2In loop to buffer TF noise under of parameters in all three models.

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.

Table 3 The ability of 1C loop to buffer TF noise under of parameters in all three models.

Table 4 The ability of 2C loop to buffer TF noise under of parameters in all three models.

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 non-ability 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 non-ability to buffer TF noise in 2-5% 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 sub-system, 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 re-used, but degrades along with target mRNA.

Due to an intrinsic complexity and non-linearity 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 non-linear, 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 one-to-one 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 non-biotic 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.

Figure 8

The target protein profiles in 1 C FFL are shown in different models for identical both initial conditions and the steady state level of a target protein. Solid line - Stop model, dashed line-Target degradation model, dotted line - Dual degradation model.

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:

$$\begin{array}{c}{k}_w=0.06,g_w=0.004,kq=0.04,g_q=0.002,k_s=0.5,h_s=h_r=200,\hfill \\ g_s=0.002,k_r=0.08,g_r=0.004,k_p=0.8,g_p=0.002,h_p=h_p=60,\hfill \\ g_{max}=0.004,k_{rs}gr\mathsf{\text{/}}200.\hfill \end{array}$$

(5)

The other coefficients used in both the Degradation and the Dual degradation models are:

$$h_g=60,g_{max}=0.004,k_{rs}=0.00002.$$

(6)

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:

$$z\left(i\right)=\mathsf{\text{exp}}\left(-0.1\right)z\left(i-1\right)+\sqrt{1-\mathsf{\text{exp}}\left(-0.2\right)\xi \left(i\right)},$$

(7)

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 Uhlenbeck-Ornstein 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 ε:

$$\epsilon =\left(N_B-N_A\right)\mathsf{\text{/}}{N}_A,$$

(8)

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.