Abstract
Background
Protein translation is a multistep process which can be represented as a cascade of biochemical reactions (initiation, ribosome assembly, elongation, etc.), the rate of which can be regulated by small noncoding microRNAs through multiple mechanisms. It remains unclear what mechanisms of microRNA action are the most dominant: moreover, many experimental reports deliver controversial messages on what is the concrete mechanism actually observed in the experiment. Nissan and Parker have recently demonstrated that it might be impossible to distinguish alternative biological hypotheses using the steady state data on the rate of protein synthesis. For their analysis they used two simple kinetic models of protein translation.
Results
In contrary to the study by Nissan and Parker, we show that dynamical data allow discriminating some of the mechanisms of microRNA action. We demonstrate this using the same models as developed by Nissan and Parker for the sake of comparison but the methods developed (asymptotology of biochemical networks) can be used for other models. We formulate a hypothesis that the effect of microRNA action is measurable and observable only if it affects the dominant system (generalization of the limiting step notion for complex networks) of the protein translation machinery. The dominant system can vary in different experimental conditions that can partially explain the existing controversy of some of the experimental data.
Conclusions
Our analysis of the transient protein translation dynamics shows that it gives enough information to verify or reject a hypothesis about a particular molecular mechanism of microRNA action on protein translation. For multiscale systems only that action of microRNA is distinguishable which affects the parameters of dominant system (critical parameters), or changes the dominant system itself. Dominant systems generalize and further develop the old and very popular idea of limiting step. Algorithms for identifying dominant systems in multiscale kinetic models are straightforward but not trivial and depend only on the ordering of the model parameters but not on their concrete values. Asymptotic approach to kinetic models allows putting in order diverse experimental observations in complex situations when many alternative hypotheses coexist.
Background
MicroRNAs (miRNAs) are currently considered as key regulators of a wide variety of biological pathways, including development, differentiation and oncogenesis. Recently, remarkable progress was made in understanding of microRNA biogenesis, functions and mechanisms of action. Mature microRNAs are incorporated into the RISC effector complex, which includes as a key component an Argonaute protein. MicroRNAs affect gene expression by guiding the RISC complex toward specific target mRNAs. The exact mechanism of this inhibition is still a matter of debate. In the past few years, several mechanisms have been reported, and some of the reports contradict to each other (for review, see [13]). The inhibition mechanisms include, in particular, the inhibition of translation initiation (acting at the level of cap40S or 40SAUG60S association steps), the inhibition of translation elongation and the premature translation termination. MicroRNAmediated mRNA decay and sequestration of target mRNAs in Pbodies have been also proposed. Moreover, some microRNAs mediate target mRNA cleavage [4], chromatin reorganization followed by transcriptional repression or activation [5,6], or translational activation [7,8].
The most frequently reported, but also much debated, is the mechanism of gene repression by microRNAs which occurs at the level of mRNA translation. At this level, several mode of actions have been suggested (see Fig. 1). Historically, the first proposed mechanism was the inhibition of translation elongation. The major argument supporting this hypothesis was the observation that the inhibited mRNA remained associated with the polysomal fraction (in which mRNAs are associated with polysomes). This observation was reproduced in different systems [913]. The idea of a postinitiation mechanism was further supported by the observation that some mRNAs can be repressed by a microRNA even when their translation is capindependent (mRNAs with an IRES or Acapped) [11,1416]. Although it was initially proposed that the ribosomes were somehow "frozen" on the mRNA, it is important to note that it is difficult to discriminate experimentally between different potential postinitiation mechanisms, such as elongation inhibition, premature ribosome dissociation ("ribosome dropoff") or normal elongation but with nascent polypeptide degradation. The last proposition (this mechanism can occur in conjunction with the two others) is supported by the fact that the mRNApolysomal association is puromycinsensitive, indicating that it depends on a peptidyltransferase activity [13,17]. However, no nascent peptide has ever been experimentally demonstrated; thus its degradation would occur extremely rapidly after the synthesis [10,11,18]. The premature ribosome dissociation is supported by the decreased readthrough of inhibited mRNA [11]. The ribosomal dropoff and/or ribosomal "slowing" are supported by the slight decrease in the number of associated ribosomes observed in some studies [10,13].
Figure 1. Interaction of microRNA with protein translation process. Four mechanisms of translation repression which are considered in the mathematical modeling are indicated: 1) on the initiation process, preventing assembling of the initiation complex; 2) on a late initiation step, such as searching for the start codon; 3) on the ribosome assembly; 4) on the translation process. There exist other mechanisms of microRNA action on protein translation (transcriptional, transport to Pbodies, ribosome dropoff, cotranslational protein degradation and others) that are not considered in this paper. Here 40S and 60S are light and heavy components of the ribosome, 80S is the assembled ribosome bound to mRNA, eIF4F is an translation initiation factor, PABC1 is the PolyA binding protein, "cap" is the mRNA cap structure needed for mRNA circularization, RISC is the RNAinduced silencing complex.
Concurrently, several reports have been published indicating an action of microRNAs at the level of initiation. An increasing number of papers reports that microRNAtargeted mRNAs shift towards the light fractions in polysomal profiles [1820]. This shows a decrease of mature translating ribosomes, suggesting that microRNAs act at the initiation step. Moreover, several reports show that microRNAmediated inhibition is relieved when translation is driven by a capindependent mechanism such as IRESmRNA or AcappedmRNA [18,20,21]. This observation was confirmed in several invitro studies [2225]. In particular, in one of those, an excess of eIF4F could relieve the inhibition, and the inhibition led to the decreased 80S in the polysomal gradient [22].
Most of the data indicating a shift towards the light polysomal fraction or the requirement for a capdependent translation are often interpreted in favour of involvement of microRNAs at early steps of translation, i.e., cap binding and 40S recruitment. However, some of them are also compatible with a block at the level of 60S subunit joining. This hypothesis is also supported by invitro experiments showing a lower amount of 60S relative to 40S on inhibited mRNAs. Moreover, toeprinting experiments show that 40S is positioned on the AUG [26]. Independently, it was shown that eIF6, an inhibitor of 60S joining, is required for microRNA action [27], but this was in contradiction with other studies [2].
Thus, the data on the exact step of translational inhibition are clearly contradictory. Taking also into account the data about mRNA degradation and Pbodies localization, it is difficult to draw a clear picture of the situation, and the exact mechanism by which microRNA represses mRNA expression is highly controversial, not mentioning the interrelations between the different mechanisms and their possible concomitant action. Several attempts to integrate the different hypotheses have been made [13,2830]. For example, it was proposed that one mechanism could act as a "primary" effect, and the other as a "secondary" mechanism, either used to reinforce the inhibition or as a backup mechanism. In others, the different mechanisms could all coexist, but occur differentially depending on some yet unidentified characteristics. For example, it has been observed than the same mRNA targeted by the same microRNA can be regulated either at the initiation or the elongation step depending on the mRNA promoter and thus on the mRNA nuclear history [31]. It was also proposed that technical (experimental) problems, including the variety of experimental systems used, may also account for these discrepancies [13]. However, this possibility does not seem to be sufficient to provide a simple and convincing explanation to the reported discrepancies.
A possible solution to exploit the experimental observations and to provide a rational and straightforward data interpretation is the use of mathematical models for microRNA action on protein translation. For many years, the process of protein synthesis is a subject of mathematical modeling with use of various approaches from chemical kinetics and theoretical physics. Many of the models created consider several stages of translation, however, most of them concentrate on the elongation and termination processes. In [3234], a nonequilibrium statistical physics description of protein synthesis was proposed. Models taking into account gene sequence were developed in [33,3537]. These models can predict the probability of that a ribosome will completely terminate a transcript, spatiotemporal organization of ribosomes in polysome, dependence of the protein synthesis rate on various factors, such as presence of slow synonymous codons in the gene sequence [33,37] and the frequency of nonsense errors [35]. Several models of the effect of microRNA on protein translation were proposed. Thus, in [38] the authors tried to determine which inhibition mechanism (via translation repression or transcript degradation) is the most abundant in mammalian cells using Bayesian modeling and microarray data. Quantitative features of sRNAmediated gene regulation were considered in [39]. A simple kinetic model of microRNAmediated mRNA degradation was proposed in [40] and compared to a temporal microarray dataset.
In this paper we will analyze two simple models of microRNA action on protein translation developed recently by Nissan and Parker [41]. They studied the microRNAdependent steady states rates of protein synthesis [41] and provided a critical analysis for the experiments with alternative mRNA cap structures and IRES elements [22,23,25]. This analysis led to a possible explanation of the conflicting results. The authors suggested that the relief of translational repression upon replacement of the cap structure can be explained if microRNA is acting on a step which is not ratelimiting in the modified system, in which case, the effect of microRNA can simply not be observed. It was claimed that it is impossible to discriminate between two alternative interpretations of the biological experiments with cap structure replacement, using sole monitoring of the steady state level of protein synthesis [41].
Two remarks can be made in this regard. Firstly, in practice not only the steady state level of protein can be observed but also other dynamical characteristics, such as the relaxation time, i.e. the time needed to achieve the steady state rate after a perturbation (such as restarting the translation process). We argue that having these measurements in hands, one can distinguish between two alternative interpretations. In this paper we provide such a method from the same models as constructed by Nissan and Parker, for comparison purposes. However, the method applied can be easily generalized for other models.
Secondly, even in the simple nonlinear model of protein translation, taking into account the recycling of ribosomal components, it is difficult to define what is the rate limiting step. It is known from the theory of asymptotology of biochemical networks [42] that even in complex linear systems the "rate limiting place" notion is not trivial and cannot be reduced to a single reaction step. Moreover, in nonlinear systems the "rate limiting place" can change with time and depend on the initial conditions. Hence, conclusions of [41] should be reconsidered for the nonlinear model, made more precise and general. The notion of rate limiting step should be replaced by the notion of dominant system.
In this paper we perform careful analysis of the Nissan and Parker's models and provide their approximate analytical solutions, which allows us to generalize the conclusions of [41] and make new checkable predictions on the identifiability of active mechanism of microRNAdependent protein translation inhibition.
The paper is organized in the following way. The Methods contain introduction, all necessary definitions and basic results of the asymptotology of biochemical reaction networks (quasiequilibrium, quasi steadystate, limiting step and dominant system asymptotics), used further in the Results. The Methods section is deliberately made rather detailed to make the reading selfsufficient. These details are necessary for reproducing the analytical calculations but not for understanding the interpretation of the modeling results. When the most important notions are introduced (such as dominant system, critical parameters), they are emphasized in bold. The Results section starts with listing model assumptions, followed by deriving semianalytical solutions of Nissan and Parker's model and interpretation of the analysis results and prediction formulations. For those readers who are interested only in the applied aspect of this work, it is possible to skip the details of deriving the analytical solutions and start reading from the "Model assumptions" section, look at the definition of the model parameters and variables and continue with 'the 'Effect of microRNA on the translation dynamics" section.
Results
Model assumptions
We consider two models of action of microRNA on protein translation process proposed in [41]: the simplest linear model, and the nonlinear model which explicitly takes into account recycling of ribosomal subunits and initiation factors.
Both models, of course, are significant simplifications of biological reality. Firstly, only a limited subset of all possible mechanisms of microRNA action on the translation process is considered (see Fig. 1). Secondly, all processes of synthesis and degradation of mRNA and microRNA are deliberately neglected. Thirdly, interaction of microRNA and mRNA is simplified: it is supposed that when microRNA is added to the experimental system then only mRNA with bound microRNAs are present (this also assumes that the concentration of microRNA is abundant with respect to mRNA). Concentrations of microRNA and mRNA are supposed to be constant. Interaction of only one type of microRNA and one type of mRNA is considered (not a mix of several microRNAs). The process of initiation is greatly simplified: all initiation factors are represented by only one molecule which is marked as eIF4F.
Finally, the classical chemical kinetics approach is applied, based on solutions of ordinary differential equations, which supposes sufficient and wellstirred amount of both microRNAs and mRNAs. Another assumption in the modeling is the mass action law assumed for the reaction kinetic rates.
It is important to underline the interpretation of certain chemical species considered in the system. The ribosomal subunits and the initiation factors in the model exist in free and bound forms, moreover, the ribosomal subunits can be bound to several regions of mRNA (the initiation site, the start codon, the coding part). Importantly, several copies of fully assembled ribosome can be bound to one mRNA. To model this situation, we have to introduce the following quantification rule for chemical species: amount of "ribosome bound to mRNA" means the total number of ribosomes translating proteins, which is not equal to the number of mRNAs with ribosome sitting on them, since one mRNA can hold several translating ribosomes (polyribosome). In this view, mRNAs act as places or catalyzers, where translation takes place, whereas mRNA itself formally is not consumed in the process of translation, but, of course, can be degraded or synthesized (which is, however, not considered in the models described further).
The simplest linear protein translation model
The simplest representation of the translation process has the form of a circular cascade of reactions [41] (see Fig. 2).
Figure 2. The simplest model of microRNA action on the protein translation. The simplest model of microRNA action on the protein translation, represented with use of Systems Biology Graphical Notation (a) and schematically with the condition on the constants (b). The two mechanisms of microRNA action (capdependent and capindependent) are depicted.
The list of chemical species in the model is the following:
1. 40S, free small ribosomal subunit.
2. mRNA:40S, small ribosomal subunit bound to the initiation site.
3. AUG, small ribosomal subunit bound to the start codon.
The catalytic cycle is formed by the following reactions:
1. 40S → mRNA:40S, Initiation complex assembly (rate k_{1}).
2. mRNA:40S → AUG, Some late and capindependent initiation steps, such as scanning the 5'UTR for the start AUG codon recognition (rate k_{2}) and 60S ribosomal unit joining.
3. AUG → 40S, combined processes of protein elongation and termination, which leads to production of the protein (rate k_{3}), and fall off of the ribosome from mRNA.
The model is described by the following system of equations [41]:
where Prsynth(t) is the rate of protein synthesis.
Following [41], let us assume that k_{3 }>> k_{1}, k_{2}. This choice was justified by the following statement: "...The subunit joining and protein production rate (k_{3}) is faster than k_{1 }and k_{2 }since mRNA:40S complexes bound to the AUG without the 60S subunit are generally not observed in translation initiation unless this step is stalled by experimental methods, and elongation is generally thought to not be rate limiting in protein synthesis..." [41].
Under this condition, the equations (1) have the following approximate solution (which becomes the more exact the smaller the ratio):
for the initial condition
From the solution (2) it follows that the dynamics of the system evolves on two time scales: 1) fast elongation dynamics on the time scale ≈ 1/k_{3}; and 2) relatively slow translation initiation dynamics with the relaxation time . The protein synthesis rate formula (3) does not include the k_{3 }rate, since it is neglected with respect to k_{1}, k_{2 }values. From (3) we can extract the formula for the protein synthesis steadystate rate Prsynth (multiplier before the parentheses) and the relaxation time t_{rel }for it (inverse of the exponent power):
Now let us consider two experimental situations: 1) the rates of the two translation initiation steps are comparable k_{1 }≈ k_{2}; 2) the capdependent rate k_{1 }is limiting: k_{1 }<< k_{2}. Accordingly to [41], the second situation can correspond to modified mRNA with an alternative capstructure, which is much less efficient for the assembly of the initiation factors, 40S ribosomal subunit and polyA binding proteins.
For these two experimental systems (let us call them "wildtype" and "modified" correspondingly), let us study the effect of microRNA action. We will model the microRNA action by diminishing the value of a kinetic rate coefficient for the reaction representing the step on which the microRNA is acting. Let us assume that there are two alternative mechanisms: 1) microRNA acts in a capdependent manner (thus, reducing the k_{1 }constant) and 2) microRNA acts in a capindependent manner, for example, through interfering with 60S subunit joining (thus, reducing the k_{2 }constant). The dependence of the steady rate of protein synthesis and the relaxation time on the efficiency of the microRNA action (i.e., how much it is capable to diminish a rate coefficient) is shown on Fig. 3.
Figure 3. Predicted change in the steadystate rate of protein synthesis and its relaxation time. Graphs illustrating the predicted change in the steadystate rate of protein synthesis (left), and its relaxation time, i.e., the time needed to recover from a perturbation to the steady state value (right). Four curves are presented. The black ones are for the wildtype cap structure, which is modeled by k_{1 }= k_{2}. The red ones are for the modified structure, when k_{1 }<<k_{2}. The main conclusion from the left graph is that if microRNA acts on a late initiation step, diminishing k_{2}, then its effect is not measurable unless k_{2 }is very strongly suppressed (as reported in [41]). The main conclusion from the right graph is that the effect of microRNA can be measurable in this case if one looks at dynamical features such the relaxation time.
Interestingly, experiments with cap structure replacement were made and the effect of microRNA action on the translation was measured [22,23]. No change in the protein rate synthesis after applying microRNA was observed. From this it was concluded that microRNA in this system should act through a capdependent mechanism (i.e., the normal "wildtype" cap is required for microRNA recruitment). It was argued that this could be a misinterpretation [41] since in the "modified" system, capdependent translation initiation is a rate limiting process (k_{1 }<< k_{2}), hence, even if microRNA acts in the capindependent manner (inhibiting k_{2}), it will have no effect on the final steady state protein synthesis rate. This was confirmed this by the graph similar to the Fig. 3a.
From the analytical solution (2) we can further develop this idea and claim that it is possible to detect the action of microRNA in the "modified" system if one measures the protein synthesis relaxation time: if it significantly increases then microRNA probably acts in the capindependent manner despite the fact that the steady state rate of the protein synthesis does not change (see the Fig. 3b). This is a simple consequence of the fact that the relaxation time in a cycle of biochemical reactions is limited by the second slowest reaction (see [42] or the "Dominant system for a simple irreversible catalytic cycle with limiting step" section in Methods). If the relaxation time is not changed in the presence of microRNA then we can conclude that none of the two alternative mechanisms of microRNAbased translation repression is activated in the system, hence, microRNA action is dependent on the structure of the "wildtype" transcript cap.
The observations from the Fig. 3 are recapitulated in the Table 1. This analysis (of course, oversimplified in many aspects) provides us with an important lesson: observed dynamical features of the translation process with and without presence of microRNA can give clues on the mechanisms of microRNA action and help to distinguish them in a particular experimental situation. Theoretical analysis of the translation dynamics highlights what are the important characteristics of the dynamics which should be measured in order to infer the possible microRNA mechanism.
Table 1. Modeling two mechanisms of microRNA action in the simplest linear model
This conclusion suggests the notion of a kinetic signature of microRNA action mechanism which we define as the set of measurable characteristics of the translational machinery dynamics (features of time series for protein, mRNA, ribosomal subunits concentrations) and the predicted tendencies of their changes as a response to microRNA action through a particular biochemical mechanism.
The nonlinear protein translation model
To explain the effect of microRNA interference with translation initiation factors, a nonlinear version of the translation model was proposed [41] which explicitly takes into account recycling of initiation factors (eIF4F) and ribosomal subunits (40S and 60S).
The model contains the following list of chemical species (see also Fig. 4):
Figure 4. Nonlinear model of microRNA action on the protein translation. Nonlinear model of microRNA action on the protein translation, represented with use of Systems Biology Graphical Notation (a) and schematically with the condition on the constants (b). The difference from the simplest model (Fig. 2) is in the explicit description of initiation factors eIF4F, and ribosomal subunits 40S and 60S recycling.
1. 40S, free 40S ribosomal subunit.
2. 60S, free 60S ribosomal subunit.
3. eIF4F, free initiation factor.
4. mRNA:40S, formed initiation complex (containing 40S and the initiation factors), bound to the initiation site of mRNA.
5. AUG, initiation complex bound to the start codon of mRNA.
6. 80S, fully assembled ribosome translating protein.
There are four reactions in the model, all considered to be irreversible:
1. 40S + eIF4F → mRNA:40S, assembly of the initiation complex (rate k_{1}).
2. mRNA:40S → AUG, some late and capindependent initiation steps, such as scanning the 5'UTR by for the start codon AUG recognition (rate k_{2}).
3. AUG → 80S, assembly of ribosomes and protein translation (rate k_{3}).
4. 80S → 60S+40S, recycling of ribosomal subunits (rate k_{4}).
The model is described by the following system of equations [41]:
where [40S] and [60S] are the concentrations of free 40S and 60S ribosomal subunits, [eIF4F] is a concentration of free translation initiation factors, [mRNA : 40S] is the concentration of 40S subunit bound to the initiation site of mRNA, [AUG] is the concentration of the initiation complex bound to the start codon, [80S] is the concentration of ribosomes translating protein, and Prsynth is the rate of protein synthesis.
The model (6) contains three independent conservations laws:
The following assumptions on the model parameters were suggested in [41]:
with the following justification: "...The amount 40S ribosomal subunit was set arbitrarily high ... as it is thought to generally not be a limiting factor for translation initiation. In contrast, the level of eIF4F, as the canonical limiting factor, was set significantly lower so translation would be dependent on its concentration as observed experimentally... Finally, the amount of subunit joining factors for the 60S large ribosomal subunit were estimated to be more abundant than eIF4F but still substoichiometric when compared to 40S levels, consistent with in vivo levels... The k_{4 }rate is relatively slower than the other rates in the model; nevertheless, the simulation's overall protein production was not altered by changes of several orders of magnitude around its value..." [41].
Notice that further in our paper we show that the last statement about the value of k_{4 }is needed to be made more precise: in the model by Nissan and Parker, k_{4 }is a critical parameter (see "Asymptotology and dynamical limitation theory for biochemical reaction networks" section in Methods). It does not affect the steady state protein synthesis rate only in one of the possible scenarios (inefficient initiation, deficit of the initiation factors).
Steady state solution
The final steady state of the system can be calculated from the conservation laws and the balance equations among all the reaction fluxes:
where "s" index stands for the steady state value. Let us designate a fraction of the free [60S] ribosomal subunit in the steady state as . Then we have
and the equation to determine x, in which we have neglected the terms of smaller order of magnitude, based on conditions (10):
From the inequalities on the parameters of the model, we have δ > 1, γ << 1 and, if k_{1 }>> k_{4}/[eIF4F]_{0 }then α << β. From these remarks it follows that the constant term γ(1 β) of the equation (13) should be much smaller than the other polynomial coefficients, and the equation (13) should have one solution close to zero and two others:
provided that α << 1  δ or α << 1  β. In the expression for x_{1 }we cannot neglect the term proportional to a, to avoid zero values in (13).
The solution x_{2 }is always negative, which means that one can have one positive solution x_{0 }<< 1 if and two positive solutions x_{0 }and x_{1 }if . However, from (12), (14) and (10) it is easy to check that if x_{1 }> 0 then x_{0 }does not correspond to a positive value of [eIF4F]_{s}. This means that for a given combination of parameters satisfying (10) we can have only one steady state (either x_{0 }or x_{1}).
The two values x = x_{0 }and x = x_{1 }correspond to two different modes of translation. When, for example, the amount of the initiation factors [eIF4F]_{0 }is not enough to provide efficient initiation (, x = x_{1}) then most of the 40S and 60S subunits remain in the free form, the initiation factor [eIF4F] being always the limiting factor. If the initiation is efficient enough , then we have x = x_{0 }<< 1 when almost all 60S ribosomal subunits are engaged in the protein elongation, and [eIF4F] being a limiting factor at the early stage, however, is liberated after and ribosomal subunits recycling becomes limiting in the initiation (see the next section for the analysis of the dynamics).
Let us notice that the steady state protein synthesis rate under these assumptions is
This explains the numerical results obtained in [41]: with low concentrations of [eIF4F]_{0 }microRNA action would be efficient only if it affects k_{2 }or if it competes with eIF4F for binding to the mRNA cap structure (thus, effectively further reducing the level [eIF4F]_{0}) With higher concentrations of [eIF4F]_{0}, other limiting factors become dominant: [60S]_{0 }(availability of the heavy ribosomal subunit) and k_{4 }(speed of ribosomal subunits recycling which is the slowest reaction rate in the system). Interestingly, in any situation the protein translation rate does not depend on the value of k_{1 }directly (of course, unless it does not become "globally" rate limiting), but only through competing with eIF4F (which makes the difference with the simplest linear protein translation model).
Equation (15) explains also some experimental results reported in [22]: increasing the concentration of [eIF4F] translation initiation factor enhances protein synthesis but its effect is abruptly saturated above a certain level.
It would be interesting to make some conclusions on the shift of the polysomal profile from the steady state solutions (14). In this model, the number of ribosomes sitting on mRNA N_{polysome }is defined by , where [mRNA] is the concentration of mRNA. However, [mRNA] is not an explicit dynamical variable in the model, it is implicitly included in other model constants, such as k_{1}, together with the effective volume of cytoplasmic space considered in the model. Nevertheless, the model can predict the relative shift of the polysome profile. In the steady state
and N_{polysome }changes in the same way as the protein synthesis steady state value.
Analysis of the dynamics
It was proposed to use the following model parameters in [41]: k_{1 }= k_{2 }= 2, k_{3 }= 5, k_{4 }= 1, [40S]_{0 }= 100, [60S]_{0 }= 25, [eIF4F]_{0 }= 6. As we have shown in the previous section, there are two scenarios of translation possible in the Nissan and Parker's model which we called "efficient" and "inefficient" initiation. The choice between these two scenarios is determined by the critical combination of parameters . For the original parameters from [41], β = 0.48 < 1 and this corresponds to the simpler onestage "inefficient" initiation scenario. To illustrate the alternative situation, we changed the value of k_{4 }parameter, putting it to 0.1, which makes β = 4.8 > 1. The latter case corresponds to the "efficient" initiation scenario, the dynamics is more complex and goes in three stages (see below).
Simulations of the protein translation model with these parameters and the initial conditions
are shown on the Fig. 5. The system shows nontrivial relaxation process which takes place in several epochs. Qualitatively we can distinguish the following stages:
Figure 5. Simulation of the nonlinear protein translation model. Simulation of the nonlinear protein translation model with parameters k_{1 }= 2, k_{2 }= 2, k_{3 }= 5, k_{4 }= 0.1, [40S]_{0 }= 100, [60S]_{0 }= 25, [eIF4F]_{0 }= 6. a) and b) chemical species concentrations at logarithmic and linear scales; c) and d) reaction fluxes at logarithmic and linear scales. By the dashed line several stages are delimited during which the dynamics can be considered as (pseudo)linear. To determine where ">>" and "<<" conditions are violated, we arbitrarily consider "much bigger" or "much smaller" as difference in one order of magnitude (by factor 10).
1) Stage 1: Relatively fast relaxation with conditions [40S] >> [eIF4F], [60S] >> [AUG]. During this stage, the two nonlinear reactions 40S + eIF4F → mRNA : 40S and AUG + 60S → 80S can be considered as pseudomonomolecular ones: eIF4F → mRNA : 40S and AUG → 80S with rate constants dependent on [40S] and [60S] respectively. This stage is characterized by rapidly establishing quasiequilibrium of three first reactions (R1, R2 and R3 with k_{1}, k_{2 }and k_{3 }constants). Biologically, this stage corresponds to assembling of the translation initiation machinery, scanning for the start codon and assembly of the first full ribosome at the start codon position.
2) Transition between Stage 1 and Stage 2.
3) Stage 2: Relaxation with the conditions [40S] >> [eIF4F], [60S] << [AUG]. During this stage, the reactions 40S + eIF4F → mRNA : 40S and AUG + 60S → 80S can be considered as pseudomonomolecular eIF4F → mRNA : 40S and 60S → 80S. This stage is characterized by two local quasisteady states established in the two network reaction cycles (formed from R1R2 and R3R4 reactions). Biologically, this stage corresponds to the first round of elongation, when first ribosomes moves along the coding region of mRNA. The small ribosomal subunit 40S is still in excess which keeps the initiation stage (reaction R1R2 fluxes) relatively fast.
4) Transition between Stage 2 and Stage 3.
5) Stage 3: Relaxation with the conditions [40S] << [eIF4F], [60S] << [AUG]. During this stage, the reactions 40S + eIF4F → mRNA : 40S and AUG + 60S → 80S can be considered as pseudomonomolecular 40S → mRNA : 40S and 60S → 80S. During this stage all reaction fluxes are balanced. Biologically, this stage corresponds to the stable production of the protein with constant recycling of the ribosomal subunits. Most of ribosomal subunits 40S are involved in protein elongation, so the initiation process should wait the end of elongation for that they would be recycled.
Stages 13 can be associated with the corresponding dominant systems [42] which are shown on Fig. 6. Below we give a more detailed analysis of stages 13 and transitions between them.
Figure 6. Dominant systems for three stages of relaxation. Dominant systems for three stages of relaxation of the model (6). Stage 1) The dominant system is a pseudolinear network of reactions. Stage 2) The dominant system is a quasisteady state approximation, where one supposes that the fluxes in two network cycles are balanced. Stage 3) The dominant system is a pseudolinear network of reactions.
Stage 1: translation initiation and assembly of the first ribosome at the start codon
The dominant system of the Stage 1 (Fig. 6a) can be modeled as a linear system of equations (notice that it is not equivalent to the system of equations that would correspond to fully monomolecular reaction network because the reaction R2 is still bimolecular despite the fact that the products of this reaction do not interact, which leads to the linear description):
where = k_{1 }· [40S], = k_{3 }· [60S] and we consider that at this stage the changes of 40S and 60S are relatively slow. This system has simple approximate solution, taking into account constraints on the parameters k_{2 }<< , ; k_{4 }<< , , k_{2}, also assuming k_{2 }<<   , and for the initial condition
we have
From this solution, one can conclude that the relaxation of this model goes at several time scales (very rapid ~ min(1/, 1/) and slow ~1/k_{4}) and that when eIF4F, mRNA:40S and AUG already reached their quasiequilibrium values, [80S] continues to grow. This corresponds to the quasiequilibrium approximation asymptotics (see the "Quasi steadystate and quasiequilibrium asymptotics" section of the Methods). At some point 80S will reach such a value that it would be not possible to consider 60S constant: otherwise the conservation law (9) will be violated. This will happen when [80S] << [60S] condition is not satisfied anymore, i.e., following our convention to consider "much smaller" as difference in one order of magnitude, at . The same consideration is applicable for another conservation law (7) in which [80S] is included, but from the time point . From [40S]_{0 }>> [eIF4F]_{0 }and [40S]_{0 }> [60S]_{0 }we have min(t', t") = t'. This means that the parameters [40S], [60S] of the (local) steady states for [eIF4F] and [AUG] should slowly (at the same rate as [80S]) change from the time point t' (variable [mRNA: 40S] does not change because its local steady state does not depend on [40S], [60S]). In other words, after t = t' the Stage 1 solution (20) should be prolonged as
From these equations, one can determine the effective duration of the Stage 1: by definition, it will be finished when one of the two conditions ([40S] >> [eIF4F], [60S] >> [AUG]) will be violated, which happens at times ~ and ~ correspondingly, hence, the second condition will be violated first (from [60S]_{0 }< [40S]_{0}).
Stage 2: first stage of protein elongation, initiation is still rapid
The Stage 2 is characterized by conditions [eIF4F] << [40S], [60S] << [AUG]. This fact can be used for deriving the quasisteady state approximation: we assume that the reaction fluxes in two network cycles (R1R2 and R3R4) are independently balanced:
Then (6) is simplified and, using the conservation laws, we have a single equation on [40S]:
where A = [40S](t) + [AUG](t) is a constant quantity conserved accordingly to the quasisteady state approximation (see "Quasi steadystate and quasiequilibrium asymptotics" section of Methods). Equation (23) can be already integrated but let us further simplify it for our analysis. Having in mind k_{4 }<< k_{3 }and assuming that at the beginning of the Stage 2 [AUG] >> , we can simplify (23) to
and, further, assuming that at the beginning of the Stage 2 we have [40S] >> [AUG] let us approximate the righthand side of the equation by a piecewiselinear function
where [40S]_{t=t" }is the amount of 40S at the beginning of the Stage 2. Then the descent of [40S](t) can be separated into linear and exponential phases:
where K_{1}, K_{2 }are linear and exponential slopes and [40S]_{s2 }is the quasisteady state value of [40S] at the end of the Stage 2:
Other dynamic variables are expressed through [40S](t) as
At some point, the amount of free small ribosomal subunit 40S, which is abundant at the beginning of the Stage 2, will not be sufficient to support rapid translation initiation. Then the initiation factor eIF4F will not be the limiting factor in the initiation and the condition [40S] >> [eIF4F] will be violated. We can estimate this time as .
Stage 3: steady protein elongation, speed of initiation equals to speed of elongation
During the Stage 3 all fluxes in the network become balanced and the translation arrives at the steady state. From Fig. 6 it is clear that the relaxation goes independently in the cycle R3  R4, where the relaxation equations are simply
where t"' is the time when the Stage 3 of the relaxation starts. This relaxation goes relatively fast, since k_{3}[AUG]_{t=t''' }is relatively big. So, during the Stage 3, one can consider the cycle R1  R4 equilibrated, with [80S] = [80]_{s}, [60S] = [60]_{s }values.
Hence, the relaxation during the Stage 3 consists in redistributing concentrations of 40S and mRNA:40S to their steady states in a linear chain of reactions R1  R2 (the value of [AUG] is relatively big and can be adjusted from the conservation law (7)). Using the pseudolinear approximation of this stage (see Fig. 6), we can easily write down the corresponding approximate relaxation equations:
where B = (1  ) ([40]_{s} [40S]_{t = t'''}). [40]_{s }and [mRNA : 40S]_{s }are the steadystate values of the corresponding variables, see (12). The values [60S]_{t = t''', }[eIF4F]_{t = t''', }[AUG]_{t = t''' }and [mRNA : 40S]_{t = t″' }can be estimated from (28), using the [40S]_{t = t''' }value. The relaxation time at this stage equals
The solution for the Stage 3 can be further simplified if k_{2 }<< k_{1}[eIF4F]_{t = t''' }or k_{2 }>> k_{1 }[eIF4F]_{t = t'''}.
Transitions between stages
Along the trajectory of the dynamical system (6) there are three dominant system each one transforming into another. At the transition between stages, two neighbor dominant systems are united and then split. Theoretically, there might be situations when the system can stay in these transition zones for long periods of time, even infinitely. However, in the model (6) this is not the case: the trajectory rapidly passes through the transition stages and jumps into the next dominant system approximation.
Three dominant approximations can be glued, using the concentration values at the times of the switching of dominant approximation as initial values for the next stage. Note that the Stage 2 has essentially one degree of freedom since it can be approximated by a single equation (23). Hence, one should only know one initial value [40S]_{t = t'' }to glue the Stages 1 and 2. The same is applied to the gluing of Stages 2 and 3, since in the end of Stage 2 all variable values are determined by the value of [40S]_{t = t'''}.
Case of always limiting initiation
As it follows from our analysis, the most critical parameter of the nonlinear protein translation model is the ratio . Above we have considered the case β > 1 which is characterized by a switch of the limiting factor in the initiation (from eIF4F at the Stages 1 and 2 to 40S at the Stage 3).
In the case β < 1 the dynamics becomes simpler and consists of one single stage: relaxation accordingly to (20) and further with correction (21) with the relaxation time ~ (the quasiequilibrium approximation corresponding to the Stage 1 works well for the whole translation process). The reason for this is that if the initiation is not efficient then the system is never in the situation of the Stage 2 conditions when the cycle R1R2 is balanced with much bigger flux than the cycle R3R4. This approximation is the more exact the smaller β value, however, the value of β should not be necessary very small. For example, for the default parameter values of the model β = 0.48, and it well reproduces the dynamics (see Fig. 7cd). From numerical experiments one can see that even for β = 0.95 the dynamics is qualitatively well reproduced. To model the Acap structure effect with very weak capacity for initiation (assembly of the initiation factors and 40S subunit), we should also consider the case
Figure 7. Comparison of the numerical and approximate analytical solutions of the nonlinear protein translation model. Examples of the exact numerical (circles) and approximate analytical (solid lines) solutions of the nonlinear protein translation model. a) For the set of parameters k_{1 }= 2, k_{2 }= 2, k_{3 }= 5, k_{4 }= 0.1; b) For parameters k_{1 }= 1, k_{2 }= 5, k_{3 }= 50, k_{4 }= 0.01; c) For the set of parameters from [41], k_{1 }= 2, k_{2 }= 2, k_{3 }= 5, k_{4 }= 1; d) Reaction fluxes for the set of parameters c). Dashed black vertical lines denote evaluated transition points between the dynamics stages. Dashed red vertical points denote the time points where [40S](t) = 10·[eIF4F](t) and [40S](t) = [eIF4F](t)/10 respectively.
for which the solution derived above is not directly applicable. However, the analytical calculations in this case can be performed in the same fashion as above. The detailed derivation of the solution is given in Additional file 1. The effect of putting k_{1 }very small on the steady state protein synthesis and the relaxation time is shown on Fig. 7.
Additional file 1. Analytical analysis of the case of very inefficient cap structure. In this text we derive an asymptotical solution for the case when k_{1 }is very small corresponding to the case of very inefficient translation initiation (for example, in the case of Acap structure replacement experiment)
Format: PDF Size: 155KB Download file
This file can be viewed with: Adobe Acrobat Reader
In a similar way all possible solutions of the equations (6) with very strong inhibitory effect of microRNA on a particular translation step can be derived. These solutions will describe the situation when the effect of microRNA is so strong that it changes the dominant system (limiting place of the network) by violating the initial constraints (10) on the parameters (for example, by making k_{3 }smaller than other k_{i}s). Such possibility exists, however, it can require too strong (nonphysiological) effect of microRNAdependent translation inhibition.
Effect of microRNA on the translation dynamics
Our analysis of the nonlinear Nissan and Parker's model showed that the protein translation machinery can function in two qualitatively different modes, determined by the ratio . We call these two modes "efficient initiation" (β > 1) and "inefficient initiation" (β < 1) scenarios. Very roughly, this ratio determines the balance between the overall speeds of initiation and elongation processes. In the case of "efficient initiation" the rate of protein synthesis is limited by the speed of recycling of the ribosomal components (60S). In the case of "inefficient initiation" the rate of protein synthesis is limited by the speed of recycling of the initiation factors (eIF4F). Switching between two modes of translation can be achieved by changing the availability of the corresponding molecules ([60S]_{0 }or [eIF4F]_{0}) or by changing the critical kinetic parameters (k_{2 }or k_{4}). For example, changing k_{4 }from 1 (Fig. 7c) to 0.1 (Fig. 7a), performs such a switch for the original parameter values from [41].
As a result of the dynamical analysis, we can assemble an approximate solution of the nonlinear system under assumptions (10) about the parameters. An example of the approximate solution is given on Fig. 7. The advantage of such a semianalytical solution is that one can predict the effect of changing the system parameters. For example, on Fig. 7b the solution is compared to an exact numerical one, where the parameters have been changed but still obey the initial constraints (10).
One of the obvious predictions is that the dynamics of the system is not sensitive to variations of k_{3}, so if microRNA acts on the translation stage controlled by k_{3 }then no microRNA effect could be observed looking at the system dynamics (being the fastest one, k_{3 }is not a critical parameter in any scenario).
If microRNA acts on the translation stage controlled by k_{4 }(for example, by ribosome stalling mechanism) then we should consider two cases of efficient (β > 1) and inefficient (β < 1) initiation. In the first case the steady state protein synthesis rate is controlled by k_{4 }(as the slowest, limiting step) and any effect on k_{4 }would lead to the proportional change in the steady state of protein production. By contrast, in the case of inefficient initiation, the steady state protein synthesis is not affected by k_{4}. Instead, the relaxation time is affected, being ~ . However, diminishing k_{4 }increases the β parameter, hence, this changes "inefficient initiation" scenario for the opposite, hence, making k_{4 }critical for the steady state protein synthesis anyway when k_{4 }becomes smaller than . For example, for the default parameters of the model, decreasing k_{4 }value firstly leads to no change in the steady state rate of protein synthesis, whereas the relaxation time increases and, secondly, after the threshold value starts to affect the steady state protein synthesis rate directly (see Fig. 8A, B). This is in contradiction to the message from [41] that the change in k_{4 }by several orders of magnitude does not change the steady state rate of protein synthesis.
Figure 8. Effect of mimicking different mechanisms of miRNA action on translation. Effect of decreasing some model parameters mimicking different mechanisms of miRNA action on translation. Relaxation time here is defined as the latest time at which any chemical species in the model differs from its final steady state by 10% A) and B) correspond to the scenario with "inefficient" initiation, with use of the model parameters proposed in [41] (k_{1 }= k_{2 }= 2, k_{3 }= 5, k_{4 }= 1, [eIF4F]_{0 }= 6, [60S]_{0 }= 25, [40S]_{0 }= 100), which gives β = 0.48 < 1. C) and D) correspond to the scenario with "efficient" initiation, with parameters (k_{1 }= 2, k_{2 }= 3, k_{3 }= 50, k_{4 }= 0.1, [eIF4F]_{0 }= 6, [60S]_{0 }= 25, [40S]_{0 }= 100), which gives β = 7.2 > 1. The absciss value indicates the degree of inhibition (decreasing) of a parameter. EH) same as AD) but for a modified cap structure, modeled by reduced k_{1 }parameter: k_{1 }= 0.01 for these curves, the other parameters are the same as on AD) correspondingly.
Analogously, decreasing the value of k_{2 }can convert the "efficient" initiation scenario into the opposite after the threshold value . We can recapitulate the effect of decreasing k_{2 }in the following way. 1) in the case of "efficient" initiation k_{2 }does not affect the steady state protein synthesis rate up to the threshold value after which it affects it in a proportional way. The relaxation time drastically increases, because decreasing k_{2 }leads to elongation of all dynamical stages duration (for example, we have estimated the time of the end of the dynamical Stage 2 as ). However, after the threshold value the relaxation time decreases together with k_{2}, quickly dropping to its unperturbed value (see Fig. 8CD). 2) in the case of "inefficient" initiation the steady state protein synthesis rate depends proportionally on the value of k_{2 }(15), while the relaxation time is not affected (see Fig. 8AB).
MicroRNA action on k_{1 }directly does not produce any strong effect neither on the relaxation time nor on the steady state protein synthesis rate. This is why in the original work [41] capdependent mechanism of microRNA action was taken into account through effective change of the [eIF4F]_{0 }value (total concentration of the translation initiation factors), which is a critical parameter of the model (see 15).
The effect of microRNA on various mechanism and in various experimental settings (excess or deficit of eIF4F, normal cap or Acap) is recapitulated in Table 2. The conclusion that can be made from this table is that all four mechanisms show clearly different patterns of behavior in various experimental settings. From the simulations one can make a conclusion that it is still not possible to distinguish between the situation when microRNA does not have any effect on protein translation and the situation when it acts on the step which is neither rate limiting nor "second rate limiting" in any experimental setting (k_{3 }in our case). Nevertheless, if any change in the steadystate protein synthesis or the relaxation time is observed, theoretically, it will be possible to specify the mechanism responsible for it.
Table 2. Modeling of four mechanisms of microRNA action in the nonlinear protein translation model
Available experimental data and possible experimental validation
It is important to underline that the Nissan and Parker's models analyzed in this paper are qualitative descriptions of the protein translation machinery. The parameter values used represent rough orderofmagnitude estimations or real kinetic rates. Moreover, these values should be considered as relative and unitless since they do not match any experimental time scale (see below). Nevertheless, such qualitative description already allows to make predictions on the relative changes of the steady states and relaxation times (see the Table 2), and in principle these predictions can be verified experimentally. Let us imagine an experiment in which it would be possible to verify such predictions. In this experiment, two time series should be compared: 1) one measured in a system in which microRNA acts on a normal "wildtype" protein translation machinery and 2) another system almost fully identical to the first one but in which one of the translation stages is modified (made slow and ratelimiting, or, opposite, very rapid). There are multiple possibilities to modify the rate of this or that translation stage. The initiation can be affected by changing the concentration of the initiation factors such as eIF4F as in [22]. The scanning stage can be affected by introducing various signals in the 5'UTR sequence of mRNA such as inframe AUG codons (see, for example, [43]). In principle, the elongation stage can be modified by introducing slow synonymous codons in the coding sequence (there even exist mathematical models of their effect [33,35,37] that can be used for the optimal experiment design). The stage of elongation termination can be influenced by varying the concentration of the corresponding release factors (ETF1 or ETF2), at least in vitro. The two time series measured after activation or introduction of microRNA should be characterized for the relative changes of steady state values and relaxation times of protein and mRNA concentrations, and, if possible, the number of ribosomes in the polysome. Also, ideally, it is desired to construct several experimental systems in which the amount of inhibition by microRNA can be gradually changed (for example, by changing the number of the corresponding seed sequences in the 3'UTR region).
To the best of our knowledge, there is no such a dataset published until so far, even partially. In several recent papers, one can find published time series of protein and mRNA concentrations or their relative changes measured after introducing microRNA. For example, the deadenylation time course is shown in [25]: translation decreases after 20 min and stops at 30 min, deadenylation begin at 30 min, goes around 1 h. In [44], the authors study the kinetics of degradation of mRNA. After adding microRNA to the system, the amount and the length of the targeted mRNA starts to decrease at around 35 hours, and decreases by 90% at 8 hours. In [45], the authors study the global change of protein after transfection of a microRNA. They described a small change at the mRNA level at 8 h after miRNA transfection, and the considerable decrease appeared only after 32 hours while the protein concentration change was apparent at the timecourse between 8 hours and 32 hours. In the in vitro system used in [22], at 15 min after incubation with microRNA there was already a 25% decrease of translation, indicating that the translational inhibition can be a relatively rapid mechanism.
These data on protein translation kinetics show that the relaxation time range could vary from several minutes to several hours and even tens of hours depending on the critical step affected, on various mRNA properties and on the whole biological system taken for the experiment (for example, the presence or absence of different effectors influencing different steps of the translation process). These data should be taken into account when constructing more realistic and quantitative models of microRNA action on protein translation.
Discussion
The role of microRNA in gene expression regulation is discovered and confirmed since ten years, however, there is still a lot of controversial results regarding the role of concrete mechanisms of microRNAmediated protein synthesis repression. Some authors argue that it is possible that the different modes of microRNA action reflect different interpretations and experimental approaches, but the possibility that microRNAs do indeed silence gene expression via multiple mechanisms also exists. Finally, microRNAs might silence gene expression by a common and unique mechanism; and the multiple modes of action represent secondary effects of this primary event [13].
The main reason for accepting a possible experimental bias could be the studies in vitro, where conditions are strongly different from situation in vivo. Indeed, inside the cell, mRNAs (microRNA targets) exist as ribonucleoprotein particles or mRNPs, and second, all proteins normally associated with mRNAs transcribed in vivo are absent or at least much different from that bound to the same mRNA in an in vitro system or following the microRNAs transfection into cultured cells. The fact that RNAbinding proteins strongly influence the final outcome of microRNA regulation is proved now by several studies [19,46,47]. The mathematical results provided in this paper suggests a complementary view on the coexistence of multiple microRNAmediated mechanisms of translation repression. Mathematical modeling suggests to us to ask a question: if multiple mechanisms act simultaneously, would all of them equally contribute to the final observable repression of protein synthesis or its dynamics? The dynamical limitation theory gives an answer: the effect of microRNA action will be observable and measurable in two cases: 1) if it affects the dominant system of the protein translationary machinery, or 2) if the effect of microRNA action is so strong that it changes the limiting place (the dominant system).
In a limited sense, this means, in particular, that the protein synthesis steady rate is determined by the limiting step in the translation process and any effect of microRNA will be measurable only if it affects the limiting step in translation, as it was demonstrated in [41]. Due to the variety of external conditions, cellular contexts and experimental systems the limiting step in principle can be any in the sequence of events in protein translation, hence, this or that microRNA mechanism can become dominant in a concrete environment. However, when put on the language of equations, the previous statement already becomes nontrivial in the case of nonlinear dynamical models of translation (and even linear reaction networks with nontrivial network structure). Our analysis demonstrates that the limiting step in translation can change with time, depends on the initial conditions and is not represented by a single reaction rate constant but rather by some combination of several model parameters. Methodology of dynamical limitation theory that we had developed [42,48], allows to deal with these situations on a solid theoretical ground.
Furthermore, in the dynamical limitation theory, we generalize the notion of the limiting step to the notion of dominant system, and this gives us a possibility to consider not only the steady state rate but also some dynamical features of the system under study. One of the simplest measurable dynamical feature is the protein synthesis relaxation time, i.e. the time needed for protein synthesis to achieve its steady state rate. The general idea of "relaxation spectrometry" goes back to the works of Manfred Eigen, a Nobel laureate [49] and is still underestimated in systems biology. Calculation of the relaxation time (or times) requires careful analysis of time scales in the dynamical system, which is greatly facilitated by the recipes proposed in [42,50]. As we have demonstrated in our semianalytical solutions, measuring the steady state rate and relaxation time at the same time allows to detect which step is possibly affected by the action of microRNA (resulting in effective slowing down of this step). To our knowledge, this idea was never considered before in the studies of microRNAdependent expression regulation. The Table 2 recapitulates predictions allowing to discriminate a particular mechanism of microRNA action.
Conclusions
The analysis of the transient dynamics gives enough information to verify or reject a hypothesis about a particular molecular mechanism of microRNA action on protein translation. For multiscale systems only that action of microRNA is distinguishable which affects the parameters of dominant system (critical parameters), or changes the dominant system itself. Dominant systems generalize and further develop the old and very popular idea of limiting step. Algorithms for identifying dominant systems in multiscale kinetic models are straightforward but not trivial and depend only on the ordering of the model parameters but not on their concrete values. Asymptotic approach to kinetic models of biological networks suggests new directions of thinking on a biological problem, making the mathematical model a useful tool accompanying biological reasoning and allowing to put in order diverse experimental observations.
However, to convert the methodological ideas presented in this paper into a working tool for experimental identification of the mechanisms of microRNAdependent protein translation inhibition, requires special efforts. Firstly, we need to construct a model which would include all known mechanisms of microRNA action. Secondly, realistic estimations on the parameter value intervals should be made. Thirdly, careful analysis of qualitatively different system behaviors should be performed and associated with the molecular mechanisms. Fourthly, a critical analysis of available quantitative information existing in the literature should be made. Lastly, the experimental protocols (sketched in the previous section) for measuring dynamical features such as the relaxation time should be developed. All these efforts makes a subject of a separate study which is an ongoing work.
Methods
Asymptotology and dynamical limitation theory for biochemical reaction networks
Most of mathematical models that really work are simplifications of the basic theoretical models and use in the backgrounds an assumption that some terms are big, and some other terms are small enough to neglect or almost neglect them. The closer consideration shows that such a simple separation on "small" and "big" terms should be used with precautions, and special culture was developed. The name "asymptotology" for this direction of science was proposed by [51] defined as "the art of handling applied mathematical systems in limiting cases".
In chemical kinetics three fundamental ideas were developed for model simplification: quasiequilibrium asymptotic (QE), quasi steadystate asymptotic (QSS) and the idea of limiting step.
In the IUPAC Compendium of Chemical Terminology (2007) one can find a definition of limiting step [52]: "A ratecontrolling (ratedetermining or ratelimiting) step in a reaction occurring by a composite reaction sequence is an elementary reaction the rate constant for which exerts a strong effect stronger than that of any other rate constant  on the overall rate."
Usually when people are talking about limiting step they expect significantly more: there exists a rate constant which exerts such a strong effect on the overall rate that the effect of all other rate constants together is significantly smaller. For the IUPAC Compendium definition a ratecontrolling step always exists, because among the control functions generically exists the biggest one. On the contrary, for the notion of limiting step that is used in practice, there exists a difference between systems with limiting step and systems without limiting step.
During XX century, the concept of the limiting step was revised several times. First simple idea of a "narrow place" (the least conductive step) could be applied without adaptation only to a simple cycle or a chain of irreversible steps that are of the first order (see Chap. 16 of the book [53] or [54] or the section "Dominant system for a simple irreversible catalytic cycle with limiting step" of this paper). When researchers try to apply this idea in more general situations they meet various difficulties.
Recently, we proposed a new theory of dynamic and static limitation in multiscale reaction networks [42,48]. This approach allows to find the simplest network which can substitute a multiscale reaction network such that the dynamics of the complex network can be approximated by the simpler one. Following the asymptotology terminology [55], we call this simple network the dominant system (DS). In the simplest cases, the dominant system is a subsystem of the original model. However, in the general case, it also includes new reactions with kinetic rates expressed through the parameters of the original model, and rates of some reactions are renormalized: hence, in general, the dominant system is not a subsystem of the original model.
The dominant systems can be used for direct computation of steady states and relaxation dynamics, especially when kinetic information is incomplete, for design of experiments and mining of experimental data, and could serve as a robust first approximation in perturbation theory or for preconditioning.
Dominant systems serve as correct generalization of the limiting step notion in the case of complex multiscale networks. They can be used to answer an important question: given a network model, which are its critical parameters? Many of the parameters of the initial model are no longer present in the dominant system: these parameters are noncritical. Parameters of dominant systems (critical parameters) indicate putative targets to change the behavior of the large network.
Most of reaction networks are nonlinear, it is nevertheless useful to have an efficient algorithm for solving linear problems. First, nonlinear systems often include linear subsystems, containing reactions that are (pseudo)monomolecular with respect to species internal to the subsystem (at most one internal species is reactant and at most one is product). Second, for binary reactions A + B → ..., if concentrations of species A and B (c_{A}, c_{B}) are well separated, say c_{A }>> c_{B }then we can consider this reaction as B → ... with rate constant proportional to c_{A }which is practically constant, because its relative changes are small in comparison to relative changes of c_{B}. We can assume that this condition is satisfied for all but a small fraction of genuinely nonlinear reactions (the set of nonlinear reactions changes in time but remains small). Under such an assumption, nonlinear behavior can be approximated as a sequence of such systems, followed one each other in a sequence of "phase transitions". In these transitions, the order relation between some of species concentrations changes. Some applications of this approach to systems biology are presented by [50]. The idea of controllable linearization "by excess" of some reagents is in the background of the efficient experimental technique of Temporal Analysis of Products (TAP), which allows to decipher detailed mechanisms of catalytic reactions [56].
Below we give some details on the approaches used in this paper to analyze the models by Nissan and Parker [41].
Notations
To define a chemical reaction network, we have to introduce:
• a list of components (species);
• a list of elementary reactions;
• a kinetic law of elementary reactions.
The list of components is just a list of symbols (labels) A_{1},...A_{n}. Each elementary reaction is represented by its stoichiometric equation
where s enumerates the elementary reactions, and the nonnegative integers α_{si}, β_{si }are the stoichiometric coefficients. A stoichiometric vector γ_{s }with coordinates. γ_{si }= β_{si} α_{si }is associated with each elementary reaction.
A nonnegative real extensive variable N_{i }≥ 0, amount of A_{i}, is associated with each component A_{i}. It measures "the number of particles of that species" (in particles, or in moles). The concentration of A_{i }is an intensive variable: c_{i }= N_{i}/V, where V is volume. In this paper we consider the volume (of cytoplasm) to be constant. Then the kinetic equations have the following form
where T is the temperature, w_{s }is the rate of the reaction s, v is the vector of external fluxes normalized to unite volume. It may be useful to represent external fluxes as elementary reactions by introduction of new component Ø together with income and outgoing reactions Ø → A_{i }and A_{i }→ Ø.
The most popular kinetic law of elementary reactions is the mass action law for perfect systems:
where k_{s }is a "kinetic constant" of the reaction s.
Quasi steadystate and quasiequilibrium asymptotics
Quasiequilibrium approximation uses the assumption that a group of reactions is much faster than other and goes fast to its equilibrium. We use below superscripts '^{f}' and '^{s}' to distinguish fast and slow reactions. A small parameter appears in the following form
To separate variables, we have to study the spaces of linear conservation law of the initial system (35) and of the fast subsystem
If they coincide, then the fast subsystem just dominates, and there is no fastslow separation for variables (all variables are either fast or constant). But if there exist additional linearly independent linear conservation laws for the fast system, then let us introduce new variables: linear functions b^{1}(c),...b^{n}(c), where b^{1}(c),...b^{m}(c) is the basis of the linear conservation laws for the initial system, and b^{1}(c),...b^{m+l}(c) is the basis of the linear conservation laws for the fast subsystem. Then b^{m+l+1}(c),...b^{n}(c) are fast variables, b^{m+1}(c),...b^{m+l}(c) are slow variables, and b^{1}(c),...b^{m}(c) are constant. The quasiequilibrium manifold is given by the equations and for small ε it serves as an approximation to a slow manifold.
The quasi steadystate (or pseudo steady state) assumption was invented in chemistry for description of systems with radicals or catalysts. In the most usual version the species are split in two groups with concentration vectors c^{s }("slow" or basic components) and c^{f }("fast intermediates"). For catalytic reactions there is additional balance for c^{f}, amount of catalyst, usually it is just a sum The amount of the fast intermediates is assumed much smaller than the amount of the basic components, but the reaction rates are of the same order, or even the same (both intermediates and slow components participate in the same reactions). This is the source of a small parameter in the system. Let us scale the concentrations c^{f }and c^{s }to the compatible amounts. After that, the fast and slow time appear and we could write c^{s }= W^{s}(c^{s}, c^{f}), , where ε is small parameter, and functions W^{s}, W^{f }are bounded and have bounded derivatives (are "of the same order"). We can apply the standard singular perturbation techniques. If dynamics of fast components under given values of slow concentrations is stable, then the slow attractive manifold exists, and its zero approximation is given by the system of equations W^{f}(c^{s}, c^{f}) = 0.
The QE approximation is also extremely popular and useful. It has simpler dynamical properties (respects thermodynamics, for example, and gives no critical effects in fast subsystems of closed systems).
Nevertheless, neither radicals in combustion, nor intermediates in catalytic kinetics are, in general, close to quasiequilibrium. They are just present in much smaller amount, and when this amount grows, then the QSS approximation fails.
The simplest demonstration of these two approximation gives the simple reaction: S + E ↔ SE → P + E with reaction rate constants and k_{2}. The only possible quasiequilibrium appears when the first equilibrium is fast: = κ^{±}/ε. The corresponding slow variable is C^{s }= c_{S }+ c_{SE}, b_{E }= c_{E }+ c_{SE }= const.
For the QE manifold we get a quadratic equation . This equation gives the explicit dependence C_{SE}(C^{s}), and the slow equation reads Ċ^{s }= k_{2}c_{SE}(C_{s}), C^{s }+ c_{P }= b_{S }= const.
For the QSS approximation of this reaction kinetics, under assumption b_{E}<< b_{S}, we have fast intermediates E and SE. For the QSS manifold there is a linear equation , which gives us the explicit expression for c_{SE}(c_{S}): (the standard MichaelisMenten formula). The slow kinetics reads . The difference between the QSS and the QE in this example is obvious.
The terminology is not rigorous, and often QSS is used for all singular perturbed systems, and QE is applied only for the thermodynamic exclusion of fast variables by the maximum entropy (or minimum of free energy, or extremum of another relevant thermodynamic function) principle (MaxEnt). This terminological convention may be convenient. Nevertheless, without any relation to terminology, the difference between these two types of introduction of a small parameter is huge. There exists plenty of generalizations of these approaches, which aim to construct a slow and (almost) invariant manifold, and to approximate fast motion as well. The following references can give a first impression about these methods: Method of Invariant Manifolds (MIM) ([57,58], Method of Invariant Grids (MIG), a discrete analogue of invariant manifolds ([59]), Computational Singular Perturbations (CSP) ([6062]) Intrinsic LowDimensional Manifolds (ILDM) by [63], developed further in series of works by [64]), methods based on the Lyapunov auxiliary theorem ([65]).
Multiscale monomolecular reaction networks
A monomolecular reaction is such that it has at most one reactant and at most one product. Let us consider a general network of monomolecular reactions. This network is represented as a directed graph (digraph) [66]: vertices correspond to components A_{i}, edges correspond to reactions A_{i }→ A_{j }with kinetic constants k_{ji }> 0. For each vertex, A_{i}, a positive real variable c_{i }(concentration) is defined.
"Pseudospecies" (labeled Ø) can be defined to collect all degraded products, and degradation reactions can be written as A_{i }→ Ø with constants k_{0i}. Production reactions can be represented as Ø → A_{i }with rates k_{i0}. The kinetic equation for the system is
or in vector form: ċ = K_{0 }+ Kc. Solution of this system can be reduced to a linear algebra problem: let us find all left (l^{i}) and right (r^{i}) eigenvectors of K, i.e.:
with the normalization (l^{i}, r^{i}) = δ_{ij}, where δ_{ij }is Kronecker's delta. Then the solution of (36) is
where c^{s }is the steady state of the system (36), i.e. when all = 0, and c(0) is the initial condition.
If all reaction constants k_{ij }would be known with precision then the eigenvalues and the eigenvectors of the kinetic matrix can be easily calculated by standard numerical techniques. Furthermore, Singular Value Decomposition (SVD) can be used for model reduction. But in systems biology models often one has only approximate or relative values of the constants (information on which constant is bigger or smaller than another one). Let us consider the simplest case: when all kinetic constants are very different (separated), i.e. for any two different pairs of indices I = (i, j), J = (i', j') we have either k_{I }>> k_{J }or k_{J }<< k_{I}. In this case we say that the system is hierarchical with timescales (inverses of constants k_{ij}, j ≠ 0) totally separated.
Linear network with totally separated constants can be represented as a digraph and a set of orders (integer numbers) associated to each arc (reaction). The lower the order, the more rapid is the reaction. It happens that in this case the special structure of the matrix K (originated from a reaction graph) allows us to exploit the strong relation between the dynamics (36) and the topological properties of the digraph. In this case, possible values of l_{i }k are 0, 1 and the possible values of r_{i }are 1, 0, 1 with high precision. In previous works, we provided an algorithm for finding nonzero components of l_{i}, r_{i}, based on the network topology and the constants ordering, which gives a good approximation to the problem solution [42,48,50].
Dominant system for a simple irreversible catalytic cycle with limiting step
A linear chain of reactions, A_{1 }→ A_{2 }→ ...A_{n}, with reaction rate constants k_{i }(for A_{i }→ A_{i+1}), gives the first example of limiting steps. Let the reaction rate constant k_{q }be the smallest one. Then we expect the following behavior of the reaction chain in time scale ≳1/k_{q}: all the components A_{1},...A_{q1 }transform fast into A_{q}, and all the components A_{q+1},...A_{n1 }transform fast into A_{n}, only two components, A_{q }and A_{n }are present (concentrations of other components are small), and the whole dynamics in this time scale can be represented by a single reaction A_{q }→ A_{n }with reaction rate constant k_{q}. This picture becomes more exact when k_{q }becomes smaller with respect to other constants.
The kinetic equation for the linear chain is
The coefficient matrix K of these equations is very simple. It has nonzero elements only on the main diagonal, and one position below. The eigenvalues of K are k_{i }(i = 1,...n  1) and 0. The left and right eigenvectors for 0 eigenvalue, l^{0 }and r^{0}, are:
all coordinates of l^{0 }are equal to 1, the only nonzero coordinate of r^{0 }is and we represent vectorcolumn r^{0 }in row.
The catalytic cycle is one of the most important substructures that we study in reaction networks. In the reduced form the catalytic cycle is a set of linear reactions:
Reduced form means that in reality some of these reaction are not monomolecular and include some other components (not from the list A_{1},... A_{n}). But in the study of the isolated cycle dynamics, concentrations of these components are taken as constant and are included into kinetic constants of the cycle linear reactions.
For the constant of elementary reaction A_{i }→ we use the simplified notation k_{i }because the product of this elementary reaction is known, it is A_{i+1 }for i < n and A_{1 }for i = n. The elementary reaction rate is w_{i }= k_{i}c_{i}, where c_{i }is the concentration of A_{i}. The kinetic equation is:
where by definition c_{0 }= c_{n}, k_{0 }= k_{n}, and w_{0 }= w_{n}. In the stationary state (ċ_{i }= 0), all the w_{i }are equal: w_{i }= w. This common rate w we call the cycle stationary rate, and
where b = ∑_{i}c_{i }is the conserved quantity for reactions in constant volume. Let one of the constants, k_{min}, be much smaller than others (let it be k_{min }= k_{n}):
In this case, in linear approximation
The simplest zero order approximation for the steady state gives
This is trivial: all the concentration is collected at the starting point of the "narrow place", but may be useful as an origin point for various approximation procedures.
So, the stationary rate of a cycle is determined by the smallest constant, k_{min}, if it is much smaller than the constants of all other reactions (43):
In that case we say that the cycle has a limiting step with constant k_{min}.
There is significant difference between the examples of limiting steps for the chain of reactions and for irreversible cycle. For the chain, the steady state does not depend on nonzero rate constants. It is just c_{n }= b, c_{1 }= c_{2 }=... = c_{n1 }= 0. The smallest rate constant k_{q }gives the smallest positive eigenvalue, the relaxation time is τ = 1/k_{q}. The corresponding approximation of eigenmode (right eigenvector) r^{1 }has coordinates: . This exactly corresponds to the statement that the whole dynamics in the time scale ≳1/k_{q }can be represented by a single reaction A_{q }→ A_{n }with reaction rate constant k_{q}. The left eigenvector for eigenvalue k_{q }has approximation l^{1 }with coordinates . This vector provides the almost exact lumping on time scale ≳1/k_{q}. Let us introduce a new variable c_{lump }= ∑_{i}l_{i}c_{i}, i.e. c_{lump }= c_{1 }+ c_{2 }+... + c_{q}. For the time scale ≳1/k_{q}we can write c_{lump }+ c_{n }≈ b, dc_{lump}/dt ≈ k_{q}c_{lump}, dc_{n}/dt ≈ k_{q}c_{lump}.
In the example of a cycle, we approximate the steady state, that is, the right eigenvector r^{0 }for zero eigenvalue (the left eigenvector is known and corresponds to the main linear balance b: ≡ 1). In the zeroorder approximation, this eigenvector has coordinates
If k_{n}/k_{i }is small for all i < n, then the kinetic behavior of the cycle is determined by a linear chain of n  1 reactions A_{1 }→ A_{2 }→ ...A_{n}, which we obtain after cutting the limiting step. The characteristic equation for an irreversible cycle, , tends to the characteristic equation for the linear chain, , when k_{n }→ 0.
The characteristic equation for a cycle with limiting step (k_{n}/k_{i }<< 1) has one simple zero eigenvalue that corresponds to the conservation law ∑c_{i }= b and n  1 nonzero eigenvalues
where δ_{i }→ 0 when .
A cycle with limiting step (41) has real eigenspectrum and demonstrates monotonic relaxation without damped oscillations. Of course, without limitation such oscillations could exist, for example, when all k_{i }≡ k > 0, (i = 1,...n).
The relaxation time of a stable linear system (41) is, by definition, τ = 1/min{Re(λ_{i})} (λ ≠ 0). For small k_{n}, τ ≈ 1/k_{τ}, k_{τ }= min{k_{i}}, (i = 1,...n  1). In other words, for a cycle with limiting step, k_{τ }is the second slowest rate constant: k_{min }<< k_{τ }≤ ....
Authors' contributions
AZ and AG have written the main body of the manuscript. NM, NN and AHB provided the critical review of miRNA mechanisms and contributed to writing the manuscript. AG and AZ developed the mathematical methodology for identifying dominant systems. AZ performed the analytical computations and numerical simulations. All authors have participated in discussing the results and model predictions. All authors read and approved the final manuscript.
Acknowledgements
We acknowledge support from Agence Nationale de la Recherche (project ANR08SYSC003 CALAMAR) and from the Projet Incitatif Collaboratif "Bioinformatics and Biostatistics of Cancer" at Institut Curie. AZ and EB are members of the team "Systems Biology of Cancer" Equipe labellis'ee par la Ligue Nationale Contre le Cancer. This work was supported by the European Commission Sixth Framework Programme Integrated Project SIROCCO contract number LSHGCT2006037900. We thank Vitaly Volpert and Laurence Calzone for inspiring and useful discussions.
References

Chekulaeva M, Filipowicz W: Mechanisms of miRNAmediated posttranscriptional regulation in animal cells.
Curr Opin Cell Biol 2009, 21:452460. PubMed Abstract  Publisher Full Text

Eulalio A, Huntzinger E, Izaurralde E: Getting to the root of miRNAmediated gene silencing.
Cell 2008, 132(1):914. PubMed Abstract  Publisher Full Text

Filipowicz W, Bhattacharyya S, Sonenberg N: Mechanisms of posttranscriptional regulation by microRNAs: are the answers in sight?
Nat Rev Genet 2008, 9(2):102114. PubMed Abstract  Publisher Full Text

Yekta S, Shih I, Bartel D: MicroRNAdirected cleavage of HOXB8 mRNA.
Science 2004, 304(5670):594596. PubMed Abstract  Publisher Full Text

Kim D, Saetrom P, Snove OJ, Rossi J: MicroRNAdirected transcriptional gene silencing in mammalian cells.
PNAS 2008, 105(42):162305. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Place R, Li L, Pookot D, Noonan E, Dahiya R: MicroRNA373 induces expression of genes with complementary promoter sequences.
PNAS 2008, 105(5):160813. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Orom U, Nielsen F, Lund A: MicroRNA10a binds the 5'UTR of ribosomal protein mRNAs and enhances their translation.
Mol Cell 2008, 30(4):46071. PubMed Abstract  Publisher Full Text

Vasudevan S, Tong Y, Steitz J: Switching from repression to activation: microRNAs can upregulate translation.
Science 2007, 318(5858):19314. PubMed Abstract  Publisher Full Text

Gu S, Jin L, Zhang F, Sarnow P, Kay M: Biological basis for restriction of microRNA targets to the 3' untranslated region in mammalian mRNAs.
Nat Struct Mol Biol 2009, 16(2):144150. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nottrott S, Simard M, Richter J: Human let7a miRNA blocks protein production on actively translating polyribosomes.
Nat Struct Mol Biol 2006, 13(12):110814. PubMed Abstract  Publisher Full Text

Petersen C, Bordeleau M, Pelletier J, Sharp P: Short RNAs repress translation after initiation in mammalian cells.
Mol Cell 2006, 21(4):53342. PubMed Abstract  Publisher Full Text

Olsen P, Ambros V: The lin4 regulatory RNA controls developmental timing in Caenorhabditis elegans by blocking LIN14 protein synthesis after the initiation of translation.
Dev Biol 1999, 216(2):67180. PubMed Abstract  Publisher Full Text

Maroney P, Yu Y, Fisher J, Nilsen T: Evidence that microRNAs are associated with translating messenger RNAs in human cells.
Nat Struct Mol Biol 2006, 13(12):11027. PubMed Abstract  Publisher Full Text

Baillat D, Shiekhattar R: Functional dissection of the human TNRC6 (GW182related) family of proteins.
Mol Cell Biol 2009, 29(15):41444155. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Karaa Z, Iacovoni J, Bastide A, Lacazette E, Touriol C, Prats H: The VEGF IRESes are differentially susceptible to translation inhibition by miR16.
RNA 2009, 15(2):24954. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lytle J, Yario T, Steitz J: Target mRNAs are repressed as efficiently by microRNAbinding sites in the 5' UTR as in the 3' UTR.
PNAS 2007, 104(23):966772. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bandres E, Cubedo E, Agirre X, Malumbres R, Zarate R, Ramirez N, Abajo A, Navarro A, Moreno I, Monzo M, GarciaFoncillas J: Identification by Realtime PCR of 13 mature microRNAs differentially expressed in colorectal cancer and nontumoral tissues.
Mol Cancer 2006, 5:29. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Pillai R, Bhattacharyya S, Artus C, Zoller T, Cougot N, Basyuk E, Bertrand E, Filipowicz W: Inhibition of translational initiation by Let7 MicroRNA in human cells.
Science 2005, 309(5740):15736. PubMed Abstract  Publisher Full Text

Bhattacharyya S, Habermacher R, Martine U, Closs E, Filipowicz W: Relief of microRNAmediated translational repression in human cells subjected to stress.
Cell 2006, 125(6):11111124. PubMed Abstract  Publisher Full Text

Kiriakidou M, Tan G, Lamprinaki S, De PlanellSaguer M, Nelson P, Mourelatos Z: An mRNA m7G cap bindinglike motif within human Ago2 represses translation.
Cell 2007, 129(6):114151. PubMed Abstract  Publisher Full Text

Humphreys D, Westman B, Martin D, Preiss T: MicroRNAs control translation initiation by inhibiting eukaryotic initiation factor 4E/cap and poly(A) tail function.
PNAS 2005, 102(47):169616. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mathonnet G, Fabian M, Svitkin Y, Parsyan A, Huck L, Murata T, Biffo S, Merrick W, Darzynkiewicz E, Pillai R, Filipowicz W, Duchaine T, Sonenberg N: MicroRNA inhibition of translation initiation in vitro by targeting the capbinding complex eIF4F.
Science 2007, 317:17641767. PubMed Abstract  Publisher Full Text

Thermann R, Hentze M: Drosophila miR2 induces pseudopolysomes and inhibits translation initiation.
Nature 2007, 447:875878. PubMed Abstract  Publisher Full Text

Wang B, Love T, Call M, Doench J, Novina C: Recapitulation of short RNAdirected translational gene silencing in vitro.
Mol Cell 2006, 22(4):53360. PubMed Abstract  Publisher Full Text

Wakiyama M, Takimoto K, Ohara O, Yokoyama S: Let7 microRNAmediated mRNA deadenylation and translational repression in a mammalian cellfree system.
Genes Dev 2007, 21(15):185762. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang B, Yanez A, Novina C: MicroRNArepressed mRNAs contain 40S but not 60S components.
PNAS 2008, 105(14):53438. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Clhendrimada T, Finn K, Ji X, Baillat D, Gregory R, Liebhaber S, Pasquinelli A, Shiekhattar R: MicroRNA silencing through RISC recruitment of eIF6.
Nature 2007, 447(7146):823828. PubMed Abstract  Publisher Full Text

Jackson R, Standart N: How do microRNAs regulate gene expression?
Sci STKE 2007, (367):re1. PubMed Abstract  Publisher Full Text

Kozak M: Faulty old ideas about translational regulation paved the way for current confusion about how microRNAs function.
Gene 2008, 423:108115. PubMed Abstract  Publisher Full Text

ValenciaSanchez M, Liu J, Hannon G, Parker R: Control of translation and mRNA degradation by miRNAs and siRNA.
Genes Dev 2006, 20:515524. PubMed Abstract  Publisher Full Text

Kong Y, Cannell I, de Moor C, Hill K, Garside P, Hamilton T, Meijer H, Dobbyn H, Stoneley M, Spriggs K, Willis A, Bushell M: The mechanism of microRNAmediated translation repression is determined by the promoter of the target gene.
PNAS 2008, 105(26):886671. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Basu A, Chowdhury D: Modeling Protein Synthesis from a Physicist's Perspective: A Toy Model.
Am J Physics 2007, 75:931937. Publisher Full Text

Dong J, Schmittmann B, Zia R: Towards a model for protein production rates.
J Stat Phys 2007, 128:2134. Publisher Full Text

Basu A, Chowdhury D: Traffic of interacting ribosomes: Effects of singlemachine mechanochemistry on protein synthesis.
Phys Rev E 2007., 75(1)
021902111

Gilchrist M, Wagner A: A model of protein translation including codon bias, nonsense errors, and ribosome recycling.
J Theor Biol 2006, 239:417434. PubMed Abstract  Publisher Full Text

Zouridis H, Hatzimanikatis V: A Model for Protein Translation: Polysome SelfOrganization Leads to Maximum Protein Synthesis Rates.
Biophysical Journal 2007, 92:717730. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chou T, Lakatos G: Clustered Bottlenecks in mRNA Translation and Protein Synthesis.
Phys Rev Lett 2004, 93:198101. PubMed Abstract  Publisher Full Text

Huang J, Morris Q, Frey B: Bayesian inference of microRNA targets from sequence and expression data.
J Comp Bio 2007, 14:550563. Publisher Full Text

Levine E, Zhang Z, Kuhlman T, Hwa T: Quantitative Characteristics of Gene Regulation by Small RNA.
PLoS Biol 2007, 5(9):e229. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Khanin R, Vinciotti V: Computational Modeling of PostTranscriptional Gene Regulation by MicroRNAs.
J Comp Biol 2008, 15(3):305316. Publisher Full Text

Nissan T, Parker R: Computational analysis of miRNAmediated repression of translation: Implications for models of translation initiation inhibition.
RNA 2008, 14:14801491. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gorban A, Radulescu O, Zinovyev A: Asymptotology of chemical reaction networks.
Chemical Engineering Science 2010, 65:23102324. PubMed Abstract  Publisher Full Text

Iacono M, Mignone F, Pesole G: uAUG and uORFs in human and rodent 5'untranslated mRNAs.
Gene 2005, 349:97105. PubMed Abstract  Publisher Full Text

Wu L, Fan J, Belasco J: MicroRNAs direct rapid deadenylation of mRNA.
PNAS 2006, 103(11):40349. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Selbach M, Schwanhausser B, Thierfelder N, Fang Z, Khanin R, Rajewsky N: Widespread changes in protein synthesis induced by microRNAs.
Nature 2006, 455(7209):5863. Publisher Full Text

George A, Tenenbaum S: MicroRNA modulation of RNAbinding protein regulatory elements.
RNA Biol 2006, 3:5759. PubMed Abstract

Jing Q, Huang S, Guth S, Zarubin T, Motoyama A, Chen FJ, abd Di Padova , Lin S, Gram H, Han J: Involvement of microRNA in AUrich elementmediated mRNA instability.
Cell 2005, 120:623634. PubMed Abstract  Publisher Full Text

Gorban AN, Radulescu O: Dynamic and static limitation in reaction networks, revisited.
Advances in Chemical Engineering 2008, 34:103173. Publisher Full Text

Eigen M: Immeasurably fast reactions. Nobel Lecture. In Nobel Lectures, Chemistry, 19631970. Amsterdam: Elsevier; 1972:170203.

Radulescu O, Gorban A, Zinovyev A, Lilienbaum A: Robust simplifications of multiscale biochemical networks.
BMC Systems Biology 2008, 2(1):86. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Kruskal M: Asymptotology. In Proceedings of Conference on Mathematical Models on Physical Sciences. Englewood Cliffs, HJ: PrenticeHall; 1963:1748.

Ratecontrolling step: IUPAC Compendium of Chemical Terminology [http://goldbook.iupac.org/R05139.html] webcite
2007.

Johnston HS: Gas phase reaction theory. New York: Roland Press; 1966.

Boyd RK: Some common oversimplifications in teaching chemical kinetics.
J Chem Educ 1978, 55:8489. Publisher Full Text

White RB: Asymptotic Analysis of Differential Equations. London: Imperial College Press & World Scientific; 2006.

Yablonsky G, Olea M, Marin G: Temporal Analysis of Products (TAP): Basic Principles, Applications and Theory.
Journal of Catalysis 2003, 216:120134. Publisher Full Text

Gorban AN, Karlin I: Invariant manifolds for physical and chemical kinetics, Volume 660 of Lect. Notes Phys. BerlinHeidelbergNew York: Springer; 2005.

Roussel M, Fraser S: On the geometry of transient relaxation.
J Chem Phys 1991, 94:71067111. Publisher Full Text

Gorban A, Karlin I, Zinovyev A: Invariant grids for reaction kinetics.
Physica A 2004, 333:106154. Publisher Full Text

Lam SH: Using CSP to Understand Complex Chemical Kinetics.
Combustion Science and Technology 1993, 89(5):375404. Publisher Full Text

Lam SH, Goussis DA: The CSP Method for Simplifying Kinetics.
International Journal of Chemical Kinetics 1994, 26:461486. Publisher Full Text

Zagaris A, Kaper H, Kaper T: Analysis of the computational singular perturbation reduction method for chemical kinetics.
J Nonlinear Sci 2004, 14:5991. Publisher Full Text

Maas U, Pope S: Simplifying chemical kinetics: intrinsic lowdimensional manifolds in composition space.
Combustion and Flame 1992, 88:239264. Publisher Full Text

Bykov V, Goldfarb I, Gol'dshtein V, Maas U: On a Modified Version of ILDM Approach: Asymptotical Analysis Based on Integral Manifolds Method.

Kazantzis N, Kravaris C: A New Model Reduction Method for Nonlinear Dynamical Systems using Singular PDE Theory. In Model Reduction and CoarseGraining Approaches for Multiscale Phenomena. Edited by Gorban A, Kazantzis N, Kevrekidis Y, Ottinger H, Theodoropoulos C. Springer; 2006:315. Publisher Full Text

Temkin O, Zeigarnik A, Bonchev D: Chemical Reaction Networks: A GraphTheoretical Approach. Boca Raton, FL: CRC Press; 1996.