Intrinsic Noise of microRNA-Regulated Genes and the ceRNA Hypothesis

MicroRNAs are small noncoding RNAs that regulate genes post-transciptionally by binding and degrading target eukaryotic mRNAs. We use a quantitative model to study gene regulation by inhibitory microRNAs and compare it to gene regulation by prokaryotic small non-coding RNAs (sRNAs). Our model uses a combination of analytic techniques as well as computational simulations to calculate the mean-expression and noise profiles of genes regulated by both microRNAs and sRNAs. We find that despite very different molecular machinery and modes of action (catalytic vs stoichiometric), the mean expression levels and noise profiles of microRNA-regulated genes are almost identical to genes regulated by prokaryotic sRNAs. This behavior is extremely robust and persists across a wide range of biologically relevant parameters. We extend our model to study crosstalk between multiple mRNAs that are regulated by a single microRNA and show that noise is a sensitive measure of microRNA-mediated interaction between mRNAs. We conclude by discussing possible experimental strategies for uncovering the microRNA-mRNA interactions and testing the competing endogenous RNA (ceRNA) hypothesis.


Introduction
MicroRNAs are short sequences of RNA (^22 base pairs) that post-transcriptionally regulate gene expression in eukaryotes by destabilizing target mRNAs [1,2]. Since their discovery almost two decades ago [3], there has been a steady increase in the number of discovered microRNAs. MicroRNAs play an important role in many biological processes, including animal development [4,5], tumor suppression [6,7], synaptic development [8,9], programmed cell death [10,11], and hematopoietic cell fate decisions [12,13]. In prokaryotes, an analogous role is played by an important class of small non-coding RNAs (antisense sRNAs) that also act posttranscriptionally to negatively regulate proteins. These antisense sRNAs can vary in size from tens to a few hundred nucleotides [14] and prevent translation by binding to the target mRNAs.
While both inhibitory microRNAs and sRNAs play similar functional roles, they act by very different mechanisms [15]. Eukaryotic MicroRNAs undergo extensive post-processing and are eventually incorporated into the RISC assembly [16,17]. The RISC complex containing the activated microRNA binds mRNAs and targets them for degradation. A single RISC complex molecule can degrade multiple mRNA molecules suggesting that microRNAs act catalytically to repress protein production. In contrast, both the mRNAs and sRNAs are degraded when bound to each other in a process that can happen spontaneously [18,19] or be mediated by extra machinery such as RNA chaperone Hfq [20][21][22]. This suggests that prokaryotic sRNAs, unlike their eukaryotic counterparts, act stoichiometrically on their targets. Stoichiometric regulation has been extensively studied theoretically and experimentally based on mass-action kinetics [23][24][25][26][27][28][29] and experimental protocols have been proposed for determination of different modeling parameters [30]. On the other hand, there exist relatively little similar work on understanding microRNAs and other catalytic non-coding RNAs [31][32][33]. Both theoretical calculations (see Supporting Information of [23]) and quantitative experiments [31] indicate that mean protein-expression of microRNA regulated genes follows a linear-threshold behavior similar to that of sRNAs in prokaryotes. In contrast, [33] argued that stoichiometric sRNAs are noisier than catalytic microRNAs. Presently, it is unclear how to reconcile these results and it raises the natural question of how the differing mechanisms employed by sRNAs and microRNAs affect the intrinsic noise profiles of regulated genes.
To answer this question, we used a generalized model of gene regulation by inhibitory non-coding RNAs to calculate the mean expression and intrinsic noise of regulated proteins. Our model is based on stochastic mass-action kinetics with tunable parameters that allow us to vary the mode of action from stoichiometric interactions such as those in prokaryotic sRNAs to highly catalytic interactions that mimic eukaryotic microRNAs. Our model is also applicable to plant microRNAs where microRNAs bind strongly to regulated mRNAs [34]. Contrary to [33], we show that in many parameter regimes the intrinsic noise properties of microRNAs and sRNAs are qualitatively similar. Finally, motivated by the competing endogenous RNA (ceRNA) hypothesis which suggests that microRNAs induce an extensive mRNA-interaction network, we extended our model to consider the case where a microRNA regulates multiple mRNAs. We calculate the intrinsic noise for these models and show that noise is an extremely sensitive measure of crosstalk between mRNAs. This suggests a new experimental strategy for uncovering microRNA-induced mRNA interactions. Recently we learned of three studies [35][36][37] completed simultaneously with our work, in which they model multiple mRNAs interacting with multiple microRNAs. Bosia et al. [35] studied robustness in a noisy model of ceRNA, while Figliuzzi et al. [36] and Ala et al. [37] both focus on steady state behavior. Figliuzzi et al. [36] studied sensitivity of ceRNA to transcription rate changes using analytical techniques whereas Ala et al. [37] combined analytical and bioinformatics results with experiments and showed cross-talk between ceRNAs. Our work complements and extends these other works to analyze the noise characteristics of microRNA interactions.

Model Description
Here we propose a model of gene regulation by non-coding RNAs based on mass-action kinetics [33]. A schematic of the model is shown in Figure 1A. Our model has four species: mRNA molecules denoted by m, functional non-coding RNAs denoted by s, mRNA-noncoding RNA complexes denoted by c, and the number of expressed proteins denoted by p. We note that s can represent either the concentration of prokaryotic small RNAs or the concentration of functional microRNAs found within the RISC complex. mRNA molecules are transcribed at the rates a m and translated at a rate a p . Free mRNAs degrade at a rate {1 m . a s represents the mean rate at which functional non-coding RNAs are formed. For prokaryotic sRNAs, this is simply the transcription rate of sRNAs. For microRNAs, this is an aggregate rate that accounts for the complicated kinetics of transcription and assembly into the RISC complex. mRNAs and noncoding RNAs form complexes c at a rate m and disassociate at rate c. Importantly, the complexes can also be degraded at a rate {1 c . This degradation can be actively regulated or conversely stem from dilution due to cell division.
Once mRNAs bind noncoding RNA and form a complex c, they can no longer be translated, resulting in decreased protein expression. To account for the differences between microRNAs and sRNAs, we have an additional parameter b which measures the ''recycling rate'' of noncoding RNA. When b& {1 c , c, the non-coding RNAs function catalytically with a single non-coding RNA molecule able to bind and degrade multiple mRNA molecules [38,39]. This is a good model for gene regulation by microRNA in eukaryotes. In contrast, when b% {1 c the noncoding RNAs function stoichiometrically. In particular, for prokaryotic sRNAs it is commonly believed that b~0 and no recycling of noncoding RNAs takes place. In plants microRNA's sequence almost perfectly matches the mRNA sequence [34]. Plant microRNAs can be modeled using intermediate values of the recycling rate b and disassociation rate c, although it is possible that they may behave stoichiometrically (b^0) as is the case with bacterial small RNAs. Finally, we note that the ratio of b to c is a measure of how much of the gene regulation happens through nucleolytic cleavage in contrast to translational repression.
We use two different approaches to mathematically model gene regulation by non-coding RNAs. We calculate both the mean expression levels as well as ''intrinsic noise'' due to stochasticity in the underlying biochemical reactions. First, we perform simple Monte-Carlo simulations for the reaction scheme outline above using the Gillespie algorithm [40]. While these simulations are exact, this method is computationally intensive making it difficult to systematically explore how noise properties depend on kinetic parameters. For this reason, we use a second approach based on linear noise approximation (LNA) [25,41,42]. The LNA approximates the exact master equations using Langevin equations. For very small particle numbers, this approximation breaks down even in determining the average particle numbers. Nonetheless, for medium and large particle numbers, the LNA is very adept at capturing qualitative behaviors of the system and we see a good agreement between LNA and exact simulations for our system for larger particle numbers and volumes. This allows us to derive analytical expressions for the noise and systematically explore how gene expression noise depends on system parameters.
with s,m,c,p being the concentration of noncoding RNA, mRNA, complex and protein respectively, and g s ,g m ,g c , and g p , being the noise in the birth-death processes of non-coding RNAs, mRNAs, complex, and proteins respectively, g m the binding noise, g c the unbinding noise, and g b the RNA recycling noise. The variance of these noise terms is given by where d ij is the Kronecker delta and overline represents timeaverage. For the noise terms we have: with n the time-averaged steady-state concentration of species n. Note we can convert between particle number and concentration by multiplying by appropriately chosen volumes (see Material and Methods). We modeled each interaction in Figure 1A as an independent Poisson process with a mean rate given by massaction kinetics [41]. Fluctuations in each interaction have been modeled by an independent Langevin term. Care must be taken to ensure the correct sign in the cross-correlations [42]. In the remainder of the paper, analytical results from the linear noise approximation are shown along with simulations using the Gillespie algorithm [40]. Both methods are in good agreement.
Finally, in what follows we will focus on protein mean and noise. However note that protein in our model is dependent on mRNA through a simple linear birth-death process. In this sense protein can be thought of as a readout of mRNA that goes through a lowpass filter set by {1 p . As a result the general behavior of protein and mRNA is very similar in our model and lead to similar qualitative results.

Mean Expression Levels
We begin by deriving the steady-state response of our system by setting the time-derivatives of the left hand side of equation (1) where we have defined the quantities Note that to derive equation (4) we have set the noise terms in equation (1) to zero. This effectively eliminates the contribution of cross-correlation terms into steady-state solutions (i.e. ms is approximated as m s). For very small particle numbers, where LNA breaks down, this approximation also breaks down. This can especially happen for small system volumes. However we will not study such parameter regimes in this paper (see Materials and Methods). Furthermore, we confine ourselves to the biologicallyrelevant parameter sets where the mRNA-microRNA association rate, m, is much larger than degradation rates of RNAs ( {1 m , {1 s ) and lambda is small (l%1). For these parameters, there is a clear linear-threshold behavior in mean protein concentration ( Figure  1B). Such behavior has been extensively studied in the context of non-coding RNAs where gene expression is often divided into three distinct regimes: a repressed regime, an expressing regime, and a crossover regime [23,25,31]. In the repressed regime (a m %wa s ), non-coding RNAs almost always bind mRNAs and prevent translation. In contrast, in the expressing regime (a m &wa s ) there are many more mRNAs than noncoding-RNAs, resulting in protein production that varies linearly with a m . Finally, there is a crossover between these two behaviors when a m^w a s .
The factor w multiplying a s renormalizes the transcription rate of the non-coding RNA to account for the fact that a single noncoding RNA can degrade multiple mRNAs. To see this note that is just the probability that a non-coding RNA is degraded when it is bound to an mRNA in a complex under the assumption complexes do not disassociate (c%b or c% {1 c ). Thus, we can think of w as the average number of mRNAs that a non-coding RNA will degrade before it is itself degraded. As expected, when b~0, w~1 and these results reduce to those derived for prokaryotic sRNAs [230 25]. Overall the experimentally observed linear threshold behavior of mean protein in quantitative single-cell experiments [31] is consistent with our theoretical calculation of the mean.
Comparing decay patterns of mRNAs repressed by microRNAs to free mRNAs provides a method to determine the existence and efficiency of microRNA binding. For example, Braun et al. [43] use immunoblotting to compare eukaryotic cells with no microRNA to cells transfected with microRNA transcripts and show a two-fold decrease in mRNA lifetime as a result of microRNA interaction. To determine if the general decay pattern of mRNAs in our model mimics their results, we studied the effect of microRNAs in the transient behavior of mRNAs. Figure 2 shows the results. To mimic the population averaging inherent to immunoblots, each curve is calculated by averaging 100 Gillespie simulations starting from the same initial conditions. Note that once the microRNA is introduced into the system (by setting m~2) the lifetime of mRNA dramatically decreases. Changing b from 0 (stoichiometric) to 10 (catalytic) further decreases the lifetime. The fluctuations in particle number close to zero are due to the finite ensemble size. In this regard our results qualitatively agree with experiments in [43] regardless of the regulatory regime studied. Furthermore we observe a fast initial decay rate followed by slower decay at later times instead of a single exponential decay. This is in agreement with the observation made in [43]. However we believe that beyond such qualitative agreement, a precise fit to the experimental data would require more complicated models. It is well known that the degradation of mRNAs includes several stages. In a recent work by Deneke et al. [44] it was shown that inclusion of multistep degradation can lead to piece-wise exponential decay of mRNAs. We speculate that such a multistep process in conjunction with the mechanism proposed here can lead to the correct quantitative decay patterns observed experimentally.

Intrinsic Noise
Gene regulation is inherently stochastic due to the small number of molecules present in the cell. Noise has important phenotypic consequences for a variety of biological phenomena [45] and for this reason it is important to characterize the intrinsic noise properties of non-coding RNAs. As is usual, we define the intrinsic noise as the variance in protein level divided by mean protein level squared, p{p ð Þ 2 p 2 , where overline represents steady state time averaging. This is a measure of relative protein fluctuations compared to their mean. Study of noise in bacterial sRNAs shows a peak in the crossover regime that emerges as a result of competition between mRNA and sRNA [25]. The switch-like behavior of the system, due to its linear-threshold nature, results in an amplification of any fluctuations in the mRNA level that is in excess of the mRNAs bound to noncoding-RNA. As a result, noise is increased in the crossover regime. We observe a similar qualitative behavior for noise in all parameter regimes of our model obtained by computing the solution to equation (1) at steady state. Figure 3 is a plot of protein noise as a function of mean protein concentration for catalytic and stoichiometric interactions, showing a peak at protein levels that correspond to the crossover regime in both cases. The height of the peak increases with the binding affinity m of a non-coding RNA for its target mRNA. On a whole, the noise profile of catalytic and non-catalytic genes is remarkably similar. There are however slight differences. The peak is slightly higher and occurs at a slightly larger mean protein level for catalytic non-coding RNAs. This suggests that the underlying reason for the noise peak is not the enzymatic mode of action of non-coding RNAs, but the fact that the number of mRNAs and the number of non-coding RNAs (appropriately normalized by w) are almost equal.
To better understand the effect of catalytic interaction on noise, we calculated protein noise versus a m for different values of b. Figure 4A shows the results for Gillespie simulation and linear noise approximation (the inset) as a function of a m wa s and b showing that the three-regime noise behavior is robust over a large range of b. Note that in all cases, the noise is peaked in the crossover regime. Aside from small discrepancies caused by the approximate nature of the analytical method as well as the finite statistical sample used in simulations, we see a good qualitative agreement between the two methods. To better visualize the robustness of noise to changes in b we have plotted the maximum of noise (at crossover regime) as a function of b in Figure 4B. Notice that the peak height initially increases as a function of b and reaches a plateau for large b upon entering the catalytic regime. However there are not any significant differences in the two interaction regimes. Finally as another test of the robustness of the model to parameter changes, we studied noise of a catalytic (i.e. eukaryotic microRNA, b& {1 c ,c) and stoichiometric interaction (i.e. prokaryotic sRNA, b~0) for the same level of mean protein produced. Figure 4C shows this comparison for noise in the repressed regime as a function of b with a m chosen such as to keep mean protein concentration constant. As can be seen, the noise decreases from its original value in the stoichiometric regime (b~0) to a slightly lower value as b is slightly increased (b^2) and any further increase in b does not affect noise dramatically. This again reveals that the system is robust to changes in b.
As a final test of the robustness of system on parameters, we studied protein mean and noise for a range of parameters using LNA. Figure 5 shows the results with first and second row corresponding to protein mean and noise respectively. Note that the linear-threshold behavior of mean and the peak at crossover regime is consistent for a wide range of parameters, showing that these qualitative features are independent of our choice of parameters, given that the number of particles are large enough for LNA to hold. The third row in figure 5 shows noise at the crossover regime (a m~w a s ) plotted for a wide range of paramters. We see a monotonic change in noise as different parameters are varied showing that changes in parameters can change the exact quantitative results but that the qualitative behavior of the intrinsic noise is similar for almost all biologically realistic choices of parameters.

Including Transcriptional Bursting
Experimental evidence suggests that mRNAs are often produced in bursts [46]. Transcriptional bursting represents another important source of stochasticity that was ignored in the analysis presented above. We can extend the model presented above to incorporate transcriptional bursting by considering the case where the genes encoding for mRNAs can be in two distinct states: an ''on'' state where mRNAs can be transcribed and an ''off '' state where transcription is not possible. For example, in eukaryotes the two states may correspond to whether the chromatin is condensed or not [47]. To model transcriptional bursting we replace the equation for mRNA production in Eq. (1) by the pair of equations where the probability that a gene is in the on state is denoted by g, k { (k z ) are the on(off) rate of the gene, and a on m is the maximum mRNA transcription rate. The average state of the gene is given by [25]. The gene noise g g is Gaussian white noise with mean zero and variance given by g g (t)g g (tz)~2k z gd(). The mRNA noise g m is now g m (t)g m (tz)~(a on m gz {1 m m)d(). All other equations remain the same as before.
In the analysis that follows, we assume that k { is fixed but k z can vary. This corresponds to the biological situation where a gene is regulated by a repressor that can turn the gene off. To compare noise for different values of b, we choose k z so that the mean protein levels remain constant. Furthermore, we concentrate on the repressed regime (see Figure 6). Here noise decreases slightly as b is increased which is very similar to the case without bursting as was shown in Figure 4B. This means that noise in the repressed regime is insensitive to bursting regardless of how catalytic the interaction is. This is true as long as the microRNA-mRNA binding rate is large compared to RNA lifetimes and mean protein levels exhibit the biologically observed linear-threshold behavior. However note that it is possible to design a very slow gene with the same mean value as here, by choosing small on and off rates, k { =k z , while keeping their ratio constant. In such cases, if the dynamics of the gene are chosen much slower than the dynamics of the rest of the network, effectively the system will switch between two different transcription rates corresponding to two different linear threshold behaviors. As a result such a system would behave as the superposition of two networks with their crossover regimes happening at different transcription rates, leading to a bimodal distribution of protein noise. It has been shown experimentally that in some biological systems transcription happens in this regime. For example Raj et al. [48] studied transcription through bursting in mammalian cells and discovered a bimodal protein distribution. Using stochastic modeling [49], they showed that their results agree with the notion of a slower gene activation process. In this work we have limited ourselves to faster gene dynamics where distributions are unimodal. However, as a future direction, it would be interesting to study the effect of noncoding RNAs on the time-scales of the process in different parameter regimes. If these time-scales remain small compared to gene activation/deactivation time-scale, protein noise will exhibit a bimodal distribution. We speculate that in this case instead of the three-regime behavior discussed here, the system will exhibit five (or four) regulatory regimes depending on the existence of (or lack of) leaky genes.

Asymptotic Formulas for Noise
Expressed. To gain further insight, we have derived asymptotic formulas for the noise in the repressed and expressing regime. In the expressing regime with large mRNA transcription, we define the small parameter E: l a m {wa s with l given by Eq. 8. In the expressing regime, E%1, and the steady-state expression levels of mRNA and non-coding RNAs take the form with q given by Eq. 7. Furthermore, the protein noise in this regime is identical to that of a transcriptionally regulated gene without post-transcriptional regulation [50] and is given by (dp where dp:p{p and b:a p m (see Materials and Methods) and it is assumed that mRNAs degrade much faster than proteins ( p & m ). In simple models of gene expression where mRNAs are produced stochastically with a Poisson rate a m , the quantity b is often called the 'burst size' and measures the average number of proteins made from a single mRNA molecule [51]. Our result suggest that in the expressing regime protein noise is independent of b and the noise profile of post-transcriptionally regulated genes is similar to those of genes regulated via transcription factors. [50].
Repressed. In the repressed regime, we now have E: l wa s {a m %1 and the average number of mRNA and noncoding RNA molecules is given by is the probability that a complex is degraded. The noise in this regime is given by where is a constant that is dependent on both b and c (see Materials and Methods). When b or {1 c &c it can be shown that ^1. Note that this condition includes the catalytic regime with b&c, {1 c as well as the stoichiometric regime b~0, {1 c &c (see Materials and Methods). Thus, we conclude that in the repressed regime, non-coding RNAs reduce noise. In particular, proteins are now produced from mRNAs in a smaller burst with typical size given by Eb%b. Since the burst size is just the average number of proteins made from an mRNA, b~a b m , we can equivalently interpret this results as changing the effective lifetime of the mRNA molecules from m to E m [23,25]. This interpretation is consistent with the stochastic simulations on mRNA lifetimes shown in Figure 2.

Scaling Behavior Near Crossover Regime
Our analysis show that the width of noise peak at the crossover regime is independent of recycling ratio b. To understand this behavior better, we studied the crossover regime for different parameter values using linear noise approximation (for comparison to simulations see Materials and Methods). Figure 7 shows the protein noise and protein mean as a function of a m wa s after rescaling both by their value at the crossover (a m~w a s ) for various recycling ratios. Notice that for c% {1 c these normalized plots of protein mean and noise show an approximate data collapse. As c is increased, this scaling behavior breaks down (Figure 7). The collapse of data for mean protein can be analytically derived given the fact that mean protein is only dependent on b,c, c through the   This scaling of the mean protein number can be better understood if we define x: a m wa s and define p(x) as the mean number of proteins corresponding to this value. Dividing this quantity by its value at x~1 (mean protein at crossover) and using equation (4) results in the following equation for the normalized mean protein: where n: qw is the probability a complex will disassociate. In the limit c%b, {1 c this parameter will be equal to 1 and the scaled mean protein level becomes independent of b causing the curves for different b to collapse on top of each other (Figure 7). It is also interesting to note that any other condition on the parameters that removes the dependence of qw on b will also have the same effect (e.g. c,b% {1 c ). Somewhat more surprisingly, near crossover the noise also shows an approximate data collapse. We suspect that this collapse is likely indicative of universality near the crossover between the repressed and expressing regimes.

mRNA Crosstalk and the ceRNA Hypothesis
Recently, the competing endogenous RNA (ceRNA) hypothesis proposed that microRNAs may play a crucial role in the cell in global gene regulation by inducing interactions between mRNA species [52]. The central mechanism underlying the ceRNA hypothesis is the idea that mRNA species may have interactions amongst themselves that are not direct but are instead indirect and mediated by competition and depletion of shared microRNA pools. The hypothesis is that these indirect mRNA interactions result in a biologically important mRNA network. However, the breadth and strength of microRNA induced interactions in eukaryotic genomes is still not well understood. For this reason, it is crucial to develop new strategies for measuring microRNA induced interactions between commonly regulated mRNAs. To this end, we asked whether the noise profile of regulated mRNAs could be used to uncover microRNA induced interactions. As a first step, we studied the simplified case where two different mRNA species are regulated by a single microRNA and compete over a common pool of microRNAs, and we focused on the effect of microRNA-induced crosstalk between mRNA species on the noise properties of regulated genes. The results presented here can be easily generalized to the case of many mRNAs interacting with many microRNAs.
A schematic of our simplified model is shown in Figure 8A. Two species of mRNAs are regulated by a common microRNA. Notice that the mRNAs do not directly interact in the model and all interactions are indirectly induced by the shared microRNA pool. We can represent this using a straight-forward generalization of the model considered earlier and variances reflecting the Poisson nature of interactions: The binding of microRNAs to mRNAs is controlled by the interaction rates m 1,2 . These rates reflect the binding affinity of microRNAs for the two mRNA species. We asked how protein noise and means change as we vary the transcription rates, a m1,2 , of the mRNAs. The remaining parameters are assumed to be the same for both mRNAs and from now on we have suppressed the indices on these parameters. MicroRNAs function catalytically and we focus on the parameter regime b& {1 c ,c. Figure 8B and C shows the mean protein levels, p 1 , and intrinsic noise of protein species 1 as a function of the transcription rates of the two mRNA genes, a m1,2 , for the case of equal binding affinities (m 1~m2 ). Notice there is a peak in the noise similar to the single-species case. Once again there is good agreement between simulation and analytic calculations based on Langevin noise and the Linear Noise Approximation (see Materials and Methods). We also examined the case where the mRNAs have different binding affinities for the microRNA, m 1~0 :2,m 2~2 . This results in an asymmetry in the noise peak but does not change the major qualitative features of our results (see Figure 8D and E) As in the single-mRNA species case, the behavior of the system can be divided into regimes depending on whether the combined normalized transcription rate of both mRNA species is bigger or smaller than the sRNA transcription rates. We find the crossover regime to happen at a m1 w 1 z a m2 w 2^a s with w i :1zb i ci (for derivation refer to Materials and Methods ). This region splits the transcription rate space (a m1 ,a m2 ) into an expressing regime, a m1 w 1 z a m2 w 2 * > a s , and a repressing regime, a m1 w 1 z a m1 w 2 =a s . Here, similar to the single species case we see a sharp peak in the noise at the crossover regime ( Figure 8F). Since the crossover regime depends on a linear combination of transcription rates, a m1 w 1 z a m2 w 2 , the existence of crosstalk can be determined by tracking the changes in crossover regime of one species as the transcription rate of another species is changed. In this regard measuring protein noise can be useful, since it has a peak that is easy to detect, whereas such a dramatic feature cannot be observed in the mean levels of proteins. This provides a tool for probing crosstalk between mRNAs in experimental setups if transcription rates can be tuned to access the crossover regime.
As an alternative experimental strategy for testing the ceRNA hypothesis, we studied mRNA decay in the context of ceRNA hypothesis, by tracking the mean mRNA number as a function of time before it reaches steady state. Figure 9A shows time-course of average mRNA number from species one as the transcription rate of mRNA species two is changed. Note that mRNA decays slower as the transcription rate of the competing mRNA is increased. This is due to access to a smaller pool of microRNAs, hence slower overall degradation. Figure 9B shows the same plots as interaction rate of the competing mRNA with microRNA is increased. Again competing over shared pool of microRNA results in a slower degradation. This suggests an alternative experimental approach for determining crosstalk through noncoding RNAs based on mRNA decay. We propose that such crosstalk can be detected by changing any paramaters of an mRNA that increases its competing strength for microRNAs (e.g. transcription rate, interaction rate) and probing changes in the decay rate of other mRNAs.

Discussion
In this work, we studied gene regulation by inhibitory noncoding RNAs. Whereas gene regulation by prokaryotic sRNAs has been extensively studied [23][24][25][26][27][28][29], there exist relatively few models of gene regulation by catalytic microRNAs [31][32][33]. Here, we used a simple kinetic model to study both the mean expression levels and intrinsic noise properties of post-transcriptional regulation by non-coding RNAs. Using a single parameter, our model interpolates between the stoichiometric behavior of prokaryotic sRNAs and the catalytic behavior characteristic of eukaryotic microRNAs. We found that both sRNAs and microRNAs exhibit a linear threshold behavior, with expressing and repressed regimes separated by a crossover regime. At the crossover, the mRNA transcription rate roughly equals the product of the non-coding-RNA transcription rate and the average number of mRNA molecules degraded by a single non-coding RNA molecule. In all cases, the intrinsic noise was smaller in the repressed regime and showed a sharp peak in the crossover regime. We found that for most parameter regimes, the intrinsic noise in the crossover regime shows an approximate data collapse, giving hints that there may be universal behavior near the transition from the repressed to expressing regime. We then generalized our calculations to study crosstalk between two mRNAs regulated by a single microRNA. We found that the intrinsic noise is an extremely sensitive measure of microRNA induced interactions between mRNAs.
Our results for the mean expression profile is consistent with recent experimental studies [31]. However, our conclusions about the intrinsic noise profiles of catalytic non-coding RNAs are different from similar work done by Hao et al. [33]. They concluded that the intrinsic noise profiles of catalytic and stochiometric interactions differed significantly. The reason for this discrepancy is that Hao et al. did not normalize the sRNA transcription rate a s by w. Consequently, they compared the crossover regions of sRNA regulated genes to repressed regions of microRNAs. As shown above, after making this normalization there is extensive data collapse of intrinsic noise profiles for both stoichiometric and catalytic genes.
One of the striking results of our calculation is the similarity between sRNA-regulated and microRNA-regulated genes. This similarity is somewhat surprising given that microRNAs and  Figure 1A. B,C) Gillespie results for protein mean (B) and noise (C) as a function of transcription rates of the two mRNAs with equal m 0 s. Parameters same as in Figure 1 with m 1~m2~2 . D,E) Gillespie results for protein mean (D) and noise (E) as a function of transcription rates of the two mRNAs with unequal m 0 s. Parameters same as in Figure 1 with m 1~0 :2,m 2~2 . Noise surface plots have been smoothed and interpolated for better visibillity. F) LNA results for noise of two nonidentical species as a function of normalized transcription rates (with w i :1zb i ci ). All the parameters are chosen to be distinct, i.e. a s~3 ,a p1~4 ,a p2~1 0, s~3 0, m1~1 0, m2~2 0, c1~1 , c2~5 , p1~3 0,b 1~1 0,b 2~1 ,c 1~1 ,c 2~1 0,m 1~0 :2,m 2~2 ,V~10. doi:10.1371/journal.pone.0072676.g008 sRNAs are found in different kingdoms (prokaryotes versus eukaryotes) and utilize distinct biochemistry and molecular machinery. Eukaryotic microRNAs use complicated nuclear machinery to export microRNA into cytoplasm and bring it to mature state by incorporating the RNA strand into the RISC complex. In contrast, prokaryotic sRNAs undergo relatively little post-processing and bind mRNAs spontaneously or via the chaperone protein Hfq. Given these extensive differences, the similarity between the expression characteristics of microRNAregulated and sRNA regulated genes are suggestive of convergent evolution.
Our calculations show that mRNAs regulated by catalytic noncoding RNAs have large peaks in the intrinsic noise. This differs significantly from results that would be derived from more traditional treatments of catalytic interactions based on the Michaelis-Menten or Hill equations. The underlying reason for this difference is that traditionally, the Michaelis-Menten equations are derived under the assumption that substrates of enzymes are in excess compared to the enzymes themselves. In contrast, here we are interested in the case where the number of enzymes (microRNAs) and number of substrates (mRNAs) are comparable. This accounts for the sharp peak in noise observed in the crossover regime in our model that is absent in traditional treatments of enzyme kinetics.
Our calculations also suggest a strategy for testing the ceRNA hypothesis [52], which posits that microRNAs induce extensive interactions between mRNA molecules. Our calculations suggest that protein noise may be a sensitive measure of the competition between the two species. Thus, it may be possible to uncover interactions between mRNA by measuring changes in the noise of regulated genes. We suspect that this will be true even when we generalize our calculations to consider the case where n mRNA species compete over the same pool of microRNAs. In this case, we hypothesize that there would still be a sharp peak in the intrinsic noise when the total appropriately normalized transcription rates of all mRNAs equals the transcription rate of microRNAs. In the future, it will be interesting to analyze these more complicated models in greater detail.
Finally we would like to note that the results discussed in this paper are based on a simplified version of the actual biochemical interactions with the goal of incorporating different mechanisms of gene regulation by noncoding RNAs into a simple mathematical framework. Adding further levels of complexity to the model can result in more realistic results, however the use of simple models can still provide insights into qualitative behavior of the system. To determine the extent to which our results hold under more realistic biological assumptions further theoretical and experimental work is necessary. The microRNAs discussed in this paper correspond to the RNA interference (RNAi) process and repress genes by binding to mRNAs. However, more recently, a small subset of microRNA processes, known as RNA activation (RNAa), has been discovered that can activate genes by binding to mRNAs [53][54][55]. The underlying mechanism for gene activation is not fully understood, but involves derepression of already repressed genes via microRNAs. In the future, it will be interesting to incorporate this biology into our model. In addition, mRNAs often undergo degradation in several stages. In this work, we model the degradation of mRNAs as a Poisson process with a single degradation constant. However it is well-known that the degradation of mRNAs is a complicated process often involving multiple steps [56,57]. Each additional step can be modeled by extra chemical interactions. These new steps do not affect the mean mRNA decay rates, and our results qualitatively match the experimental results [43]. However these new interactions have additional noise that in general cannot be summed together as a single Poisson process. As a result our assumption of a single Poisson process for degradation underestimates noise. In this regard our results set a lower bound on the noise.
Finally, we note that prokaryotic and eukaryotic cells are likely to be described by very different parameter sets within our model. In this paper we largely focused on changing the recycling ratio b as a measure of catalytic versus stoichiometric interactions among mRNAs and noncoding RNAs. In doing so we were keeping the rates of production and decay of proteins and mRNAs constant. However these rates can vary significantly between prokaryotes and eukaryotes [58][59][60][61][62]. Nevertheless, as shown in Figure 5, our qualitative results are insensitive to most of the detailed kinetic parameters of regulation. This suggests our main conclusions should hold despite the differences in gene regulation between prokaryotes and eukaryotes.

Gillespie Simulations and LNA
We can convert between concentration and particle number by multiplying by a volume, V . As is standard in the field, we have used notation where V has been dropped from all the equations in the paper. The dependence on volume can be easily recovered by dimensional analysis. All Gillespie simulations (including ceRNA simulations) were done at the volume (V~10). The LNA results are also dependent on volume and are calculated with V~10. As a result, throughout the paper, particle numbers from Gillepie simulations and concentrations from LNA have been accordingly scaled by volume to account for their correct dimensions. The choice of volume has been such that the problem is computationally and analytically tractable. For larger volumes the number of particles increases, leading to a decrease in fluctuations. As a result the approximate methods produce very accurate results. However the simulations become more computationally demanding. On the other hand, for smaller particle numbers LNA breaks down [63,64]. It has been shown that in this case the peak in noise at crossover shifts from its expected theoretical value [25], and even mean protein levels do not match with simulation. This discrepancy increases as volume (and hence particle number) is decreased. Hence the results presented in this paper will qualitatively hold for some smaller systems, but they cannot be generalized to arbitrarily small systems. (dp) 2~p 2 2 p za 2 where l 4 :{ {1 p . Now by use of partial fractions and integration we get (dp) 2~p 2 (a p mz {1 p p)z

Single Species Scaling Results
As noted in the main text and Figure 7 we see a scaling at the crossover regime using linear noise approximation. We see similar results using Gillespie algorithm as shown in Figure 10. The slight discrepancy between the two methods is partly due to the approximate nature of LNA and partly due to finite statistical sample size used in Gillespie. As a result the maximum in Gillespie algorithm is always slightly off from a m~w a s .

Single Species Asymptotic calculations
Expressing Regime. In the expressing regime with large mRNA transcription rate we demand E: l a m {wa s %1. Expanding m and s in terms of E we will get: To find the eigenvalues of the transfer matrix A we will expand it in terms of E: after some calculation we find the following closed form for eigenvalues of A: where Figure 10. Scaling at Crossover Regime (Simulation Results). Parameters same as Figure 1 with m~2. Data has been smoothed using a moving average method. A) Protein noise normalized to its value at a m~w a s (more explicitly (dp) 2 (p) 2 = (dp) 2 (p) 2 jam where l's are the eigenvalues of A, and B is the matrix of eigenvectors, according to: we saw a good agreement between these results and Gillespie simulations. The two methods showed at most a deviation of 30% from each other (see Figure 11 for more information). After some calculations we derive the following polynomial for sRNA concentration, s:

ceRNA Asymptotics
Furthermore, for ease of notation, in what follows we will use the following definitions: Expressing Regime. In this regime A T &G:max(B T , KB T ,B 1 B 2 ) and we can simplify the polynomial equation by defining E:G=A T and multiplying it by E while keeping coeffiecients to first order: with D: Note that GD is of order O(E 0 ) and does not require E expansion. The equation for s has the following asymptotic solutions: