Numerical modelling of microRNA-mediated mRNA decay identifies novel mechanism of microRNA controlled mRNA downregulation.

Post-transcriptional control of mRNA by micro-RNAs (miRNAs) represents an important mechanism of gene regulation. miRNAs act by binding to the 3' untranslated region (3'UTR) of an mRNA, affecting the stability and translation of the target mRNA. Here, we present a numerical model of miRNA-mediated mRNA downregulation and its application to analysis of temporal microarray data of HepG2 cells transfected with miRNA-124a. Using the model our analysis revealed a novel mechanism of mRNA accumulation control by miRNA, predicting that specific mRNAs are controlled in a digital, switch-like manner. Specifically, the contribution of miRNAs to mRNA degradation is switched from maximum to zero in a very short period of time. Such behaviour suggests a model of control in which mRNA is at a certain moment protected from binding of miRNA and further accumulates with a basal rate. Genes associated with this process were identified and parameters of the model for all miRNA-124a affected mRNAs were computed.

INTRODUCTION

Micro-RNAs (miRNAs) have become a major focus of research in molecular biology over the last decade. miRNAs are a family of small, non-coding RNAs, ∼21-nt long that regulate gene expression in a sequence specific manner (1–3). They were originally identified in Caenorhabditis elegans (4). Since then, hundreds of miRNAs have been identified in all metazoan genomes. miRNAs show diverse expression patterns and may play a regulatory role in various developmental or physiological processes. They are estimated to comprise up to 5% of animal genes, making them one of the most abundant classes of regulatory genes. Functional studies indicate that miRNAs participate in the regulation of almost every cellular process, and changes in their expression are observed in many human diseases, including cancer (5–8). It has been estimated that almost 30% of the protein coding genes of an organism are subject to miRNA-mediated control (9,10).

miRNAs control gene expression post-transcriptionally by regulating translation or stability of their target mRNAs. miRNA are processed from precursor molecules, which are transcribed from independent miRNA genes or from portions of introns of protein-coding RNA polymerase II transcripts (2). The pri-miRNAs are processed in two steps, catalysed by the RNA III-type endonucleases Drosha and Dicer, and mature miRNAs are then exported to the cytoplasm. There, they interact with mRNAs or the translational machinery, resulting in downregulation of the expression of their target mRNAs. In plants, miRNAs base-pair with mRNAs with nearly perfect complementarity, triggering endonucleolytic mRNA cleavage by an RNAi-like mechanism (11). In animal cells, they participate in the formation of micro-ribonucleoproteins (miRNAP), which bind to the 3′ UTR of a target mRNA and initiate deadenylation and decay of the target mRNA (12). Alternatively, the miRNAP can repress translation initiation at either the cap recognition stage or the 60 s subunit joining stage (13,14). Those mRNAs that are repressed by deadenylation or at the translation initiation step are moved to P-bodies for degradation or storage. Current knowledge about the mechanisms of miRNA-mediated mRNA downregulation has been discussed in several reviews (1–3,15–17).

Most of the computational efforts aiming to understand miRNA-mediated regulation have focused on identification of miRNA genes (18) and prediction of target mRNAs (9,19–21). Kinetic models of posttranscriptional regulation by small RNAs (22–24) or miRNAs have been proposed (10,24,25). Levine et al. (25) suggested a linear model of gene regulation by two classes of small RNAs in Escherichia coli and demonstrated that this regulation model follows a threshold linear response. The authors of this study introduced a two-step model, where the miRNA binds to the mRNA and initiates a secondary process that leads to degradation of the inactive mRNA. They showed that the mRNA levels, and the corresponding protein levels, could be modulated by target-specific parameters. Furthermore, they determined that the behaviour of this system could be modulated for some miRNA targets by global effectors. It has also been shown that regulation by small RNAs is advantageous when fast response is required, and that regulation by small RNAs is frequently used for additional control of gene expression (23). Xie et al. (10) analysed the role of miRNA in negative feedback regulation of gene expression. By altering the decay of mRNAs, miRNAs can control the dynamics of mRNA gene expression in a negative feedback loop. Recently, Khanin and Vincioti (24) introduced a model of miRNA-based control of mRNA degradation. The model extends the previous linear kinetics of Levine (22) to a more realistic saturative Michaelis–Menten type of kinetic. Parameters of this model are computed using an experimentally measured microarray time series of mRNA levels in response to transfection of miRNA. Parameters are computed using a sophisticated maximum likelihood optimization, developed originally for reconstruction of transcription factor activity (26). This work provides the first attempt to quantify the kinetics of free miRNA concomitantly with the effect on the downregulation of target mRNAs.

In this article, we introduce a model of the dynamics of post-transcriptional control of mRNA expression by binding of miRNA. The model considers both the kinetics of mRNA downregulation and the dynamics of the miRNA decay. Parameters of the model are estimated for a set of mRNAs downregulated by a miRNA from the miRNA-124a transfection experiment of Wang et al. (27). Simulation of the influence of the miRNA on mRNA degradation suggests a novel mechanism of miRNA control of target mRNA concentration in the cell.

MATERIALS AND METHODS

In modelling miRNA-mediated post-transcriptional silencing of mRNA, we start from the simplified situation shown in Figure 1. As mentioned in the introduction, miRNAs are synthesized from precursor pre-miRNAs, which are processed into a mature miRNA in two steps by Drosha and Dicer (1,2,17). For the purpose of modelling, without loss of generality, it can be assumed that the rate of synthesis in these steps is given entirely by the amount of miRNA gene transcript. Thus, the pri-miRNA processing steps can be combined into the rate constant driving the expression of the miRNA gene. Mature miRNA binds to target mRNA and the complex is then either degraded or stored in P-bodies, with the same net effect on total amount of free mRNA (1–3,15,17,28). The mechanism of mRNA downregulation was modelled by considering the concentrations of mRNA in three different forms, free mRNA, mRNA bound in the mRNA–miRNA complex, and mRNA localized to P-bodies (25). Repression of translation of the mRNA is caused by binding of the miRNA to its target in both of the latter cases, and it is not possible to distinguish between these two processes. Both processes lead to mRNA downregulation, thus both control mechanisms can be combined into one rate constant.

Figure 1.

Simplified scheme of miRNA-mediated silencing of a target mRNA. Both mRNAs and miRNAs are transcribed from their gene precursors. In the second step, miRNAs bind to complementary regions of their target mRNAs, allowing the complex to be transported to p-bodies and degraded. This process effectively reduces the amount of transcribed target gene mRNA as a function of miRNA local concentration. Synthesis of active miRNA molecules consists of several steps, which are omitted here, and are discussed further in the introduction.

In modelling of gene expression, it is commonly assumed that the rate of expression of a given target gene is controlled by the amount of transcription factor and its binding strength to the promoter region. The kinetic equation for expression of a given mRNA can be written as where represents the kinetics of mRNA expression dependent on the amounts of transcription factors z and the transcription-specific parameters b. Generally, the transcription rate depends on one or more transcription factors whose influence is weighted by the parameters vector b. In this work, transcription is assumed to be constant and independent of other genes or miRNAs (). This assumption is consistent with the design of the experiment of Wang et al. (27) analysed here, in which the cell line was transfected with miR-124a. The second term represents the degradation rate, controlled by the constant kd. In the absence of the miRNA, the degradation rate is constant, and together with the mRNA production rate, leads to a steady state of mRNA abundance. In the case when the miRNA binds to the mRNA (Figure 1), the degradation rate is no longer constant and becomes dependent on the amount of miRNA. The kd is therefore not a constant, but is dependent on the amount of the binding miRNA, s, and a set of parameters, a: where kd0 is a basal degradation rate constant and kd is the degradation rate given by the binding of the miRNA. This mass-balance equation has been used, in its linear form, for modelling of gene regulation by small RNAs (22) and in a more general saturative Michaelis–Menten type form (24).

It has been shown previously (29–31) that the rate of biological control processes often follows a sigmoidal function with an initial period where the reaction to an input signal is delayed, followed by linear-like increase growing until saturation. With this assumption, Equation (2) can be written as Parameter a is proportional to the strength of binding of the miRNA to its target and a can be considered a reaction delay parameter (32–34). Rewriting Equations (1)–(3) gives a differential equation for the rate of miRNA-mediated repression of mRNA levels: Steady state solution (dy/dt = 0) leads to for steady state levels, y, of the target mRNA.

If the miRNA concentration does not change with time (s = const), the expression level of the mRNA can be fixed to a constant value. As mentioned previously, at the beginning (t = 0) of the transfection experiment performed by Wang et al., the amount of mRNA is not yet affected by the miRNA. Therefore, y(t = 0) in Equation (1) can be fixed to y. With this assumption, Equations (1)–(3) can be solved for the target mRNA. Using Equation (6), the control parameters can be fitted independently for each target gene.

As with any other cellular macromolecules, the steady-state levels of miRNAs need to be properly controlled to ensure normal development and the miRNA molecules must be degraded in a controlled manner. Although much attention has been paid to the repression of mRNAs mediated by miRNA, much less attention has been given to the control of the miRNAs themselves. Ramachandran et al. (35) showed that in Arabidopsis, SDN gene products degrade mature miRNAs. The Wang et al. experiment analysed here shows, for all target genes selected as potentially controlled by the transfected miR-124a, an increase of the mRNA targets at the last time points. Such behaviour can be explained by a decrease in miRNA levels to an amount below the limit necessary to repress the target mRNAs. Such disproportion leads to restoration of the amounts of mRNA which are not sufficiently blocked by binding of such low amounts of miRNA. Therefore, the miRNA levels cannot be considered constant, but instead vary over time.

Considering a first-order type of degradation kinetics the level of miRNA is given by a mass balance equation where s represents the rate of production of the miRNA and a represents degradation with a rate constant a. Assuming that the miRNA is introduced into the cell line where it is absent (as in the experiment of Wang et al., where miR-124a was transfected to a cell line where it is normally not present), the miRNA temporal profile is given by where s represents the initial amount of miRNA at t = 0.

Similar to previous works (22,24), we ignore the possibility that the degradation of miRNAs is controlled non-linearly over time, and consider parameter a as a constant.

When the miRNA levels are time-dependent, Equations (1)–(3) cannot be solved analytically and their parameters must be estimated by fitting the model function to the experimental data.

RESULTS

Analysis of the temporal microarray dataset of mRNA changes after transfection of HepG2 cell with miRNA124a

The numerical model was applied to the analysis of a microarray dataset that recorded temporal changes in mRNA levels after transfection of HepG2 cells with miRNA124a (miR-124a) (27). The level of gene expression was measured for 54 675 genes at seven time points (4, 8, 16, 24, 32, 72 and 120 h after transfection). Negative control samples transfected as a reference were also recorded for each gene. The RMA-normalized data were downloaded from the NVBIA GEO database (GSE6207, http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?acc=GSE6207&targ=self&form=html&view=quick).

Potential targets of miR-124a in the Wang et al. set were selected using the TargetScan (9,19,20) web service (http://www.targetscan.org). Altogether, 773 conserved target records were identified. A target was defined as downregulated if the ratio between control and sample intensities was lower than 0.2 in at least one of the time points of the expression profile. This criterion identified 42 mRNA records, representing 25 individual genes (see Supplementary Table S1). This requirement is quite stringent and definitely did not identify all mRNA genes which are downregulated by miR-124a.

Parameters y0, kd0, kd1, a2 and a3 were computed by fitting the temporal profile of individual mRNAs by minimizing the sum of square differences between the experimental and the computed profiles. s, representing the initial amount of miR-124a immediately after its introduction to the cell, was arbitrarily set to the value of the most highly expressed mRNA in the set of the 42 selected profiles. Expression levels of mRNAs at t = 0 were extrapolated linearly from the values of the following time points.

Analysis of values of the individual parameters showed that a3 has a rather small variance (a3 = 0.016 ± 0.003 h−1 for α = 0.05) among all 42 analysed mRNA profiles. This observation is in good agreement with the experimental design analysed here. In the assumptions used in constructing the model, we excluded the possibility of further processing of the miRNA, implying that the miRNA degradation rate is constant and the same for all mRNAs. The constant value of a3 obtained for all selected mRNAs is in good agreement with this assumption. The obtained value of 0.016 is close to the predictions made by Khanin and Vinciotti (24) for the same experimental dataset (a3 = 0.024 or 0.02, respectively). Therefore, the value of a3 was fixed to 0.016 and the optimization process was run again for the remaining parameters for all selected mRNA profiles.

Temporal miRNA influence

For the 42 selected mRNA profiles, we analysed the temporal influence of the miRNA on mRNA downregulation given by the term

In most cases, this function follows a gradually decreasing sigmoid, as shown in Figure 2a, which was expected. For a specific set of genes, the profile of the miRNA influence was radically different (Figure 2b). The course of the miRNA influence resembled a threshold function. For the first period, the mRNA influence is maximal (equal to 1) and at a certain moment it drops to 0, with a very short transition period (Figure 2b). In order to investigate this property for a broader set, we selected all profiles (from the set of 773 potential target records) where the decrease in the profile, compared to control, was higher than the arbitrarily set value of 0.5. This procedure selected 339 mRNA profiles (219 individual genes), 11 of which were excluded due to high experimental error (marked red in Supplementary Table S1). Using the model, the parameters y, were optimized to fit the experimental mRNA profiles (solid line in Figure 2, see Supplementary Table S1) and the miRNA influence profile (dashed line in Figure 2), given by the Equation (9), was computed for all profiles in this set. The steepness of the drop of the miRNA influence profile can be characterized by the minimal value of its first derivative (for convenience, to get a positive value, the minima were multiplied by –1). In order to find other mRNA profiles of the type given in Figure 2b, we computed the first derivatives for all miRNA influence profiles and recorded their minima. The mRNAs were sorted according to increasing value of the sign-inverted minima. The plot of the sorted values showed a knee at the approximate value of 0.16 (Figure 3), dividing the whole set into two disjoint parts. The first part represented profiles given in Figure 2a, and the second part represented those with profiles characteristic for Figure 2b. The knee in the distribution diagram was quite sharp, indicating that there is a short transition between the two types of control, and that these two types are distinct from each other. This finding suggests that there exists either an analogue control, characterized by the curves in Figure 2a, or a digital switch, characterized by Figure 2b. Altogether, 98 (80 individual genes) of the switch-like type (Figure 2b) genes, out of 339 miR-124a downregulated mRNA profiles (219 individual genes), were identified (marked yellow in Supplementary Table S1).

Figure 2.

Two typical temporal profiles of mRNAs down-regulated by miR-124a. The mRNA profile is drawn with circles. The solid line represents a best fit to the data using the model. The influence of the miRNA given by Equation (9) (right axis) is indicated by a dashed line. Left axis—relative gene expression. (A) An ‘analogue’ mode of control represented by the profile of ATP6V0E1. (B) Digital-like mode of control represented by the profile of TEAD1.

Figure 3.

Sign inverted first derivatives of computed miRNA influence profile [Equation (9)] for selected mRNAs sorted according to the increasing value. Vertical line divides mRNAs to two types—analogue and digital like.

DISCUSSION

The bi-valued character of miRNA influence on mRNA accumulation (Figure 2b, dashed line) suggests the following model of control. In the first period, miRNA binds to mRNA and mediates its decay. At a certain moment, the mRNA becomes protected from binding of the miRNA and the miRNA influence drops to zero. Consequently, the mRNA decays with its basal rate only. The accumulation rate of the mRNA is then solely a difference between the rate of synthesis and the basal degradation rate. The mRNA profile shape associated with such behaviour (Figure 2b) is characteristic and different from the ‘ordinary’ shape given in Figure 2a, which is the expected model of miRNA-mediated mRNA downregulation. The difference between the individual profiles of this type is given only by the position of the transition phase. In most cases, the transition phase lies within the first 40 hours after transfection. The characteristic ‘V’ shape of the mRNA profile was also identified as an individual cluster of mRNA profiles by Khanin et al. (24) (as cluster 10 in Figure 2 of their article). Khanin et al. used a different algorithm for identification of potential miR-124a targets, and thus analysed a slightly different set of potential targets. In both cases, a set of profiles with the characteristic shape was identified, indicating that this group of mRNAs is conserved. Khanin et al. in their analysis used Michaelis–Menten type of kinetics and they did not analyse the influence of miRNA over time which showed here to be crucial in the identification of the ‘digital’ type of control found here. Their method focused more on identification of potential targets of miRNA than elucidation of the mechanism of their control.

Analysis of our data revealed that proteins that bind the 3′UTR of mRNAs (PTBP1, Quaking, Vigilin and ARPP-19) were highly represented in this group (Supplementary Table S2). Recently, it has been shown that association of RBPs to 3′UTRs of mRNA can prevent miRNA binding to these regions and protect the mRNA from miRNA-mediated repression (36,37). Thus, we hypothesized that for a specific group of mRNAs, the miRNA-dependent/miRNA-independent switch may be mediated by the binding of RBPs to the 3′UTR, preventing miRNA binding to those mRNAs. However, to act as a switch, RBP binding to the mRNA and mRNA protection should also display the switch-like behaviour, presumably through a sudden increase in the levels of RBPs. Importantly, we have identified four mRNAs encoding RBPs (Quaking, PTBP1, Vigilin and ARPP-19) that are targeted by miR-124a and display switch-like dependency on miR-124a, making them good candidates for RBPs that may protect mRNAs in a switch-like manner. Evidences for this hypothesis require additional experimental work which is out of scope of this article.

The miRNA-mediated degradation of target mRNAs assumes binding of the mature miRNA to the mRNA, which enhances the degradation rate or transport to P-bodies. The contribution of miRNA to the degradation of mRNA is generally understood to be proportional to the intracellular level of the mature miRNA. Thus, the contribution of miRNA to mRNA degradation can be viewed as graded or analogue when the output (mRNA degradation) is proportional to the input (miRNA level). In addition to the graded response, the data presented here suggest the existence of a switch-like, or digital, response in which the contribution of miRNA has a bi-valued character. In the first period, the miRNA mediates target mRNA decay. At a certain moment, however, the miRNA influence drops to zero. This switch-like behaviour may have profound consequences for mRNA control. Assuming a constant gene transcription rate, if the miRNA operates in digital mode ‘on’, the downregulation of mRNA is fully under the control of the miRNA. In contrast, if the miRNA operates in the ‘off’ mode, the accumulation of mRNA is simply controlled by the difference between the mRNA synthesis and basal decay rates, resulting in a rapid increase in the level of the mRNA. On the other hand, in the non-switch mode, the mRNA is under the control of decreasing amounts of miRNA, and consequently, the eventual increase of mRNA levels is much more delayed, if it occurs at all (compare Figure 2a and b). Thus, the switch–like behaviour may allow for rapid control of mRNA levels in response to an external signal, unresponsive to declining amounts of miRNA. In addition, this behaviour may counteract the ‘noisy’ expression of miRNAs themselves. Such behaviour suggests a model of control in which at the moment given by the minimum of the mRNA expression profile of the suspected genes the mRNA is protected from miRNA binding. mRNA is no longer recognized by the degradation machinery and can accumulate with a basal rate. The basal rate corresponds to the rate computed from the second part of the expression profile (Figure 2b). Analysis of the switch-like genes suggests that this mechanism could be associated with RNA-binding proteins, which could act as the RNA protecting agents. This hypothesis remains to be proven experimentally.

CONCLUSIONS

In summary, we developed a numerical model of miRNA-mediated downregulation of mRNA expression. Applying this model to a temporal series of microarray data identified a novel mechanism of miRNA-mediated control of mRNA expression. Our results suggest that miRNA-mediated degradation of mRNAs may proceed either by a graded response or in a digital-like manner, where the influence of the miRNA is stopped at a certain moment. We hypothesize that the ‘miRNA ON/miRNA OFF’ switch is mediated by RNA-binding proteins that bind to the 3′UTRs of target mRNAs and prevent miRNAs from binding these regions, halting the subsequent mRNA degradation. This mechanism would allow for rapid accumulation of transcripts coding for functionally related proteins, permitting cells to execute cell-specific or tissue-specific gene expression programs.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

FUNDING

The Czech Science Foundation, No. 310/07/1009 and 303/09/0475, Institutional Research Concept AV0Z50200510 and the grant of the EC Integrated Project ActinoGEN, LSHM-CT-2004-005224; JE Purkynje fellowship from the Academy of Sciences of the Czech Republic grant IAA500200716 to T.V. Funding for open access charge: Czech Science Foundation (No. 310/07/1009).

Conflict of interest statement. None declared.