Skip to main content


RNAs competing for microRNAs mutually influence their fluctuations in a highly non-linear microRNA-dependent manner in single cells

Article metrics



Distinct RNA species may compete for binding to microRNAs (miRNAs). This competition creates an indirect interaction between miRNA targets, which behave as miRNA sponges and eventually influence each other’s expression levels. Theoretical predictions suggest that not only the mean expression levels of targets but also the fluctuations around the means are coupled through miRNAs. This may result in striking effects on a broad range of cellular processes, such as cell differentiation and proliferation. Although several studies have reported the functional relevance of this mechanism of interaction, detailed experiments are lacking that study this phenomenon in controlled conditions by mimicking a physiological range.


We used an experimental design based on two bidirectional plasmids and flow cytometry measurements of cotransfected mammalian cells. We validated a stochastic gene interaction model that describes how mRNAs can influence each other’s fluctuations in a miRNA-dependent manner in single cells. We show that miRNA–target correlations eventually lead to either bimodal cell population distributions with high and low target expression states, or correlated fluctuations across targets when the pool of unbound targets and miRNAs are in near-equimolar concentration. We found that there is an optimal range of conditions for the onset of cross-regulation, which is compatible with 10–1000 copies of targets per cell.


Our results are summarized in a phase diagram for miRNA-mediated cross-regulation that links experimentally measured quantities and effective model parameters. This phase diagram can be applied to in vivo studies of RNAs that are in competition for miRNA binding.


MicroRNAs (miRNAs) are small non-coding post-transcriptional repressors of gene expression [1]. They exert important regulatory functions on both protein-coding and non-coding genes and are often involved in pivotal biological processes, like developmental biology or the molecular pathogenesis of several diseases [25]. It is commonly believed that miRNAs play central roles in conferring robustness to biological processes against environmental fluctuations [610]. The common assumption that at any given time one miRNA molecule can interact at most with one target mRNA [11] suggests a whole new layer of post-transcriptional cross-regulation, lately named the competing endogenous RNA (ceRNA) effect [12]. This theory proposes that the amount of a gene product may be tuned by varying the concentration of another transcript that shares with it the same miRNAs. Qualitative experiments based on observing induced variations in the level of transcripts show indeed that these could be coupled, due to the interaction with a common pool of miRNAs [1316]. The discovery that the miRNA–target interaction is compatible with a titration mechanism [17] supports the emergence of hypersensitivity regions [18, 19] where miRNA targets should be highly correlated and their relative stoichiometry tightly controlled [2022]. Many biochemical competition phenomena that are qualitatively like the one studied here have been extensively studied in the past [2325]. However, the relevance of competition at the post-transcriptional level is still largely debated [26]. Indeed, on the one hand, absolute quantification experiments in primary hepatocytes and liver suggested that the ceRNA effect is unlikely to affect significantly gene expression and metabolism [27]. On the other hand, differential susceptibility based on endogenous miRNA/target pool ratios provide a physiological context for target competition in vivo [28]. Moreover, recent studies on topics ranging from development [29] to cancer [3034] show how the competition for miRNA binding can actively regulate key biological processes. Crosstalk among mRNAs may, thus, be regulated depending on miRNA and mRNA relative abundances and may exhibit non-negligible target correlation profiles [20, 21].

Here, we experimentally explore these features and address the relevance of the relative mRNA–miRNA stoichiometric composition. In particular, we set out to validate a stochastic gene interaction model [20] that describes how mRNAs sharing common miRNA regulatory elements (MREs) can influence each other’s fluctuations in a miRNA-dependent way. Thus, we designed two bidirectional plasmids, each with a two-color fluorescent reporter system, enabling the simultaneous tracking of gene expression in the presence and absence of MREs. This allows us to quantify the correlations between the expression of the encoded proteins under different conditions. We found that there is an optimal range of parameters (in terms of effective transcription rates and miRNA interaction strengths) for which cross-regulation is maximal among miRNA targets. We show that such cross-regulation arises both at the level of mean protein concentrations (like [35]) and, for the first time to the best of our knowledge, at the level of fluctuations and correlations.

We show that the optimal cross-regulation regime is compatible with low numbers of mRNA molecules. In particular, the bulk quantification of mean exogenous transcripts per cell reveals that in our experiments, the crosstalk is highest in a physiological regime of order 10 to 1000 molecules per cell [36, 37]. Interestingly, in agreement with the model, we found that the same mechanism may induce bimodal population distributions with distinct high and low expression states of the targets.


Stochastic titration model for crosstalk

We propose a stochastic model for the miRNA-mediated target cross-regulation [17, 20, 23] (see Fig. 1 a). Through the formulation of a chemical master equation (see “Methods” and Additional file 1 for details of the model), the model describes the mean amount and fluctuations of two mRNAs r 1 and r 2, which are both targets of the same miRNA s. Both r 1 and r 2 can be translated into proteins (p 1 and p 2, respectively) only when not bound to the miRNA. The mRNA–miRNA complex can be degraded as a whole, while the miRNA can be recycled with probability 1−α. Since the two targets r 1 and r 2, and thus p 1 and p 2, are coupled through their common regulatory miRNA s, the pool of available mature miRNAs is the limiting factor in a system of potentially interacting targets.

Fig. 1

Model and predictions. a Sketch of the minimal model of miRNA–target interactions. One miRNA s and two targets r 1 and r 2 are independently transcribed with rates k s , \(k_{r_{1}}\), and \(k_{r_{2}}\), respectively. Each transcript can then degrade with rate g s , \(g_{r_{1}}\), or \(g_{r_{2}}\), respectively. Each miRNA s can interact with targets r 1 or r 2 with effective binding rates g 1 or g 2. α measures the probability of miRNA recycling. If not bound to a miRNA, targets r 1 and r 2 can be translated into proteins p 1 and p 2, respectively, which could then degrade with rates \(g_{p_{1}}\) and \(g_{p_{2}}\). bd Predictions from the stochastic model of interactions sketched in (a) as a function of p 0 (which is the constitutive value of p 1 when g 1 tends to 0) in terms of b the mean amount of p 1 free molecules, c the p 1 coefficient of variation \({CV}_{p_{1}}\), and d the Pearson correlation coefficient between p 1 and p 2. In (bd), the red curve is the reference curve for a given set of parameters while the red line identifies the threshold. Blue and green curves show how the red curve would move when increasing the interaction strength with the second target g 2 or the pool of miRNA via the miRNA transcription rate k s , respectively. e Schematic representation of the two bidirectional plasmids coding for the four fluorophores. miRNA microRNA, UTR untranslated region

A Gaussian approximation of the master equation [20] (see Additional file 1) allows us to evaluate mean values (〈x〉), noise (coefficient of variation CV x =σ x /〈x〉), and Pearson correlation coefficients (ρ x,y =(〈x y〉−〈x〉〈y〉)/σ x σ y ) for each molecular species x represented in Fig. 1 a (with x{r 1, r 2, p 1, p 2, s}); see Fig. 1 bd, respectively. The parameters g 1 and g 2, which are proportional to the miRNA–mRNA association rate, determine qualitatively the shapes of the functions generated by the model. In Fig. 1 bd, these curves are plotted as a function of protein constitutive expression p 0 (i.e., the value of p 1 or p 2 when g 1 or g 2 tends to 0). When one of these parameters tends to zero (say g 2), its corresponding target (r 2) is not interacting with the miRNA. The other target r 1 (and, thus, p 1) is repressed until a threshold level of r 0 (and, thus, p 0) is exceeded (Fig. 1 b) [17]. The threshold is established by miRNA regulation and its location can be adjusted by regulating both the pool of miRNAs (via the miRNA transcription rate k s ) and the pool of targets (via the target transcription rates \(k_{r_{1}}\) and \(k_{r_{2}}\)) [17, 20]. The increase of g 1 sharpens the transition between threshold and escape regimes. From the point of view of r 1, g 2 (proportional to the association constant of the second target) governs the concentration of free miRNA available within the cell. Increasing g 2 (while keeping all other parameters fixed) pushes the threshold to lower values of expression (lower r 0 and p 0) and globally increases r 1 (and p 1). r 2 behaves as a sponge for the miRNA, and increasing g 2 is equivalent to sponging away the miRNA available to target r 1. When all the miRNAs have been sponged away by r 2 (high value of g 2), then r 1 is not regulated anymore. At intermediate conditions in which miRNA is not completely sponged away by one of the targets, finely tuned cross-regulation between targets is possible. The mathematical model, thus, suggests experiments for testing this hypothesis and quantifying the crosstalk, modulated by g 1, g 2, and the amount of miRNA present in the cell.

Experimental set-up for unraveling cross-regulation

To investigate the predicted miRNA-mediated cross-regulation in single mammalian cells, we transfected two different two-color fluorescent reporters (sketched in Fig. 1 e) in the HEK 293 cell line. Both constructs consist of bidirectional promoters driving two genes whose products are fluorescent proteins. The first construct expresses the fluorescent proteins mCherry and enhanced yellow fluorescent protein (eYFP) [17], while the second construct expresses mCerulean and mKOrange. The 3 untranslated region (UTR) of both mCherry and mCerulean was engineered to contain a fixed number N of MREs for miR-20a (with N=0,1,4,7), a miRNA endogenously expressed by the HEK 293 cell line [38, 39] and related to cell proliferation.

mCherry and mCerulean are, therefore, proxies for the two targets in the model (p 1 and p 2). The 3 UTRs of eYFP and mKOrange were left unchanged in order to measure the transcriptional activity of the reporters in single cells (they are proxies for p 0). The constructs, thus, allow simultaneous monitoring of protein levels with (mCherry and mCerulean) and without (eYFP and mKOrange) miRNA regulation. The absence of a control on the number of plasmids per cell allows us to explore the variation of target transcription levels by simply sorting the cells on their eYFP or mKOrange fluorescence level.

For single-construct transfections, we observed a threshold effect [17]. Briefly, when no MREs are present, the levels of expression of mCherry (mCerulean) and eYFP (mKOrange) are proportional. In cells with one or more miR-20a site on mCherry (mCerulean), the mCherry (mCerulean) level does not increase until a threshold level of eYFP (mKOrange) is exceeded (see Additional file 1: Figure S1). This indicates that protein production is highly repressed below the threshold established by miRNA regulation and responds sensitively to target mRNA input close to it. Cotransfections of both constructs with different MRE numbers and measurements of fluorescence with a flow cytometer enabled the quantification of the crosstalk between mCherry and mCerulean as a function of N. Different combinations of MRE numbers mimic the variation of the model parameters g 1 and g 2. To capture the cross-regulation quantitatively, we measured the joint distributions of mCherry (p 1) and eYFP (p 0) levels given mCerulean (p 2) in single cells positive to the fluorophores. We, therefore, binned the data according to their eYFP levels and calculated the mCherry and mCerulean mean levels as well as their standard deviations and Pearson correlation coefficients in each eYFP bin. To show that our results are unbiased with respect to the cell line we used, we repeated the experiments in HeLa cells. These data are presented in Additional file 1.

The extent of cross-regulation is determined by relative miRNA/target stoichiometry

To quantify the crosstalk by modulating g 1 and g 2, we expressed different numbers of MREs on mCherry and mCerulean (N=1,4,7) and compared them when N=0. These cotransfections allowed us to follow the expression of one target (mCherry) while tuning the amount of free miRNA via the second target (mCerulean). As predicted by the model (Fig. 1 b), it is possible to identify two different effects: (i) the appearance of a threshold on mCherry while increasing the number N of MREs on its 3 UTR and keeping N=0 on mCerulean (Fig. 2 a and Additional file 1: Figure S2a) and (ii) a global increase of the mCherry mean fluorescence and a shift in the threshold while increasing N on mCerulean (Fig. 2 b and Additional file 1: Figure S2b). Threshold effects such those shown in Fig. 2 are a typical landmark for non-linear behavior. More specifically, the strength of the miRNA–target interaction dictates the departure from linear behavior. In the absence of an miRNA–target interaction, the mean number of transcripts (and proteins) miRNA-mediated regulation breaks this linear dependence by inducing highly non-linear threshold effects.

Fig. 2

Titration-induced threshold determines the optimal crosstalk. a, b mCherry mean fluorescence (a proxy for p 1 in the model, Fig. 1 b) is plotted against eYFP (a proxy for the constitutive expression p 0 in the model). Error bars are evaluated on the biological replicates. Continuous lines are model fits. The gray curves in (a) and (b) are the model prediction with the parameters fitted from the data and miRNA/target effective interaction strength g 1. The black arrow points to the model-predicted threshold. A threshold (or non-linear behavior) emerges when increasing mCherry MRE (a) while it disappears when increasing mCerulean MRE (b). The onset of the threshold is very close to the origin of the plot, indicating a relatively small amount of free miRNA. The intensity of crosstalk (measured in terms of fold-repression F with respect to the unregulated fluorophores) depends on the particular combination of MRE on both exogenous targets (ce). F is the ratio between the value of mCherry in the absence of miR-20a MREs and its value in the presence of MREs for each eYFP bin and for each N on mCerulean. Purple and cyan circles in legends represent the plasmids coding for the mCherry and mCerulean fluorophores. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

In Fig. 2 a, b, circles are data points with error bars over the experimental replicates while continuous lines are model fits (see Additional file 1 for details). The gray curves in Fig. 2 a and b are the model prediction with the parameters fitted from the data and g 1. Notice that the onset of the threshold is very close to the origin of the plot, indicating a relatively small amount of free miRNA. mCherry tends to the unregulated case (mCherry is linearly proportional to eYFP) on increasing the number of MREs on mCerulean. This result is well summarized by the fold-repression F between regulated and unregulated mCherry mean fluorescence (Fig. 2 ce). F is the ratio between the value of mCherry in the absence of miR-20a MREs and its value in the presence of MREs for each eYFP bin and for each N on mCerulean. Increasing the number of MREs on mCherry increases its repression, and F is highest when mCerulean has N=0 MREs while it tends to one on increasing the eYFP expression or the number of MREs on mCerulean. In particular, near the threshold, F shows a maximum whose value depends both on mCherry and mCerulean MREs. F could be indirectly considered as a measure of cross-regulation between the two targets.

These data show that the cross-regulation is maximal near the threshold and for intermediate levels of repression (in our case when mCerulean is between one and four MREs).

miRNA increase shifts the optimal cross-regulation region

To assess the cross-regulation dependence on the availability of miRNA, we transfected 100 nM of pre-miR for miR-20a together with the bidirectional constructs. In our model, this is equivalent to increasing the basal miRNA transcription rate k s . We analyzed the cases with N=4 for mCherry and N=0,1,4,7 for mCerulean. In agreement with the model predictions (see Fig. 1 b), we observed a shift of the threshold towards higher eYFP levels (Fig. 3 a) together with a global increase in the fold-repression (Fig. 3 b) and a resulting shift of the optimal crosstalk region towards a higher number of MREs. In Fig. 3 a, triangles and circles are data points for cotransfections with pre-miR20a and negative controls, respectively, while continuous lines are model fits (see Additional file 1 for details). The gray curve in Fig. 3 a is again the model prediction with the parameters fitted from the data and g 1.

Fig. 3

miRNA increase shifts the maximal crosstalk region. a mCherry mean fluorescence (a proxy for p 1 in the model, Fig. 1 b) is plotted against eYFP (a proxy for the constitutive expression p 0 in the model). Blue triangles and red circles are data from cotransfection with pre-miR20a and negative controls, respectively. Error bars are evaluated on the biological replicates. The gray curve is the model prediction with the parameters fitted from the data and miRNA/target effective interaction strength g 1. The black arrow points to the model-predicted threshold. According to the model, increasing the pool of available miRNAs (transfecting pre-miRNAs) shifts the threshold to higher constitutive expression values. b Different combinations of miR-20a MREs lead to different levels of fold-repression and crosstalk. Triangles and circles in the plot are data from transfections with pre-miR20a and negative controls, respectively. Purple and cyan circles in legends represent the plasmids coding for mCherry and mCerulean fluorophores, respectively. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

We then quantified by quantitative PCR the mean absolute number of exogenous targets in three subpopulations of cells (bulk measurements), sorted according to their eYFP intensity (low, medium, or high) both in the presence and absence of pre-miR for when N=4 on mCherry and N=1 on mCerulean (Table 1). We found that mCherry and mCerulean ranged from 40 to 400 and from 10 to 200 mean molecules per cell, respectively, without pre-miR and both from 10 to about 250 mean molecules in the presence of pre-miR (see Additional file 1 for details). Mature miR-20a, both endogenous and in cells transfected with 100 nM of its pre-miR, was quantified as well. We found about 1250 molecules per cell of mature miR-20a in the untransfected cells and 163 times more mature molecules in the pre-miR transfected cells (see Additional file 1 for details).

Table 1 Absolute quantification of exogenous targets

These data show that in our system, the maximal cross-regulation region is dependent on the relative number of both miRNAs and targets and it is compatible with a low number of mRNA molecules.

Titration induces increased cell-to-cell variability and bimodality

It is well known that the intrinsic noise of an unregulated gene product decreases when its expression level increases [40]. The effect of miRNA regulation introduces an extra source of noise (extrinsic noise). Our mathematical model predicts that, at fixed levels of expression, the total noise (intrinsic plus extrinsic) of a miRNA-regulated gene product (say \({CV}_{p_{1}}\)) should increase or decrease upon enhancing miRNA–target interaction strengths g 1 or g 2, respectively (see Fig. 1 c), compared to the unregulated case (when g 1→0). The increase of \({CV}_{p_{1}}\) is due to the coupling of the intrinsic noise of the target r 1 (and, thus, p 1) to the extrinsic noise of the miRNA s induced by the titration reactions. Increasing the interaction strength of the second target g 2 partially decouples r 1 and s, thus, reducing the noise \({CV}_{p_{1}}\). In particular, the model predicts the onset of a local maximum in the noise profile of the miRNA target versus its level of constitutive expression for high g 1 near the threshold. Indeed, near the threshold, also the coupling between the two targets r 1 and r 2 (and, thus, p 1 and p 2) becomes non-negligible (maximal cross-regulation) and contributes to the total noise. The local maximum is a vivid manifestation of the so-called retroactivity phenomenon, i.e., of how binding and titration can introduce correlations between intrinsic and extrinsic noise [41]. The intrinsic noise of one target is coupled to the extrinsic noise of the miRNA and in turn to the extrinsic noise of the other miRNA target.

Experimentally, we show that: (i) upon increasing N on mCherry (i.e., g 1), the total noise of mCherry globally increases as a function of eYFP (Fig. 4 a and Additional file 1: Figure S2c) and (ii) upon increasing N on mCerulean (i.e., g 2), the total noise of mCherry globally decreases (Fig. 4 b and Additional file 1: Figure S2d). For high levels of repression (high N on mCherry and low N on mCerulean), the mCherry CV eventually shows a local maximum near the threshold (Additional file 1: Figure S2c, d). A low level of noise indicates unimodal distributions while an increase in noise indicates increased cell-to-cell variability and may indicate bimodal population distributions with distinct high and low expression states [42]. We then checked if this was the case and found that bimodality on mCherry is present near the threshold for a high miRNA–target interaction (N=4,7 on mCherry and N=0,1 on mCerulean); see the histograms in Fig. 4 and Additional file 1: Figure S3. In particular, for N=7 on mCherry and N=0 on mCerulean, two very discernible phenotypes appear. This suggests the binary response is directly linked to the variability in the level of repression the miRNA exerts on the target [20]. Near the threshold, where the numbers of free target and miRNA molecules are small and similar, stochastic fluctuations become decisive for the cell fate. A small fluctuation in the number of miRNAs or targets will indeed produce cells with a highly repressed or unrepressed target product depending on the particular miRNA repression strength exerted on the target.

Fig. 4

Retroactivity increases cell-to-cell variability. a, b mCherry total noise, quantified by its coefficient of variation (CV, a proxy for \({CV}_{p_{1}}\) in Fig. 1 c), is plotted against eYFP (a proxy for the constitutive expression p 0 in the model). The black arrow identifies the model-predicted threshold shown in Fig. 2. Error bars are evaluated on the biological replicates. CV globally increases on increasing the number of mCherry MREs (a) while it decreases on increasing the number of mCerulean MREs (b). The competition between these two strengths results in lowering the noise even if the expected repression from the rough number of mCherry MREs is high. Histograms in the lower panels show mCherry data distributions for the shaded regions in (a) and (b). A strong miRNA target repression strength increases cell-to-cell variability with the eventual appearance of different phenotypes (bimodal distributions). Purple and cyan circles in legends represent the plasmids coding for mCherry and mCerulean fluorophores, respectively. a.u. arbitrary units, CV coefficient of variation, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

Our data show that miRNA–target titration reactions introduce non-trivial couplings between miRNA and targets (retroactivity) that possibly result in an increase in noise and bimodal cell population distributions near the threshold.

Titration-induced retroactivity causes target synchronization

The model predicts a maximum in the correlation between the two target products p 1 (mCherry) and p 2 (mCerulean) near the threshold (see Fig. 1 d). We investigated the strength of this prediction, distinguishing between correlations dependent on the experimental setting (mainly transient cotransfections and partial sharing of regulatory elements in the promoter) and correlations induced by the competition for miRNA binding, which can potentially lead to synchronized fluctuations. We, thus, defined the ratio of the Pearson correlation coefficients (the ratio of the Pearson coefficients between mCherry and mCerulean possessing different MREs for the same measure in the absence of MREs). We measured this ratio for eYFP below, around, and above the threshold (Fig. 5 ac, respectively), and observed that the competition for miRNA binding introduces correlations ranging from 4- to 12-fold higher than the basal level of correlation. The correlation between targets is a measure of the extrinsic noise component of \({CV}_{p_{1}}\) induced by the miRNA titration reactions and gives insights into the possible synchronization of the targets.

Fig. 5

Fold Pearson and p values. The Pearson ratio is measured for three different values of eYFP basal expression: below threshold (a), around threshold (b), and above threshold (c). p values are reported for each combination of miRNA MREs on the two plasmids. The regions inside the blue perimeters are statistically significant with p<0.01. As predicted by the model, the correlation is maximal around the threshold and could be even 12-fold higher than in the unregulated case. Blue-delimited areas are regions whose Pearson ratio (i.e., the ratio of the Pearson coefficients between mCherry and mCerulean possessing different MREs for the same measure in the absence of MREs) is statistically relevant with respect to the corresponding unregulated case. eYFP enhanced yellow fluorescent protein, miRNA microRNA, MRE miRNA regulatory element

Our results show that it is possible to have weakly or highly correlated targets for precise transcriptional programs. The regime of synchronized fluctuations, which is due to the titration-induced retroactivity, is determined by the number of MREs on both targets and is maximal around the threshold for intermediate miRNA repression strengths.

Interplay between transcriptional activity and miRNA–target interaction strength

Although many mRNAs have more than one MRE for a given miRNA, most of them have just one (see Additional file 1). However, even if there is only one MRE, these targets are typically expressed in multiple copies per cell and, thus, have the potential to titrate away the available miRNA molecules and crosstalk with other targets, if expressed at sufficiently high levels. This suggests testing to see if increasing the number of molecules for a reporter with only one MRE has the same effect on the target competitor as increasing the number of MREs on the same reporter.

The model predicts that either an increase (respectively, a decrease) of the transcription rate of target 2 (\(k_{r_{2}}\)) or of its interaction strength (g 2) causes the same qualitative effect on the system: a decrease (respectively, an increase) of the number of free miRNAs available to bind the first target (r 1). However, the average number of transcripts (〈r 1〉) functionally depends in a different manner on the parameters \(k_{r_{2}}\) and g 2. This implies that the two effects, albeit qualitatively similar, are not equivalent (see Fig. 6 ac).

Fig. 6

Interplay between transcriptional activity and miRNA–target interaction strength. The figure shows model predictions and experimental results obtained when investigating the effect on one target (say p 1) of the interplay between the second target (say r 2 and, thus, p 2) and the miRNA. The interplay between r 2 and miRNA is tuned both via the transcription rate \(k_{r_{2}}\) of r 2 and via the interaction strength g 2 between r 2 and the miRNA. p 1 is plotted against p 0 on a increasing the transcription rate \(k_{r_{2}}\) of r 2, b increasing the interaction strength g 2 between miRNA and r 2 when \(k_{r_{2}} > k_{s}\) (excess of targets), and c increasing the interaction strength g 2 between miRNA and r 2 when \(k_{r_{2}} < k_{s}\) (excess of miRNA). The model prediction for cases depicted in (a) and (b) are qualitatively very similar. d mCherry mean fluorescence (a proxy for p 1 in the model) is plotted against eYFP (a proxy for the constitutive expression p 0 in the model). The dashed black line corresponds to the unregulated case while the blue data points correspond to the reference case with four MREs on mCherry and one MRE on mCerulean. Either increasing the copy number of mCerulean (a proxy for \(k_{r_{2}}\) in the model), black data points, or the number of MREs on its sequence (a proxy for g 2 in the model), red data points, has the effect of decreasing the amount of miRNA available to target mCherry (which globally increases). e Fold-repression with respect to the unregulated case plotted against eYFP. Error bars are evaluated on the biological replicates. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, miRNA microRNA, MRE miRNA regulatory element

In particular, increasing \(k_{r_{2}}\) has the effect of increasing the basal number of available targets r 2 and simply shifts the p 1 threshold to lower p 0 values. That is, the p 1 curve is shifted towards the left and we see the p 1 curve approaching the unregulated case (see Fig. 6 a). On the other hand, increasing g 2 mimics an increase in the binding efficiency between the target r 2 and the miRNA and results in a decrease in the probability of the miRNA binding to r 1. It is, moreover, possible to define two different regimes, depending on the basal transcription rate of r 2. If the r 2 transcription rate is smaller than the miRNA transcription rate, then on increasing g 2, p 1 will tend to a limiting curve that is different from the unregulated case (see Fig. 6 c). That is, even for a high interaction between miRNA and r 2, there will always be an excess of miRNAs so that r 1 (and p 1) shows a threshold-like behavior. If instead the r 2 transcription rate is bigger than the miRNA transcription rate, then r 1 (and p 1) will tend to the unregulated case when increasing g 2 (see Fig. 6 b) due to an initial target surplus. This last case is qualitatively like increasing the r 2 transcription rate.

We experimentally tested this hypothesis as follows. For a fixed total amount of DNA transfected per cell plate (5 µg), (i) we transfected 1 µg of mCherry reporter with four MREs and 4 µg of mCerulean reporter with one MRE and (ii) we transfected 1 µg of mCherry reporter with four MREs, 1 µg of mCerulean reporter with four MREs, and 3 µg of empty backbone vectors. We then compared the results with the unregulated case (zero MREs on both reporters) and the case with 1 µg of mCherry with four MREs and 1 µg of mCerulean with one MRE. This experiment relies on the assumption that an increase in the quantity of the transfected reporter is a proxy for an increase in the transcription rate of the gene coded in the reporter (i.e., \(k_{r_{1,2}}\)), while an increase of the number of MREs (i.e., 1, 4, and 7 binding sites) corresponds to an increase of the miRNA–target interaction strength g 1,2.

Our results, shown in Fig. 6 d and e, suggest that increasing the number of MREs on a reporter molecule or increasing the number of reporter molecules with only one MRE have qualitatively similar, although quantitatively different, effects in terms of crosstalk (Fig. 6 e). In particular, they also suggest an excess of targets, since the two curves (red and black) in Fig. 6 d are similar. Since from our target quantification the number of molecules per cell, close to the threshold, is about 150 and the number of mature miR-20a’s is about 1000 molecules per cell, to be in excess of targets the other endogenous miR-20 targets should be active in sequestering the miRNA.

Discussion and conclusions

Previous studies pointed out the functional relevance of RNA competition for miRNA binding, thus addressing the potential role of RNAs in regulating the distribution of miRNA molecules on their targets [2934, 43, 44]. In this work, we show, both with stochastic modeling and single-cell experiments that validate the model, that RNAs competing for miRNAs influence their relative fluctuations as well and that this happens in a miRNA-dependent manner. Our results offer a detailed feature map for characterizing the post-transcriptional mRNA cross-regulation (see Fig. 7). Besides the general consistency with previous population-based qualitative results [16] and the agreement with a titration-based mechanism of miRNA–target interaction [17, 20, 21], the stochastic analysis allowed us to characterize curve trends for fluctuations and correlations of two targets of the same miRNA as a function of their expression level. The detailed picture points out that crosstalk between targets is quantitatively relevant only in conditions of intermediate miRNA repression and small amounts of target molecules (of order 10–1000), in agreement with a cell-population-based study by Bosson and coworkers [28].

Fig. 7

Phase diagram for mCherry (the target product p 1). The figure shows how the crosstalk between targets and bimodality on mCherry behave on varying the effective miRNA interaction strength and the mean numbers of target mRNA molecules. The effective miRNA interaction strength on target r 1 (and, thus, p 1) is measured theoretically through the ratio g 2/g 1 and experimentally with different combinations of miRNA binding sites on both synthetic constructs

We stress that the numbers we obtained must be considered more like orders of magnitude than exact numbers, since the way we estimated them, although correct, is indirect [17, 45]. These values are a lower bound on the actual number of molecules per cell since the RNA extraction yield is lower than 100%. Nonetheless, our numbers are compatible with those in [17] and qualitatively in agreement with what we see at the level of fluorescence (i.e., the amount of exogenous RNAs increases on increasing the fluorescence of the eYFP reporter). This finding is in contrast with that in [35] where, with a similar experimental design, the authors observed that the exogenous targets are expressed as the most expressed endogenous genes. Despite that both our study and [35] used the HEK 293 cell line, the evaluated endogenous miRNAs are different: miR-20a in our case and miR-21 in [35]. According to [38], miR-20a is the most expressed miRNA in HEK 293 whereas miR-21 is not even in the list of the 20 most expressed miRNAs in that cell line. This implies that if the relative stoichiometry between miRNA and targets matters with respect to the effectiveness of competition, the discrepancy with what was observed in [35] may be related to the endogenous amount of miR-21. Indeed, we showed that the effectiveness of the competition is limited to particular stoichiometric conditions and that these conditions are achieved near the threshold produced by the miRNA–target titration reactions. Although cross-regulation between targets is not zero even for a high level of target fluorescence (high mRNA expression), the region of maximal sensitivity of the system is limited.

The miRNA–target titration reactions introduce correlations among the intrinsic noise of one target and the extrinsic noise of the miRNA and in turn the extrinsic noise of the other miRNA target. The model predicts that this phenomenon of noise coupling, first introduced in [41] as a consequence of binding reactions, is large near the threshold, when the amount of free miRNA is comparable to the number of free targets. The idea that binding reactions introduce correlations is defined as retroactivity (see [46] for a review) and described at the mean field level in [47]. Our results are a vivid manifestation of this retroactivity phenomenon, which per se can drastically affect noise transmission in cellular networks and eventually impede a modular description of biochemical networks [41]. However, when the retroactivity is high, two or more targets may be highly cross-correlated, thus our results suggest that optimal levels of expression of genes and of miRNAs with respect to maximizing retroactivity may control the relative fluctuations of targets that have to interact or bind in complexes with a precise stoichiometry [22].

On the other hand, we found that strong miRNA repression together with low target crosstalk is sufficient to induce bimodality, i.e., the appearance of two distinct populations of cells with low and high target expression states. Bimodally expressed genes are found in different contexts, ranging from breast cancer [48] to immune cells [49], and usually each mode of the bimodal distribution corresponds to a different physiological condition (for example, a normal or a disease state). Our results suggest that if the bimodally expressed gene is a miRNA target, the system could be locked in one of the two states, changing the miRNA–target interaction strength through the expression of other competitors.

Our experimental strategy allowed us to study the interplay between MRE multiplicity and the effect of differentially expressing one of the two targets. The stochastic model predicts, and our results confirm, that the effect of increasing the transcription rate of one of the two targets reduces the share of free miRNAs available to repress the other target. Qualitatively, the same effect could be achieved by increasing the miRNA–target interaction strength of one of the two targets. Although from a purely mathematical point of view the dependence of the average number of proteins of the first target on the transcription rate and on the miRNA–target interaction strength of the second one are different, there are ranges of parameter values for which this distinction is not so sharp.

According to Sætrom and colleagues, two miRNA molecules can cooperate in the repression of their target when seed binding sites are 13–35 nucleotides apart [50], a phenomenon that has also been studied theoretically [51, 52]. In contrast, MREs very near to each other may result in the interference of miRNA binding. Experimentally, we modulated the interaction strength through the number of MREs on the constructs. Our constructs have bulged miR-20a sites separated by four-nucleotide-long spacers [17] and the seed binding sites are 18 nucleotides from each other, which should imply the cooperative miRNA repression of the targets. However, there is no significant difference in the repression of mCherry with N=4 and N=7 MREs if there are no MREs on mCerulean (Fig. 2 a) and there is little shift in the mCherry fluorescence when the number of MREs on mCerulean increased from one to four (Fig. 2 b). A possible key to interpreting our results is that upon increasing the number of MREs per construct what is really increasing is the probability of miRNA binding to the target and not the number of miRNAs simultaneously bound to a single target. However, as recently shown in [53], the very mechanism of the miRNA–target interaction is not yet completely understood, and we expect new findings in the near future that will put our mathematical modeling on more solid ground.

It is tempting to speculate that gene expression thresholding could be an important feature of cell fate decisions. However, while titrative regulatory mechanisms of miRNAs may easily switch whole gene networks on or off depending on the relative stoichiometry of miRNAs and targets, the scope of this mechanism has to be determined case by case. Indeed, theoretically, there is in principle no limit on the number of genes that can be cross-regulated by one gene or via the expression of a common miRNA [20] (see Additional file 1), given that it is just a matter of parameters to tune to the appropriate value (for example, the transcription rate of one target or miRNA). In physiological conditions, however, where there is no fuel turnover, the number of genes involved in the crosstalk may be limited and case specific.

Taken together, our results suggest as well that cross-regulation is relevant for molecules present in small amounts. Molecular species physiologically present in the order of 10–1000 molecules per cell, such as transcription factors or signaling molecules [37, 54], are more likely to be affected by cross-regulation. Although our experimental setting is artificial, it provides a deep exploration of the parameter space. A physiological system of miRNAs and targets could indeed experience only a small subset of the features so far described, rendering the characterization of crosstalk difficult. Without needing to over-express any endogenous gene, we were instead able to mimic the variation of relevant model parameters (miRNA and target transcription rates and miRNA–target interaction strength). Our phase diagram (Fig. 7) links quantitative measurements (effective miRNA repression and number of mRNA molecules) with model parameters, suggesting the possibility of moving around in the phenotype space tuning quantities such as the accessibility of binding sites or the affinity between miRNA and targets.

The functional importance of miRNA-mediated target cross-regulation has been shown in a number of cases, both diseased and physiological. For example, the pseudogene PTENP1 regulates the expression of the tumor suppressor PTEN in a miRNA-dependent manner, eventually modulating cancer cell growth [43, 55]. Our findings suggest the possibility that when PTENP1 crosstalks with PTEN, their corresponding RNAs are close to the maximum-crosstalk region of the phase diagram. Tuning the miRNA–PTEN interaction strength through the expression of the miRNA-competitor PTENP1 would then move the system closer to the bimodality region, with the two modes of the distribution being the cancer and the normal cells. At this point, the cell population can switch from one subpopulation to the other and the two competitors are highly coupled. Subsequent downregulation of PTENP1 would then lock the system in one of the two states and move it toward the low-crosstalk region of the phase diagram, characterized by a monomodal distribution of the cell population. Another relevant case involves the long non-coding RNA linc-MD1, which competes for binding miRNAs with the two transcription factors MAML1 and MEF2C. Thus, it regulates the differentiation of myoblasts in normal muscle development [44]. As a last example, the 3 UTR of CD44 competes for miRNA binding with the mRNA of CDC42, whose corresponding protein regulates the cell cycle [56]. We think that our findings give theoretical support and valuable tools for making significant progress in the understanding of these and other miRNA-mediated target cross-interactions.


A detailed description of the modelling and data analysis procedures is available in Additional file 1 (Supplementary Material).

Reporter plasmid construction

The set of fluorescent reporters coding for eYFP and mCherry was obtained from Addgene (#31463, #31464, #31465, and #31466, deposited by Phil Sharp Lab) and are the same as those used in [17]. The second set of fluorescent reporters were cloned into pBI-CMV1 (Clontech). A nuclear localization sequence (NLS) (ATGGGCCCTAAAAAGAAGCGTAAAGTC) was appended to mCerulean-N1 (Addgene #27795, deposited by Steven Vogel Lab [57]) by PCR and then inserted into the main vector with ClaI and BamHI. mKOrange-NLS (Addgene #37346, deposited by Connie Cepko Lab [58]) was cloned into the vector using EcoRI blunt and BamHI. miR-20a regulatory elements were appended to the 3 UTR of mCerulean with the same strategy applied in [17]: the N=1 bulged miR-20 binding site (TACCTGCACTCGCGCACTTTA) was appended by PCR and for both constructs, CCGG spacers separate subsequent miR-20a regulatory elements.

Transient transfections

We performed two different methods of transfection, with Lipofectamine (data in Figs. 2, 3 and 4 in the main text and in Additional file 1: Figure S3) and with CaCl2 (data in Additional file 1: Figures S1 and S2). Our results are independent of the method of transfection.

Lipofectamine transfection method HEK 293 TeT-Off cells (Clontech) below passage 6 were plated in G418 (Gibco) 200 µg/ml media in six-well dishes the day before transfection. Reporter plasmids were transfected with Lipofectamine 2000 (Invitrogen) following the manufacturer’s specifications. miR-20a pre and negative controls (Ambion) were cotransfected at the indicated concentrations. The media was changed 24 h after transfection. Assays were performed 48 h after transfection.

CaCl 2 transfection method HEK 293 TeT-Off cells (Clontech) below passage 6 were plated in 100×20 mm (Falcon BD) dishes the day before transfection. The cells were transfected using CaCl2 protocols [59]. The media was changed 24 h after transfection. Assays were performed 48 h after transfection.

Flow cytometry

Cells were harvested 48 h after transfection (cell confluency 90 %) and run on a CyanADP (Beckman Coulter) flow cytometer. For each sample, at least 0.5×106 cells were acquired. The raw FACS data were analyzed with the Summit3.1 software (Beckman Coulter) to gate cells according to their forward and side scatter profiles and to define the intensity of fluorescent signals emitted by the four reporters in each cell. These values were normalized for background fluorescence by subtracting the mean plus two standard deviations of the fluorescent signal measured in the untransfected control cells. The data were then binned according to their eYFP values.

Data processing

For each cell, we collected the raw fluorescent intensities and then we subtracted the background fluorescence levels estimated from non-transfected cell reads. The background-corrected data were then binned into consecutive equally spaced intervals of eYFP intensities. Each bin contains a subpopulation of the order of 104 cells. Eventually, each bin is characterized by the mean fluorescent intensity, by the coefficient of variation for each fluorophore, and by the Pearson correlation coefficient between mCherry and mCerulean (computed for each subpopulation). The same procedure was then applied to the three different biological replicates. The final values we report in data plots (Figs. 2 a, b, 3 a and 4 a, b in the main text and Additional file 1: Figures S1 and S2) are the mean of these measurements over the biological replicates for each bin, and the error bar is their standard deviation. It is probably worth stressing that the error bar computed in this way is a proxy for the fluctuation of the mean over the biological replicates rather than the data dispersion, the latter being a factor \(\sqrt {N_{\text {cib}}}\) larger than the former, where N cib is the number of cells in a bin (see Additional file 1: Figure S4). The error bars for Figs. 2 ce and 3 b are instead computed with a standard jackknife procedure.

Fluorescence-activated cell sorting

Cells were transfected with the N=4 eYFP-mCherry and N=1 mkOrange-mCerulean reporters with and without pre-miR-20a 100 nM (Ambion). 48 h after transfection, three cell populations were sorted according to their eYFP fluorescence (low, medium, or high YFP expression) using a BD FACS Aria III (Becton Dickinson) cell sorter. Cell pellets were washed and snap-frozen before RNA isolation.

Empirical observables and Pearson correlation coefficient ratio

We defined the empirical average of a given observable O over an ensemble of cells as \(\langle O \rangle = \sum _{i\in \text {cell}} O_{i}/N_{\text {cell}}\). The Pearson ratio is defined as the ratio of the Pearson correlation coefficient (ρ x,y =(〈xy〉−〈x〉〈y〉)/σ x σ y ) between mCherry and mCerulean with a given combination of MREs to the same measure in the absence of MREs. We evaluated the ratio for each eYFP bin (below, around, and above the threshold) for at least three different biological replicates. We then estimated the p values of each ratio with respect to the distributions having as standard deviation the error on the biological replicates and as mean values the Pearson ratio for mCherry and mCerulean with N=0 MRE for the three eYFP intervals.


  1. 1

    Flynt AS, Lai EC. Biological principles of microRNA-mediated regulation: shared themes amid diversity. Nat Rev Genet. 2008; 9(11):831–42. doi:10.1038/nrg2455.

  2. 2

    Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004; 116(2):281–97.

  3. 3

    Bushati N, Cohen SM. MicroRNA functions. Annu Rev Cell Dev Biol. 2007; 23:175–205. doi:10.1146/annurev.cellbio.23.090506.123406.

  4. 4

    Alvarez-Garcia I, Miska EA. MicroRNA functions in animal development and human disease. Development. 2005; 132(21):4653–62. doi:10.1242/dev.02073.

  5. 5

    Esquela-Kerscher A, Slack FJ. Oncomirs – microRNAs with a role in cancer. Nat Rev Cancer. 2006; 6(4):259–69. doi:10.1038/nrc1840.

  6. 6

    Li X, Cassidy JJ, Reinke CA, Fischboeck S, Carthew RW. A microRNA imparts robustness against environmental fluctuation during development. Cell. 2009; 137(2):273–82. doi:10.1016/j.cell.2009.01.058.

  7. 7

    Inui M, Martello G, Piccolo S. MicroRNA control of signal transduction. Nat Rev Mol Cell Biol. 2010; 11(4):252–63. doi:10.1038/nrm2868.

  8. 8

    Osella M, Bosia C, Corá D, Caselle M. The role of incoherent microRNA-mediated feedforward loops in noise buffering. PLoS Comput Biol. 2011; 7(3):1001101.

  9. 9

    Bosia C, Osella M, El Baroudi M, Corá D, Caselle M. Gene autoregulation via intronic microRNAs and its functions. BMC Syst Biol. 2012; 6:131.

  10. 10

    Schmiedel JM, Klemm SL, Zheng Y, Sahay A, Blüthgen N, Marks DS, et al.MicroRNA control of protein expression noise. Science. 6230; 348:128–32.

  11. 11

    Arvey A, Larsson E, Sander C, Leslie CS, Marks DS. Target mRNA abundance dilutes microRNA and siRNA activity. Mol Syst Biol. 2010; 6(1):363. doi:10.1038/msb.2010.24.

  12. 12

    Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi P. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language?Cell. 2011; 146(3):353–8.

  13. 13

    Sumazin P, Yang X, Chiu HS, Chung WJ, Iyer A, Llobet-Navas D, et al.An extensive microRNA-mediated network of RNA–RNA interactions regulates established oncogenic pathways in glioblastoma. Cell. 2011; 147(2):370–81.

  14. 14

    Tay Y, Kats L, Salmena L, Weiss D, Tan SM, Ala U, et al.Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. Cell. 2011; 147(2):344–57.

  15. 15

    Karreth F, Tay Y, Perna D. In vivo identification of tumor-suppressive PTEN ceRNAs in an oncogenic BRAF-induced mouse model of melanoma. Cell. 2011; 147(2):382–95.

  16. 16

    Ala U, Karreth F, Bosia C, Pagnani A, Taulli R, Léopold V, et al.Integrated transcriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments. Proc Natl Acad Sci. 2013; 110(18):7154–9.

  17. 17

    Mukherji S, Ebert M, Zheng G, Tsang J, Sharp P, van Oudenaarden A. MicroRNAs can generate thresholds in target gene expression. Nat Genet. 2011; 43(9):854–9.

  18. 18

    Elf J, Paulsson J, Berg O, Ehrenberg M. Near-critical phenomena in intracellular metabolite pools. Biophys J. 2003; 84:154–70.

  19. 19

    Buchler N, Louis M. Molecular titration and ultrasensitivity in regulatory networks. J Mol Biol. 2008; 384:1106–19.

  20. 20

    Bosia C, Pagnani A, Zecchina R. Modelling competing endogenous RNA networks. PLoS ONE. 2013; 8(6):66609.

  21. 21

    Figliuzzi M, Marinari E, De Martino A. MicroRNAs as a selective channel of communication between competing RNAs: a steady-state theory. Biophys J. 2013; 104(5):1203–13.

  22. 22

    Riba A, Bosia C, El Baroudi M, Ollino L, Caselle M. A combination of transcriptional and microRNA regulation improves the stability of the relative concentrations of target genes. PLoS Comput Biol. 2014; 10(2):1003490.

  23. 23

    Levine E, Zhang Z, Kuhlman T, Hwa T. Quantitative characteristics of gene regulation by small RNA. PLoS Biol. 2007; 5(9):229. doi:10.1371/journal.pbio.0050229.

  24. 24

    Cookson NA, Mather WH, Danino T, Mondragón-Palomino O, Williams RJ, Tsimring LS, et al.Queueing up for enzymatic processing: correlated signaling through coupled degradation. Mol Syst Biol. 2011; 7(1):561. doi:10.1038/msb.2011.94.

  25. 25

    Hussein R, Lim HN. Disruption of small RNA signaling caused by competition for Hfq. Proc Natl Acad Sci. 2011; 108(3):1110–5. doi:10.1073/pnas.1010082108.

  26. 26

    Thomson D, Dinger M. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016; 17(5):272–83.

  27. 27

    Denzler R, Agarwal V, Stefano J, Bartel DP, Stoffel M. Assessing the ceRNA hypothesis with quantitative measurements of miRNA and target abundance. Mol Cell. 2014; 54(5):766–76.

  28. 28

    Bosson AD, Zamudio JR, Sharp PA. Endogenous miRNA and target concentrations determine susceptibility to potential ceRNA competition. Mol Cell. 2014; 56:347–59.

  29. 29

    Xu J, Feng L, Han Z, Li Y, Wu A, Shao T, et al.Extensive ceRNA–ceRNA interaction networks mediated by miRNAs regulate development in multiple rhesus tissues. Nucleic Acid Res. 2016; 44:9438–51.

  30. 30

    Maldotti M, Incarnato D, Neri F, Krepelova A, Rapelli S, Anselmi F, et al.The long intergenic non-coding RNA CCR492 functions as a let-7 competitive endogenous RNA to regulate c-Myc expression. Biochim Biophys Acta. 2016; S1874-9399(16):30131–6.

  31. 31

    Sun C, Li S, Zhang F, Xi Y, Wang L, Bi Y, et al.Long non-coding RNA NEAT1 promotes non-small cell lung cancer progression through regulation of miR-377-3p-E2F3 pathway. Oncotarget. 2016; 7(32):51784–814. doi:10.18632/oncotarget.10108.

  32. 32

    Ma CC, Xiong Z, Zhu GN, Wang C, Zong G, Wang HL, et al.Long non-coding RNA ATB promotes glioma malignancy by negatively regulating miR-200a. J Exp Clin Cancer Res. 2016; 35(1):90.

  33. 33

    Chou J, Wang B, Zheng T, Li X, Zheng L, Hu J, et al.MALAT1 induced migration and invasion of human breast cancer cells by competitively binding miR-1 with cdc42. Biochem Biophys Res Commun. 2016; 472(1):262–9.

  34. 34

    Liu T, Zu CH, Wang SS, Song HL, Wang ZL, Xu XN, et al.Pik3c2a mRNA functions as a miR-124 sponge to facilitate cd151 expression and enhance malignancy of hepatocellular carcinoma cells. Oncotarget. 2016; 7(28):43376–89. doi:10.18632/oncotarget.9716.

  35. 35

    Yuan Y, Liu B, Xie P, Zhang MQ, Li Y, Xie Z, et al.Model-guided quantitative analysis of microRNA-mediated regulation on competing endogenous RNAs using a synthetic gene circuit. Proc Natl Acad Sci. 2015; 112(10):3158–63. doi:10.1073/pnas.1413896112.

  36. 36

    Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al.Global quantification of mammalian gene expression control. Nature. 2011; 473:337–42.

  37. 37

    Marinov GK, Williams BA, McCue K, Schroth GP, Gertz J, Myers RM, et al.From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing. Genome Res. 2014; 24:496–510.

  38. 38

    Mayr C, Bartel DP. Widespread shortening of 3 UTRs by alternative cleavage and polyadenylation activates oncogenes in cancer cells. Cell. 2009; 138:673–84.

  39. 39

    Rüegger S, Großhans H. MicroRNA turnover: when, how, and why. Trends Biochem Sci. 2012; 37(10):436–46.

  40. 40

    Kaern M, Elston TC, Blake WJ, Collins JJ. Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005; 6:451–64.

  41. 41

    Tanase-Nicola S, Warren PB, ten Wolde PR. Signal detection, modularity, and the correlation between extrinsic and intrinsic noise in biochemical networks. Phys Rev Lett. 2006; 97:068102.

  42. 42

    Blake W, Kaern M, Cantor CR, Collins JJ. Noise in eukaryotic gene expression. Nature. 2003; 422:633–7.

  43. 43

    Poliseno L, Salmena L, Zhang J, Carver B, Haveman W, Pandolfi P. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature. 2010; 465:1033–8.

  44. 44

    Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, et al.A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell. 2011; 147:358–69.

  45. 45

    Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al.Real-time quantification of microRNAs by stem–loop RT-PCR. Nucleic Acids Res. 2005; 33(20):179.

  46. 46

    Pantoja-Hernández L, Martínez-García JC. Retroactivity in the context of modularly structured biomolecular systems. Front Bioeng Biotechnol. 2015; 3:85.

  47. 47

    Del Vecchio D, Ninfa AJ, Sontag ED. Modular cell biology: retroactivity and insulation. Mol Syst Biol. 2008; 4:161.

  48. 48

    Bessarabova M, Kirillov E, Shi W, Bugrim A, Nikolsky Y, Nikolskaya T. Bimodal gene expression patterns in breast cancer. BMC Genomics. 2010; 11(Suppl 1):S8.

  49. 49

    Shalek AK, Satija R, Adiconis X, Gertner RS, Gaublomme JT, Raychowdhury R, et al.Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature. 2013; 498:236–40.

  50. 50

    Sætrom P, Heale BSE, Snøve Jr O, Aagaard L, Alluin J, Rossi JJ. Distance constraints between microRNA target sites dictate efficacy and cooperativity. Nucleic Acids Res. 2007; 35(7):2333–42.

  51. 51

    Lai X, Schmitz U, Gupta SK, Bhattacharya A, Kunz M, Wolkenhauer O, et al.Computational analysis of target hub gene repression regulated by multiple and cooperative miRNAs. Nucleic Acids Res. 2012; 40(18):8818–34.

  52. 52

    Schmitz U, Lai X, Winter F, Wolkenhauer O, Vera J, Gupta SK. Cooperative ene regulation by microRNA pairs and their identification using a computational workflow. Nucleic Acids Res. 2014; 42(12):7539–52.

  53. 53

    Broughton JP, Lovci MT, Huang JL, Yeo GW, Pasquinelli AE. Pairing beyond the seed supports microRNA targeting specificity. Mol Cell. 2016; 64:320–33.

  54. 54

    Li G, Burkhardt D, Gross C, Weissman J. Quantifying absolute protein synthesis rates reveals principles underlying allocation of cellular resources. Cell. 2014; 157:624–35.

  55. 55

    Guo X, Deng L, Deng K, Wang H, Shan T, Zhou H, et al.Pseudogene PTENP1 suppresses gastric cancer progression by modulating PTEN. Anticancer Agents Med Chem. 2016; 16(4):456–64.

  56. 56

    Jeyapalan Z, Deng Z, Shatseva T, Fang L, He C, Yang B. Expression of cd44 3 -untranslated region regulates endogenous microRNA functions in tumorigenesis and angiogenesis. Nucleic Acids Res. 2010; 39(8):3026–41.

  57. 57

    Koushik SV, Chen H, Thaler C, Puhl HL, Vogel SS. Cerulean, Venus, and Venus Y67C FRET reference standards. Biophys J. 2006; 91(12):99–101.

  58. 58

    Beier KT, Samson ME, Matsuda T, Cepko CL. Conditional expression of the TVA receptor allows clonal analysis of descendents from Cre-expressing progenitor cells. Dev Biol. 2011; 353(2):309–20.

  59. 59

    Jordan M, Schallhorn A, Wurm F. Transfecting mammalian cells: optimization of critical parameters affecting calcium-phosphate precipitate formation. Nucleic Acids Res. 1996; 24(4):596–601.

  60. 60

    Data and scripts for the quantification of transcripts. doi:10.5281/zenodo.166289. Accessed 9 Jan 2017

  61. 61

    FACS FCS raw data and scripts for the cotransfection experiments. [doi: Accessed 9 Jan 2017].

Download references


The authors would like to thank C Campa, E Fraenkel, I Franco, M Osella, and A Weisse for carefully reading the draft, and U Ala, F Balzac, M Caselle, M Falcone, PP Pandolfi, P Provero, and A Riba for helpful discussions.


LC was supported by a fellowship of Fondazione Umberto Veronesi. C Baldassi and RZ acknowledge the European Research Council for grant 267915. FDC is supported by Telethon Foundation grants GGP12095 and GGP13081, Associazione Italiana per la Ricerca sul Cancro (AIRC) grant IG17527, and the EPIGEN consortium.

Availability of data and materials

The data and scripts for the quantification of transcripts can be found at GitHub [60] with a GNU General Public License as published by the Free Software Foundation, either version 3 or later. The FACS FCS raw data and scripts for the cotransfection experiments can be found on Figshare [61] and they have an MIT license.

Authors’ contributions

C Bosia, AP, and RZ derived the stochastic model. C Bosia, FDC, ET, AP, and RZ designed the experiments. ET, FDC, FC, and RZ contributed reagents, materials, and analysis tools. C Bosia, FS, LC, DB, and ET performed the experiments. C Bosia, C Baldassi, and AP analyzed the data. C Bosia, FS, LC, FDC, AP, and RZ wrote the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Ethical approval was not needed for the study.

Author information

Correspondence to Carla Bosia.

Additional information

Carla Bosia and Francesco Sgrò are co-first authors

Andrea Pagnani and Riccardo Zecchina are co-last authors

Additional file

Additional file 1

Supplementary Material. Detailed description of the modeling and data analysis procedures. (PDF 1106 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Bosia, C., Sgrò, F., Conti, L. et al. RNAs competing for microRNAs mutually influence their fluctuations in a highly non-linear microRNA-dependent manner in single cells. Genome Biol 18, 37 (2017) doi:10.1186/s13059-017-1162-x

Download citation


  • Post-transcriptional cross-regulation
  • MicroRNA target synchronization
  • Bimodality
  • Single cell
  • Stochastic modelling