Slowly Produced MicroRNAs Control Protein Levels*

Proteins are the primary agents of function in biological systems, and their levels are critical control elements, reflecting the interplay between transcription, translation, and protein degradation. Here, we consider the role of microRNAs (miRNAs) in the post-transcriptional regulation of protein synthesis. To determine their impact on protein concentration, we constructed a mechanistic model consisting of four state variables and nine kinetic parameters that account for transcript sequestration and degradation via miRNA-mRNA complex formation. Our dynamical model predicts that, even when present in low copy number, miRNAs can exert potent effects on protein concentration. Sensitivity analysis of the steady-state solution indicates that miRNA synthesis commonly acts to fine-tune protein concentrations. However, the same analysis shows that for a small subset of miRNA-mRNA pairs characterized by slowly produced miRNAs, the miRNA synthesis rate is the dominant control element. Our model equations provide a tool to evaluate the importance of particular miRNAs on their target proteins and promote the development of miRNA-based therapies that target proteins associated with cancer, inflammation, and metabolic disorders.

Proteins perform most of the processes that make life possible and death inevitable. Proteins are involved in providing the structure of the cell and its machinery, compartmentalization, signaling, metabolism, and control of gene expression. A critical control point is the expression level of a protein, which may be considered as output of a complex system of processes beginning with the gene, involving gene transcription, mRNA translation, and protein degradation. MicroRNAs (miRNAs) 2 are short 21-26-nucleotide strands of mRNA that post-transcriptionally regulate gene expression. miRNA transcripts bind to the 3Ј UTR of mRNA transcripts and suppress translation by either physically preventing reading of the transcript by the ribosome or triggering transcript degradation via decapping or endonuclease cleavage (1). Studies estimate that miRNAs may regulate the expression of 30 -90% of human genes, with single miRNAs targeting multiple genes and single genes being regulated by multiple miRNAs (2)(3)(4). miRNAs play an increasingly appreciated role in a diverse range of biological responses, including development, apoptosis, oncogenesis, immunological response, and neurodegeneration (5)(6)(7)(8). In particular, aberrant miRNA expression acts both to suppress tumorigenesis as well as to promote oncogenic activity (9). miRNA may be involved in more than translational suppression and transcript stability; they may also affect transcript processing in the nucleus (10). These noncoding RNAs may be considered as a bridge between the control networks for gene and protein expression. Although numerous models have been proposed to identify the regulatory effects of miRNA (11)(12)(13)(14)(15), we know of no case where sensitivity analysis was performed on such a model to make claims about the relative importance of model parameters and their corresponding biological contributions.
A model that predicts the effect of a particular miRNA transcript on a given protein-coding gene would be highly informative in establishing its role in pathologic states and enhancing their candidacy as therapeutic targets. Nissan and Parker (15) developed a deterministic model of translational suppression, focusing on the individual steps involved in transcript translation, such as time to recognize the start codon and subunit binding. They reasoned that all reactions are elementary, allowing for the assumption of mass action kinetics (15). Here, we use the mass action assumption to build a model consisting of a system of four ordinary differential equations that predicts the behavior of four state variables (the concentrations of mRNA, miRNA, miRNA-mRNA complexes, and protein) over time. While focusing on less-specific processes than Nissan and Parker's (15) model, our model has the benefit of taking nine kinetic parameters that can be experimentally measured (16,17): the rates of mRNA synthesis, mRNA degradation, miRNA synthesis, miRNA degradation, protein synthesis, protein degradation, miRNA-mRNA complex association, miRNA-mRNA complex dissociation, and miRNA-mRNA complex degradation. In our model, a given miRNA transcript and its target mRNA transcript are defined by a full set of these nine kinetic parameters.
Having both a complete set of kinetic parameters and the initial concentrations for the state variables allows us to predict the dynamics of the four state variables. Each kinetic parameter was sampled randomly from a wide range of possible rates, with each set of the nine kinetic parameters representing a theoretically possible miRNA-mRNA pair. Numerically solving the equations gives a prediction of species concentrations over time. All solutions generated by randomly sampled parameter sets result in similar long term behavior characterized by an eventual settling down to a steady state. The steady-state equation for each state variable can be expressed solely in terms of kinetic parameters (supplemental part A). Sensitivity analysis was performed using the steady-state solution and parameter sets corresponding to randomly sampled miRNA-gene pairs to determine the frequency with which a given parameter contributes the most to variations in protein concentration. Our model predicts that the concentrations of most proteins with a regulatory miRNA partner will be affected mainly by the rates of protein synthesis, protein degradation, and mRNA synthesis. However, we also predict the existence of a subset of pairs that are especially vulnerable to miRNA. Pairs in this subset are shown to span a wide range of physiologically relevant protein concentrations.

MATERIALS AND METHODS
Sensitivity analysis was performed using MATLAB (Release 2009b, MathWorks, Natick, MA) coupled with the SIMLAB environment (version 2.2, Joint Research Centre, European Commission). Numerical solutions of the dynamical equations were generated using ode15s, a MATLAB function that implements the variable order method for solving sets of stiff equations (36).

RESULTS
Model Equations-For a given miRNA transcript and its target gene, the state variables of our model are the concentrations of four different quantities: mRNA, miRNA, miRNA-mRNA complexes, and protein. Each of these quantities, with associated kinetic parameters, is described below by an ordinary differential equation. We assume that the reactions are elementary, allowing us to write our equations in terms of species concentrations and kinetic rates, a simplifying assumption that is often made in models involving binding. For the sake of simplicity, we focus on miRNA-gene pairs that are one-to-one, in which a single miRNA transcript interacts exclusively with a single gene. Additionally, we assume that the rates in our model are constant. Gene activation may vary, but this assumption is necessary for a first-generation model and is common in such models (15). The model equations are as follows, where brackets indicate concentration, and Greek letters indicate kinetic parameters. The concentrations and rate parameters are in terms of "copies" (i.e. proportional to the number of molecules of miRNA, mRNA, and protein and number of units of miRNA-mRNA complexes). Fig. 1 provides an overview of the processes for which Equations 1-4 account. Equation 1 describes the rate of change of mRNA transcript concentration, [X mRNA ]. Positive contributions come from the rate of mRNA synthesis, ␣ mRNA , and the rate at which miRNA-mRNA complexes dissociate into individual mRNA and miRNA transcripts, ␤ C . Negative contributions come from the rate constant associated with miRNA-mRNA complex formation, ␣ C , and the rate of mRNA transcript degradation, mRNA . Equation 2 describes the rate of change of miRNA transcript concentration, [X miRNA ]. Positive contributions come from the rate of miRNA synthesis, ␣ miRNA , and the rates at which miRNA-mRNA complexes dissociate, ␤ C , and degrade, C . Negative contributions come from the rate constant of complex synthesis, ␣ C , and the rate of miRNA degradation, miRNA . Equation 3 describes the rate of change of miRNA-mRNA complex concentration, [X C ]. Positive contributions come from the rate constant for complex formation, ␣ C , and negative contributions come from the rates of complex dissociation, ␤ C , and complex degradation, C . Equation 4 describes the rate of change of protein concentration, [X P ]. Protein concentration increases with the protein synthesis rate, ␣ P , and decreases with the protein degradation rate, P .
The nine kinetic parameters of the model are summarized in Table 1. In our analysis of the model, the parameter space (i.e. all possible sets of values for the nine parameters) is defined by considering each kinetic parameter in the interval 10 Ϫ3 min Ϫ1 to 10 2 min Ϫ1 (for simplicity we omit "copies" from the units in Table 1). This choice is motivated by available data on human proteins that are synthesized at rates spanning the range from 10 Ϫ2 min Ϫ1 to 10 min Ϫ1 (18,19) and is meant to include most physiologically relevant values. Experimental values for biologically similar kinetic parameters, measured in bacteria, fall within our range and have been used in other biochemical reaction models (20 -22). A given miRNAgene pair is characterized by a set of nine parameters in this space. Because no data is available on how real miRNA-gene FIGURE 1. Diagram of miRNA-directed protein down-regulation. For a given gene, DNA is transcribed into mRNA at a rate ␣ mRNA . mRNA transcripts are translated into protein at a rate ␣ P , whereas transcripts degrade at a rate mRNA and protein degrades at a rate P . Similarly, miRNA transcripts are synthesized at a rate ␣ miRNA and degraded at a rate miRNA . miRNA transcripts also bind to mRNA transcripts forming miRNA-mRNA complexes at a rate ␣ C , which can dissociate back into miRNA and mRNA transcripts at a rate ␤ C or degrade at a rate P . Complex degradation destroys the mRNA transcript, freeing the miRNA to bind to additional mRNA transcripts.
miRNAs Control Protein Levels FEBRUARY 11, 2011 • VOLUME 286 • NUMBER 6 pairs are distributed over the parameter space, we take randomly sampled sets of parameters to represent theoretically possible miRNA-gene pairs.
Model Dynamics-To investigate general properties of the model, we studied the time-dependent equations. Simulations were performed by independently sampling each kinetic parameter logarithmically over the range from 10 Ϫ3 min Ϫ1 to 10 2 min Ϫ1 . We sampled logarithmically to ensure that each order of magnitude was probed equally. For simplicity, we numerically solved Equations 1-4 by setting the initial concentration of each state variable to 100. Although the results are qualitatively similar for different thresholds used to define the time to approach the steady state, we choose to define the transient time for each state variable as the time taken for the rate of change of concentration to drop below 10 Ϫ2 min Ϫ1 and remain below this threshold for all subsequent times. Fig. 2 shows the resulting transient time calculated for 10 5 independent simulations. Despite the coupling established through Equations 1-4, these results illustrate that the transient time is generally different for different species, with protein and complex concentration trailing the convergence of the others for most miRNA-mRNA pairs. Incidentally, because the dynamical equations depend linearly on the kinetic parameters, increasing all parameters by a common factor has the effect of reducing the transient times by a factor 1/.
In the absence of perturbations, protein concentrations reach equilibrium or steady-state values. That is, a balance will be attained between protein synthesis and degradation, resulting in no net change in protein concentration. Setting Equations 1-4 equal to zero and solving them for the steadystate protein concentration yields a unique solution, which depends on all nine kinetic parameters involved in the dynamical equations. Notably, steady-state protein concentrations are determined entirely by the kinetic parameters, which is a property also exhibited by the steady-state solutions of the three other species (supplemental part A). Although the conditions for this solution to be a global attractor (i.e. a set of the state space toward which a dynamical system evolves over time) remain to be established, local stability was observed in 100% of 10 6 simulations performed by logarithmically sampling the kinetic parameters from 10 Ϫ3 to 10 2 . Details are provided in supplemental part B. Some preliminary insights into the relationship between the protein concentration [X P ]* and miRNA synthesis rate ␣ miRNA can be obtained by writing Equation 5 as a function of ␣ miRNA , absorbing the other kinetic parameters into two constants, and ⑀. This yields the following expression, and ⑀ ϭ mRNA miRNA ͑␤ C ϩ C ͒ ␣ C C (Eq. 8)   (Fig. 3A) and ⑀ (Fig. 3B) also affect the steady-state concentration. But note that, for given and ⑀, the derivative of [X P ]* with respect to ␣ miRNA is Ϫ/(⑀ ϩ ␣ miRNA ) 2 , indicating that the miRNA synthesis will exert stronger regulatory effect at small rates (this is also evident from Fig. 3 when this figure is considered in linear scale). This observation can be generalized through formal analysis of the steady-state solutions, which is carried out next to investigate the sensitivity of protein concentration to variations in each kinetic parameter of the model. Sensitivity Analysis (SA)-SA is a technique for determining the robustness of the predictions of a mathematical model to the variation of its parameters, but it can also be used to discover which parameters create the greatest change in model predictions (23). For our purposes, SA will tell us the relative contribution of each kinetic parameter to variations in protein concentration. Many methods have been developed to perform SA, and the ideal method for a particular model depends upon characteristics of the system, including the degree of nonlinearity and number of parameters. We performed SA first using a method based upon partial derivatives, which is the most amenable to interpretation, and then using a variance-based method to further validate our results.
The SA based on partial derivatives is implemented for a function of n input parameters f ϭ f(x 1 , x 2 , . . . , x n ) by calculating n quantities, The quantity S i can be interpreted as a sensitivity index. It represents the instantaneous rate of change of the function f relative to the ith parameter when the remaining nϪ1 parameters are fixed. When S i is larger than S j , function f is more sensitive to parameter x i than to parameter x j . Thus, a variation in x i results in a greater change in the magnitude of f than would the same variation in x j . Because perturbing x i can more efficiently impact f, we say that x i is "more important" than x j . This aspect of our method is illustrated in supplemental Fig. S1.
Calculating the partial derivatives of [X P ]* with respect to each of the nine parameters of the model, we obtain that the conditions for the sensitivity index S ␣miRNA associated with ␣ miRNA to be larger than all the others are given by the following inequalities, where C ϭ mRNA ϩ D, and D ϭ C ␣ C ␣ miRNA ( miRNA (␤ C ϩ C )) Ϫ1 . Equation 11 follows from the comparison between S ␣miRNA and S mRNA , and Equation 10 follows from the comparison between S ␣mRNA and all of the other seven sensitivity indexes. Now, noting that the kinetic parameters are positive and that D/CϽ1, we obtain that for the rate of miRNA synthesis to be the most important model parameter, the rate of miRNA synthesis must be smaller than seven of the eight other kinetic parameters (Equation 10), and the rate of miRNA degradation must be smaller than a function of the kinetic parameters associated with the mRNA-miRNA complex (Equation 11). Under the conditions defined by Equations 10 and 11, the only kinetic parameter not directly upper bounding ␣ miRNA is the rate of mRNA degradation, mRNA , which only appears through the term C. But according to Equation 10, the smaller the parameter mRNA , the larger ␣ miRNA can be relative to the mRNA synthesis rate and the kinetic parameters of the protein. Thus, the rate of miRNA synthesis is necessarily constrained by all other parameters  FEBRUARY 11, 2011 • VOLUME 286 • NUMBER 6 when the protein is most sensitive to its variations. It is interesting that the condition of being most sensitive to miRNA synthesis also imposes a constraint on the rate of miRNA degradation (Equation 11).

miRNAs Control Protein Levels
To determine the frequency with which the different kinetic parameters are most important within our parameter space, we logarithmically sampled all parameters from 10 Ϫ3 to 10 2 to generate 10 6 independent miRNA-mRNA pairs. As shown in Fig. 4, the rates of mRNA and protein synthesis along with the rate of protein degradation are each the most important parameter in approximately 22.5% of the cases. This suggests that in general miRNA fine-tunes protein concentrations, a conclusion that is in agreement with previous models (15,24) and previous experimental work (25)(26)(27). However, the rate of miRNA degradation is found to be the most important parameter in 12% of the cases and the rate of miRNA synthesis in 3.1% of the cases. The latter is very significant for the possible development of miRNA-based therapies. Fig. 5 shows the distribution of kinetic parameters and concentrations corresponding to these 3.1% cases in which ␣ miRNA is the most important parameter. In this figure, the distributions are represented by the median, box plots at 25-75% and whisker plots at 5-95% generated using the same set of simulations used to generate Fig. 4. The distribution of ␣ miRNA is biased toward small values compared with the other parameters, with the next smallest ones being mRNA and miRNA (Fig. 5A), as expected from the mathematical analysis above. The higher sensitivity to smaller miRNA synthesis rates shown in Fig. 5A is also consistent with the more pronounced dependence of protein concentration on smaller values of ␣ miRNA predicted in connection with Equation 6. On the other hand, the distributions of the parameters ␣ mRNA , ␣ C , C , ␣ P , and P are all on average larger than the corresponding ones when all miRNA-mRNA pairs are considered (cf. Fig. 5B). Fig. 5 (C and D) shows the statistics for the concentrations of mRNA, miRNA, complexes, and protein. Compared with the set of all miRNA-mRNA pairs, those for which ␣ miRNA is most important have smaller miRNA concentrations, as expected. But they also have larger concentrations of mRNA, miRNA-mRNA complexes, and protein. Importantly, for all four species, there is significant overlap between the box plots for all pairs and the pairs for which ␣ miRNA is the most important kinetic parameter, providing evidence that miRNA synthesis can be the dominant control factor for a wide range of physiologically relevant conditions.
We have further validated our results using a variancebased SA method (28,29). This method provides a systematic way of measuring the relative contribution of variations in multiple model input parameters to variations in model output parameters (i.e. model predictions), an otherwise unwieldy task for more than a few parameters (see supplemental part D, for an in-depth explanation). Importance frequencies obtained via the variance-based SA are summarized in supplemental Fig. S2. The general trends are similar to those obtained for the derivative-based method (Fig. 4), except for the fact that the frequencies with which the rates of protein and miRNA degradation are the most important parameters are higher and lower, respectively, for the variance-based method. Some differences are in fact expected due to the one-factor-at-atime paradox (30), which relates to the fact that allowing only one variable to vary at a time while fixing all others significantly reduces the fraction of the parameter space that can be explored. That is, the slight difference observed between the results in Fig.  5 and supplemental Fig. S2 is likely to be due to differences in the distribution of importance frequencies in the parameter sets not explored by the derivative-based method. Most importantly, although the derivative-based method is biased against the rate of protein degradation and in favor of miRNA degradation, the overall trends with which we are concerned remain the same for the variance-based method. In particular, both methods predict approximately the same fraction of miRNA-mRNA pairs to be most sensitive to variations in miRNA synthesis, thus providing further evidence that our results are robust. . Importance frequencies of the kinetic parameters. Each bar indicates the frequency with which the corresponding parameter has the largest sensitivity index. These statistics were obtained from 10 6 miRNA-mRNA pairs randomly generated by logarithmically sampling each parameter between 10 Ϫ3 and 10 2 min Ϫ1 .

DISCUSSION
Changes in miRNA levels have been linked to an expanding list of processes and diseases in both animals and plants. Because multiple miRNAs can affect a single mRNA and a single miRNA can target multiple mRNA transcripts, the network of possible interactions has a high degree of complexity. As protein levels are also subject to other conditions (e.g. transcription or protein degradation rates), which themselves are tightly regulated by a number of controllers, the cause and effect of a miRNA is extremely difficult to predict with a high level of confidence. Furthermore, little is known about what controls the biosynthesis and stability of miRNAs (31). A necessary first step to understand this complexity is the modeling of the impact on protein concentration resulting from pairwise interactions between miRNA and mRNA transcripts.
Here, we have developed a dynamical model of protein synthesis that predicts the effects of a miRNA transcript on the level of the protein it targets. Although previous models of miRNA describe modulation of protein levels at steady-state, ours is the first to predict the existence of a subset of miRNAgene pairs for which the rate of miRNA synthesis is the primary controller of protein levels.
Sensitivity analysis of the steady-state protein concentration revealed that, for most miRNA-gene pairs, protein concentration is most sensitive to perturbations in the rates of protein degradation, mRNA synthesis, and protein synthesis, whereas the rate of miRNA synthesis is less important. This suggests that miRNA most often acts to finely tune protein levels, in agreement with previous model predictions (15,24). Statistical analysis of important parameter sets for a given parameter provides a novel method for characterizing it. In simulations where all parameters were sampled logarithmically from 10 Ϫ3 min Ϫ1 to 10 2 min Ϫ1 , the miRNA synthesis rate, ␣ miRNA , was the most important parameter 3.1% of the time. Studying this subset gives insight into the conditions under which the rate of miRNA synthesis is important and thus when altering the rate would have profound effects on protein concentration. From Fig. 5, we learn that miRNA synthesis plays an influential role in regulating protein concentration when the rate of synthesis is slow and hence slight perturbations in the rate of synthesis may provide a way for the up-or down-regulation of the expression of particular proteins.
We thus posit that miRNA-mRNA pairs with slow miRNA synthesis rates could serve as particularly sensitive targets in therapies aimed at modifying the miRNA synthesis rate, where decreasing the rate can lead to a marked increase in steady-state protein levels and increasing the rate can lead to a marked decrease. Such therapies could be designed to mediate the effects of aberrantly expressed miRNAs. miRNA-155 is a good example of the oncogenic potential of a single miRNA, yet little is known about its transcriptional regulation. Its overexpression leading to aberrant protein concentrations has been shown to induce hematologic malignancies including B cell lymphoma and acute myeloid leukemia in humans and transgenic mice (32)(33)(34). By reducing the effective number of circulating miRNA transcripts, one could counteract the effects of overexpression. This can be done by decreasing the rate of miRNA synthesis, increasing the rate of miRNA degradation, or simply sequestering miRNA by forcing binding to another partner (e.g. a complementary antagomir). On the other hand, there are miRNAs, such as let-7, which serve as tumor suppressors (35), and their levels may be manipulated to moderate the expression of genes they target. Although we focus on altering protein level via the miRNA synthesis rate, our model can assess the sensitivity of protein to all kinetic parameters; thus, it can be used not only to determine when altering the rate of miRNA synthesis is effective, but also to indicate which other kinetic parameters could be targeted.
The significance of the predictions of our model is corroborated by the fact that our aggregate rates are possible to measure experimentally. These rates implicitly incorporate the effects of other regulatory factors that are not explicitly modeled. For example, several of the rates included explicitly in Nissan and Parker's (15) model are folded into our rates for miRNA synthesis and degradation. We expect that our model will stimulate greater investigation into the synthesis and degradation rates of miRNAs in higher eukaryotes. Validation of our model with such experimental data from mammalian systems will permit the precise prediction of the impact of the role of miRNAs in physiology and therapeutics.