Introduction

The control of gene expression noise is a challenge faced, at some degree, by cells in all organisms, as each post-transcriptional step of the genetic regulatory chain, including translation alone1, can potentially amplify transcription noise. Biological functionality however often requires finely-tuned protein levels. Cells therefore employ a variety of strategies to ensure that protein noise is buffered2,3,4,5,6,7,8,9. Regulatory RNAs, and microRNAs (miRNAs) in particular, are thought to play a major role in this respect10,11.

miRNAs comprise a large number of short, endogenously expressed non-coding RNA species that are significantly conserved among invertebrates and vertebrates and whose expression is strongly tissue-specific12. They act primarily as negative controllers of gene expression, by silencing translation and/or catalyzing mRNA destabilization after sequence-specific binding to their targets. They can however also bind non-coding RNA species like pseudogenes and long non-coding RNAs (lncRNAs)13,14. In some cases, ‘sponging’ of miRNAs by lncRNAs has been found to contribute significantly to the adjustment of miRNA levels in the cell15. Overall, miRNAs appear today as key regulators of a very broad class of RNA molecules.

Recent work, both experimental and in silico, has started to uncover the complex and highly heterogeneous network of miRNA-RNA interactions, showing that hundreds of genes may be directly repressed by individual miRNAs, albeit modestly in some cases16,17,18,19,20,21. The search for a possible functional rationale for the breadth of miRNA-RNA couplings has increasingly focused on containment of gene expression noise20,21,22,23,24,25. Several studies have investigated the effects of miRNA-mediated regulation on RNA and protein levels24,25,26. Most recently, integrated theoretical and experimental work24 has shown that miRNA control stabilizes the levels of lowly expressed genes while having the opposite effect on highly expressed ones. This scenario, which mirrors that obtained for miRNA-mediated regulation of transcript levels27,28, suggests that other mechanisms may confer robustness at high expression levels or, more generally, across the entire expression range. Here we argue that the recently hypothesized ‘ceRNA effect’29 might constitute such a mechanism.

The ceRNA effect (whereby ceRNA stands for ‘competing endogenous RNAs’) consists, in essence, of a positive effective interaction that can arise between transcripts that are targeted by the same miRNA species due to competition29,30,31,32. The ability of the ceRNA mechanism to mediate RNA cross-talk has been investigated in detail in silico27,28,33, while experimental validations have been obtained by over-expressing miRNAs or targets16 and in specific conditions related to disease or differentiation15,34,35. Its relevance in standard physiological conditions is less clear16,17. On one hand, effective ceRNA cross-talk appears to depend strongly on the kinetics of miRNA-target interactions, on gene silencing mechanisms and on the relative abundance of regulators and targets, and therefore may require considerable fine tuning27,28,36,37. On the other, miRNAs can achieve optimal or nearly optimal regulation of competing transcripts in a broad range of parameter values by exploiting heterogeneities in e.g. binding kinetics or target degradation pathways36. A natural question is whether proteins, i.e. the functional products of coding transcripts, may also benefit from miRNA-mediated control when the respective mRNAs are competing to bind miRNAs. Interestingly, studies of protein-protein interaction (PPI) networks have revealed that interacting proteins tend to be regulated by miRNA clusters, i.e. by co-expressed miRNAs37,38,39,40,41. More specifically, it was found that (a) individual or co-expressed miRNAs frequently target the mRNAs of several components of protein complexes38, (b) direct miRNA targets and their partners show significant modularity at the level of the corresponding proteins in the human PPI network39, and (c) the products of transcripts targeted by the same miRNAs are more connected in the human PPI network than expected by chance40. These observations suggest that miRNAs may be implicated in the fine control of interacting proteins, making the ceRNA effect a rather natural mechanism through which such control could be exerted.

Here we quantify the impact of ceRNA competition on protein fluctuations and on protein complex formation by mathematical modeling and in silico experiments. Within a minimal stochastic description that includes both transcriptional and post-transcriptional control, we show that:

  1. a

    ceRNA cross-talk can stabilize the level of highly expressed proteins (with respect to the case in which no competition takes place);

  2. b

    the ceRNA effect alters the correlation pattern of co-regulated interacting proteins, particularly by turning its sign from negative to positive;

  3. c

    miRNA recycling enhances the suppression of protein expression noise through the ceRNA effect across the entire range of expression levels.

These results have significant implications. First, they suggest that ceRNA cross-talk may be crucial for the fine tuning of protein levels, thereby pointing to a further explanation for the abundance of lncRNAs and pseudogenes (i.e. of miRNA ‘sponges’) in the human transcriptome. Secondly, they indicate that a positive correlation between co-regulated subunits of a protein complex may provide, for a limited but significant set of cases, the simple and direct proof of active ceRNA cross-talk in standard physiological conditions that has been so far lacking.

Results

Model description and basic properties

As a basis, we consider the model of a ceRNA network studied previously in27,28,36, with the addition of protein synthesis and a protein complex formation step (see Fig. 1). In short, a miRNA species negatively regulates the expression of two ceRNAs, whose levels are denoted respectively as mT (for ‘target’) and mC (for ‘competitor’). Both serve as substrates for protein synthesis and the respective products (pT and pC) can interact to form a complex (CP). By turning specific parameters on or off the above basic model allows to tackle different situations. A detailed description of the model, including the various parameters, is given in the Methods. Its dynamics, formalized in terms of stochastic differential equations, can be simulated using the Gillespie algorithm42, whereas analytical estimates for quantities like correlation functions can be obtained via the Linear Noise Approximation43 (LNA); see Methods for details.

Figure 1: Model schematics and main parameters.
figure 1

A single miRNA species negatively controls the expression of the ceRNA species ceRNAT (‘target’, level mT) and ceRNAC (‘competitor’, level mC), to which it can bind with rates and , respectively. The functional products of the ceRNAs, proteins pT and pC, can eventually interact to form a complex CP. The key control parameter is the target’s transcription rate bT. See Fig. 7 for a detailed scheme that includes all processes.

Full size image

The basic features of miRNA-mediated regulation that emerge from this model at stationarity with and without ceRNA competition (and in absence of PPI) are summarized in Fig. 2. In absence of ceRNA competition (Fig. 2A–C), upon varying the synthesis rate of the target (bT) while keeping all other parameters fixed, the expression of the target’s functional product pT undergoes a crossover from a repressed regime with low copy numbers to a free (unrepressed) regime in which its level increases roughly linearly with bT. The crossover gets sharper as the miRNA-target interaction strength increases. The ability of miRNAs to generate threshold-linear expression profiles via molecular titration was pointed out already in several studies27,28,44,45. In the crossover region, pT displays strong sensitivity to small changes in bT, as testified by peaks in the coefficient of variation (CV) that get more marked and shift to lower values of protein concentration as increases, see Fig. 2C. In this regime, called ‘susceptible’ in Figliuzzi et al.28, miRNA and mRNA levels are nearly equimolar27,28,34. Contrasting this behaviour with the standard Poissonian CV obtained for an unregulated protein (, black line in Fig. 2B and C) one sees that miRNA control generically buffers expression noise for lowly expressed genes while it amplifies noise for highly expressed genes, in agreement with recent experimental work quantifying protein expression noise in miRNA-regulated genes24. When a competitor is present (Fig. 2D–F), one observes that, upon modulating bT, the level of pC starts changing when bT is around the crossover region, reflecting the effective positive coupling between ceRNAs mT and mC known as the ‘ceRNA effect’. The emergence and major features of the ceRNA effect at the level of transcripts have been characterized in refs 27,28,33,36,46,47.

Figure 2: Dependence of mean protein expression levels and relative fluctuations (CV) on the transcription rate bT of the target.
figure 2

Panels (B) and (C) describe the case of a simple miRNA-regulated target, shown in panel (A). In (B), increases in the direction of the arrow (specifically, for orange, purple, blue, black curves respectively). Panels (E) and (F) describe the case of a target regulated through ceRNA competition, depicted in (D), for and . Note that no PPI is considered in this case.

Full size image

We shall now analyze in more detail the influence of miRNA sponging and ceRNA competition on protein expression noise and on the PPI. Following36,48,49, the effectiveness of the regulatory channel linking an input variable, bT in this case, to an output variable O (e.g., the target protein level pT or the level of the protein complex CP) will be characterized by its capacity Imax, defined as the maximum mutual information between bT and O that can be achieved by changing the input distribution p(bT) while keeping the conditional distribution p(O|bT) (the ‘channel’, which stochastically returns a value of O upon presenting input bT) fixed:

where is the output distribution. We will work under the assumption that p(O|bT) is Gaussian and its variance is small for each bT (‘Small Noise Approximation’), in which case the above problem has been shown to have a simple analytical solution49,50 (see the protocol for computing capacities in Methods). Note that the input variable is constrained to vary between fixed bounds and and that, correspondingly, the output varies between Omin and Omax.

ceRNA cross-talk enhances the stability of highly expressed proteins

Recent experimental work24 has revealed that, while miRNA regulation suppresses protein noise for lowly expressed genes, at high expression levels noise is generically larger. This effect can be surprisingly reversed through ceRNA competition. In particular, by titrating miRNA molecules away from their target, a competitor can enhance the target’s stability. This is shown in Fig. 3A, where the CVs of proteins that are post-transcriptionally unregulated, regulated by a miRNA and regulated by the ceRNA effect are compared. One sees that the relative fluctuations of the target at high expression levels can be substantially reduced by the onset of ceRNA competition with respect to the miRNA-regulated case without significantly affecting the low expression regime. In addition, ceRNA regulation generates a fluctuation scenario similar to the Poissonian one that characterizes an unregulated protein, showing that competition buffers noise essentially by de-repressing the target.

Figure 3: ceRNA competition can stabilize highly expressed proteins.
figure 3

(A) Coefficient of variation of pT as a function of the mean protein level for a post-transcriptionally unregulated protein (black line, and ), a miRNA-regulated protein (red line, and ) and a ceRNA-regulated protein (blue line, and ). (B) Capacity of the target’s expression channel as a function of the miRNA-competitor interaction strength. Color code same as in panel A. (C) Derepression size ΔC of the competitor as a function of the miRNA-competitor interaction strength in the case of ceRNA regulation (same parameters as panel B).

Full size image

In Fig. 3B we show how the capacity of the protein expression channel changes with the strength of the interaction between the competitor, ceRNAC, and the miRNA. For small , Imax(pT, bT) expectedly tends to the value obtained for a miRNA-regulated protein as the effect of the competitor gets weaker and weaker. For large , on the other hand, the competitor tends to sponge all miRNAs away from the target, therefore leading to a capacity close to that of a simple transcriptional control unit. For intermediate values of , instead, Imax peaks, signalling that the expression of pT can be tuned more efficiently than by transcriptional or miRNA-mediated regulation alone. That the ceRNA effect is at origin of this behaviour can be checked by measuring the derepression size of ceRNAC, ΔC, defined as the difference between the largest and smallest steady state values of mC that are obtained by changing the input variable bT between its smallest and largest allowed values and :

Based on the ceRNA effect features shown in Fig. 2E, a large derepression size indicates that the ceRNA effect is active. Clearly, see Fig. 3C, ΔC is markedly larger for intermediate values of , thereby pointing to the onset of ceRNA crosstalk.

miRNA recycling increases the potential of miRNA-mediated gene regulatory circuits to suppress gene expression noise

miRNA-ceRNA complexes can be processed in different ways, including via unbinding and complex degradation (‘stoichiometric decay’)12,51,52. However, if miRNAs have a perfect complementarity with the target, catalytic cleavage of the latter can be induced53,54,55,56,57, after which miRNAs are likely to be available again for interaction with a new target (‘miRNA recycling’). This catalytic channel of complex processing may effectively increase the population of free miRNA regulators with respect to targets, thereby enabling a more subtle control of gene expression36. Figure 4 elucidates its role in controlling protein levels. By increasing the efficiency of mRNA cleavage at the target node (i.e. for faster rates of catalytic processing κT, see Methods), protein expression noise improves significantly with respect to the case of slow catalytic processing. Notice that the same effect is obtained in a simpler miRNA-regulated element. In specific, while for a ceRNA regulated target the stability of the protein levels improves by about 20%, for a miRNA-regulated target the improvement may be more relevant – up to about 40%.

Figure 4: Noise reduction in miRNA-mediated circuits due to the effective miRNA-recycling at the target node.
figure 4

Fast and slow (high and low κT, see Methods) miRNA recycling scenarios at the target node are shown, respectively, by dashed and solid lines for a miRNA-regulated protein (red, and ) and a ceRNA-regulated protein (blue, and ). The black curve describes the case of an unregulated target ( and ).

Full size image

The ceRNA effect alters the correlation pattern of co-regulated interacting proteins

Competition to bind regulatory miRNAs can lead to a positive correlation among co-regulated transcripts, who can respond both to fluctuations in miRNA levels and, at fixed miRNA level, to changes in ceRNA levels27,28,29,36. This is shown in Fig. 5A, where one also sees that the corresponding (non-interacting) proteins have a similar correlation pattern, as measured by the Pearson coefficient between the levels of free molecules. Proteins that form a complex are instead negatively correlated by the PPI in absence of co-regulation at the level of transcripts (in which case transcripts are uncorrelated), as shown in Fig. 5B. The negative correlation is caused by the protein complex binding kinetics alone. Interestingly, the ceRNA effect can reverse this scenario, i.e. in presence of upstream miRNA competition interacting proteins can become positively correlated (Fig. 5C). In particular, the Pearson coefficient ρ(pT, pC) is largest close to the regime where the ceRNA effect is strongest and ρ(mT, mC) peaks, while it becomes negative when the ceRNA-ceRNA correlation weakens.

Figure 5: Correlation patterns for interacting proteins measured by the Pearson coefficient ρ as a function of the mean target level.
figure 5

(A) Case of non-interacting proteins translated from competing transcripts ( and ). (B) Case of interacting proteins translated from post-transcriptionally unregulated transcripts (). (C) Case of interacting proteins translated from competing transcripts ( and ). Lines (blue and orange) correspond to analytical results obtained by the Linear Noise Approximation (see Methods), markers (black and purple) to results from stochastic simulations by the Gillespie algorithm.

Full size image

Several proteins forming binary complexes and known to share miRNA regulators display a positive correlation. For instance, both subunits of the CCND1:CDK4 complex are regulated by miR-545 in human58. Experiments in laryngeal squamous cells have revealed that, under CCND1 over-expression, the expression level of CDK4 increases59. Likewise, ITGA6 down-regulation results in ITGB4 decrease in cells where the ITGA6:ITGB4 complex is present60. Both ITGA6 and ITGB4 are known to share common miRNA regulators in humans38. Based on the above results, a positive correlation between subunits of a protein complex may therefore be explained in terms of miRNA-mediated cross-talk at the level of the respective transcripts.

Incidentally, we note that, as shown in Fig. 5C, the onset of a positive correlation between target (pT) and competitor (pC) proteins occurs only in a narrow window of parameters around the quasi-equimolar27 or susceptible28 regime. It is somewhat striking that the experimental studies of regulated binary complexes discussed here58,59,60 report a positive correlation between subunits. This suggests that, at least for these systems, kinetic parameters could be globally tuned by evolution so as to fit this narrow window. A potential added advantage of such a scenario lies, as we shall now see, in the buffering of protein complex noise that can accompany it.

Co-regulation by a miRNA can effectively control the level of a binary protein complex

If protein fluctuations get correlated due to the ceRNA effect, it might be possible to exploit the same mechanism to fine tune the level of the protein complex Cp. Figure 6 shows that the capacity of the protein complex synthesis channel is weakly modulated by the miRNA-ceRNA binding rate that quantifies the strength of miRNA regulation on the competitor (see Fig. 3B). On the other hand, optimal regulation is achieved at intermediate values of , where Imax is slightly above the value obtained for low values of and (corresponding to the case of unregulated transcript). Interestingly, though, the channel capacity seems to become generically larger as grows, suggesting that, under stronger competition, a more efficient tuning of the complex level may be achieved in a broader range of values of .

Figure 6: Capacity of the protein complex synthesis channel, Imax(Cp, bT), as a function of the miRNA-ceRNA binding strengths of the target () and the competitor ().
figure 6

To allow for comparisons, simulations were performed at fixed output variation range (with ) and variable bT.

Full size image

Discussion

A wide variety of biological functions has been assigned over time to non-coding RNA molecules. Most interestingly, perhaps, they can regulate gene expression at the post-transcriptional stage. Eukaryotic miRNAs are post-transcriptional micromanagers of gene expression61 that exert regulatory roles in situations as different as neuronal regulation62,63, brain morphogenesis64, muscle cell differentiation65, stem cell division66, glucose and lipid metabolism67 as well as in many disease states. The identification of viral miRNAs furthermore suggests that viruses may use them to interfere with the host’s gene expression68. Therapeutics targeting miRNA levels therefore appear to be particularly promising tools. Their viability however has to be evaluated against the broad microenvironment in which miRNAs operate. It is now clear that miRNAs and their targets are linked in a complex transcriptome-scale interaction network and that the effective strength of miRNA-induced repression depends tightly on miRNA-RNA binding and unbinding free energies (that are predicted to show a wide spectrum69) and on molecular levels. In such a context, understanding how perturbations probing one node may propagate to other nodes, thereby affecting gene expression as a whole, is far from trivial. Moreover, pseudogenes and long non-coding RNAs can compete with coding transcripts to bind miRNAs14, and the effects induced by competition can be especially hard to quantify. Indeed, the establishment of an effective positive interaction between different targets of the same miRNAs (whose effectiveness is supported by several perturbation-based experimental studies) may potentially contribute significantly to both the overall gene expression profile and, dynamically, to the re-shaping of the proteome in response to perturbations. While many of its features are well understood27,28,36,47,70, its significance in vivo is debated16,17,18. On the other hand, miRNAs have been shown to be able to confer precision to expression levels either by themselves or in combination with spefic motifs25,26. Recent experiments have also shown that, in cells with low expression of the reporter carrying a binding site for the miRNA protein, protein noise was reduced compared to a control, while fluctuations were increased for highly expressed reporters24.

Our study establishes a link between different miRNA-mediated functions by showing that buffering of protein expression noise can be enhanced by competition. In specific, competition by itself reduces noise on highly expressed genes, while catalytic processing of the miRNA-ceRNA complex provides a further degree of freedom through which noise can be controlled. As was found to be the case for other emergent properties of miRNA-based regulation28,36, these effects are enhanced by kinetic heterogeneities, providing further insight into the possible functions served by the evolutionarily-selected, broad distribution of rates found in miRNA-ceRNA networks. In addition, for interacting proteins whose mRNAs are regulated by the same miRNA (as frequently found in actual regulatory networks38), our study unveils several new features. First, we argued that, due to the strong ceRNA effect, a positive correlation may be established between the two sub-units of a complex, reversing the negative sign that is induced by complex formation in absence of co-regulation by a common modulator. This is a direct consequence of ceRNA cross-talk and may provide a first, important signature of ceRNA cross-talk in vivo. Finally, we have shown that protein complexes may be synthesized with higher precision if its components undergo miRNA-mediated post-transcriptional regulation, although the effect can be modest.

The influence of miRNAs on protein complexes has been experimentally investigated in the context of diseases. For instance, analyses of miRNA-mediated dysregulation of functionally related proteins during prostate cancer progression had identified miRNA-1 and miRNA-16 as master regulators of prostate cancer, since they regulate hubs of the underlying PPI network – the SMAD4 and HDAC proteins71. Other studies have revealed the regulatory role of miR141–200c in the epithelial-to-mesenchymal transition due to the orchestrated regulation of the CtBp/Zeb complex38. Likewise, miR-200 has been identified as a powerful marker of the epithelial phenotype of cancer cells72, while miR141–200c targets β-catenin, the downstream effector of the Wnt proliferation pathway73, and miR-545 targets complex-forming CCND1 and CDK4 with potentially suppressive effects on the proliferation of lung cancer cells58. The present work established an in silico framework through which such cases could be quantitatively analyzed by probing the roles of the various parameters that regulate miRNA regulation, competition, cross-talk and PPI. Experimental validation in perturbation-based experiments would provide further understanding on the global effectiveness of miRNAs as controllers of cellular protein profiling.

Methods

Stochastic model

The dynamics of the model, presented in detail in Fig. 7, can be described by the stochastic equations

Figure 7: Schematic representation of the model.
figure 7

Two proteins (pT and pC) translated from 2 distinct ceRNAs (mT and mT) that are regulated by the same miRNA (μ). Proteins pT and pC associate and dissociate to a protein complex Cp with a rate k+ and k− respectively, miRNA binds (unbinds) to the ceRNAs mT and mC with the rates and respectively forming ceRNA:miRNA complexes cT and cC. Species mT, mC, μ, pT, pC are synthesized (degraded) with the rates bT, bC, β, αT, αC (dT, dC, δ, δT, δC) correspondingly. Protein complex Cp undergoes spontaneous degradation with a rate δp. Finally, cT and cC decay catalytically (ceRNA cleavage and miRNA recycling) with the rates κT and κC respectively. Figure-specific values of the kinetic parameters are reported in Table 1.

Full size image
Table 1 Parameters values used to obtain the figures.
Full size table

where i ∈ {C, T}. The different ξ-terms stand for the contributions to the noise coming from different processes. Each of these terms is assumed to have zero mean and variance given by

where the over-line stands for an average over time in the steady state. From the stationarity conditions, one finds in particular

By selectively setting some of the parameters appearing above to zero one easily obtains the equations describing the different circuits studied here.

Linear Noise Approximation

Let us denote the vector of molecular levels for the different species involved by x = {mT, mC, μ, cT, cC, pT, pT} and by its steady state value. Assuming that the divergence from the steady state, , is small and expanding Eq. (3) around , at the leading order one gets

where A is the matrix of the first order derivatives evaluated at the steady state, while ξ represents a white noise with zero mean and covariance matrix 〈ξa(t)ξb(t′)〉 = Γabδ(t − t′), where the indices a and b range over the components of x (non-zero elements of the covariance matrix are given by Eq. (4)). From Eq. (6) it follows that43

where B′s and λ′s denote the eigenvectors and eigenvalues of A. The above formula can be used to estimate correlations and Pearson coefficients of all molecular species involved in the system.

Stochastic simulations

We have employed the Gillespie algorithm (GA), a classical stochastic simulation method that computes the population dynamics of N well-mixed molecular species interacting through one of M reactions42. The dynamics of a biochemical system is obtained based on the probability P(R, τ) of an event of reaction R to take place in the next time interval of size τ. The latter can be calculated for every set of molecular populations given the chemical reaction rates42. In brief, the GA works as follows:

Step 1 (initialization): Set up initial populations for all molecular species,

Step 2: Draw a pair (R, τ) from P(R, τ),

Step 3: Update molecular populations according to the selected reaction R and advance time by τ,

Step 4: Repeat Steps 2 and 3 until a pre-determined termination time is reached.

See Gibson et al.74 for a more detailed presentation of the method.

A C++ implementation of the algorithm for the network shown in Fig. 7 is available at https://github.com/araksm/Protein-Expression-Noise.

Protocol for computing capacities

Capacities of the different regulatory channels we consider were computed as follows

Step 1: Calculate the output noise σO(bT) (variance of the output O) obtained at stationarity by the GA for any given value of the input variable bT in the range

Step 2: Compute the optimal input distribution

with , where is the mean expression level of the output for the input bT. Following49,50, in the limit of small noise the channel’s capacity coincides with the mutual information obtained when the input variable bT is distributed according to Popt(bT).

Step 3: Generate an input signal according to Popt(bT), record corresponding output signal O(bT) and calculate the joint probability distribution P(bT, O).

Step 4: Calculate the capacity Imax(bT, O) from the data obtained in Step 3 via

where , while Omin and Omax denote maximum and minimum output expressed levels correspondingly.

Additional Information

How to cite this article: Martirosyan, A. et al. ceRNA crosstalk stabilizes protein expression and affects the correlation pattern of interacting proteins. Sci. Rep. 7, 43673; doi: 10.1038/srep43673 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.