ceRNA crosstalk stabilizes protein expression and affects the correlation pattern of interacting proteins

Gene expression is a noisy process and several mechanisms, both transcriptional and posttranscriptional, can stabilize protein levels in cells. Much work has focused on the role of miRNAs, showing in particular that miRNA-mediated regulation can buffer expression noise for lowly expressed genes. Here, using in silico simulations and mathematical modeling, we demonstrate that miRNAs can exert a much broader influence on protein levels by orchestrating competition-induced crosstalk between mRNAs. Most notably, we find that miRNA-mediated cross-talk (i) can stabilize protein levels across the full range of gene expression rates, and (ii) modifies the correlation pattern of co-regulated interacting proteins, changing the sign of correlations from negative to positive. The latter feature may constitute a potentially robust signature of the existence of RNA crosstalk induced by endogenous competition for miRNAs in standard cellular conditions.

(a) ceRNA cross-talk can stabilize the level of highly expressed proteins (with respect to the case in which no competition takes place); (b) the ceRNA effect alters the correlation pattern of co-regulated interacting proteins, particularly by turning its sign from negative to positive; (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 in 27,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 m T (for 'target') and m C (for 'competitor'). Both serve as substrates for protein synthesis and the respective A single miRNA species negatively controls the expression of the ceRNA species ceRNA T ('target' , level m T ) and ceRNA C ('competitor' , level m C ), to which it can bind with rates + k T and + k C , respectively. The functional products of the ceRNAs, proteins p T and p C , can eventually interact to form a complex C P . The key control parameter is the target's transcription rate b T . See Fig. 7 for a detailed scheme that includes all processes.
Scientific RepoRts | 7:43673 | DOI: 10.1038/srep43673 products (p T and p C ) can interact to form a complex (C P ). 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 algorithm 42 , whereas analytical estimates for quantities like correlation functions can be obtained via the Linear Noise Approximation 43 (LNA); see Methods for details.
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 (b T ) while keeping all other parameters fixed, the expression of the target's functional product p T undergoes a crossover from a repressed regime with low copy numbers to a free (unrepressed) regime in which its level increases roughly linearly with b T . The crossover gets sharper as the miRNA-target interaction strength + k T increases. The ability of miRNAs to generate threshold-linear expression profiles via molecular titration was pointed out already in several studies 27,28,44,45 . In the crossover region, p T displays strong sensitivity to small changes in b T , as testified by peaks in the coefficient of variation (CV) that get more marked and shift to lower values of protein concentration as + k T increases, see Fig. 2C. In this regime, called 'susceptible' in Figliuzzi et al. 28 , miRNA and mRNA levels are nearly equimolar 27,28,34 . Contrasting this behaviour with the standard Poissonian CV obtained for an unregulated protein ( = + k 0 T , 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 genes 24 . When a competitor is present (Fig. 2D-F), one observes that, upon modulating b T , the level of p C starts changing when b T is around the crossover region, reflecting the effective positive coupling between ceRNAs m T and m C 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. We shall now analyze in more detail the influence of miRNA sponging and ceRNA competition on protein expression noise and on the PPI. Following 36,48,49 , the effectiveness of the regulatory channel linking an input variable, b T in this case, to an output variable O (e.g., the target protein level p T or the level of the protein complex C P ) will be characterized by its capacity I max , defined as the maximum mutual information between b T and O that can be achieved by changing the input distribution p(b T ) while keeping the conditional distribution p(O|b T ) (the 'channel' , which stochastically returns a value of O upon presenting input b T ) fixed: . Note that no PPI is considered in this case.
is the output distribution. We will work under the assumption that p(O|b T ) is Gaussian and its variance is small for each b T ('Small Noise Approximation'), in which case the above problem has been shown to have a simple analytical solution 49,50 (see the protocol for computing capacities in Methods). Note that the input variable is constrained to vary between fixed bounds b T min and b T max and that, correspondingly, the output varies between O min and O max . ceRNA cross-talk enhances the stability of highly expressed proteins. Recent experimental work 24 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.
In Fig. 3B we show how the capacity of the protein expression channel changes with the strength + k C of the interaction between the competitor, ceRNA C , and the miRNA. For small + k C , I max (p T , b T ) expectedly tends to the value obtained for a miRNA-regulated protein as the effect of the competitor gets weaker and weaker. For large + k C , 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 + k C , instead, I max peaks, signalling that the expression of p T 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 ceRNA C , Δ C , defined as the difference between the largest and smallest steady state values of m C that are obtained by changing the input variable b T between its smallest and largest allowed values b T max and b T min : Based on the ceRNA effect features shown in Fig. 2E, a large derepression size indicates that the ceRNA effect is active. Clearly, see 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 induced [53][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 expression 36 . 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%. 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 levels [27][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 ρ(p T , p C ) is largest close to the regime where the ceRNA effect is strongest and ρ(m T , m C ) peaks, while it becomes negative when the ceRNA-ceRNA correlation weakens.
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 human 58 . Experiments in laryngeal squamous cells have revealed that, under CCND1 over-expression, the expression level of CDK4 increases 59 . Likewise, ITGA6 down-regulation results in ITGB4 decrease in cells where the ITGA6:ITGB4 complex is present 60 . Both ITGA6 and ITGB4 are known to share common miRNA regulators in humans 38 . 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 (p T ) and competitor (p C ) proteins occurs only in a narrow window of parameters around the quasi-equimolar 27 or susceptible 28 regime. It is somewhat striking that the experimental studies of regulated binary complexes discussed here 58-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 C p . 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 + k T , where I max is slightly above the value obtained for low values of + k T and + k C (corresponding to the case of unregulated transcript). Interestingly, though, the channel capacity seems to become generically larger as + k C grows, suggesting that, under stronger competition, a more efficient tuning of the complex level may be achieved in a broader range of values of + k T .

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 expression 61 that exert regulatory roles in situations as different as neuronal regulation 62,63 , brain morphogenesis 64 , muscle cell differentiation 65 , stem cell division 66 , glucose and lipid metabolism 67 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 expression 68 . 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 spectrum 69 ) 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 miRNAs 14 , 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 understood 27,28,36,47,70 , its significance in vivo is debated [16][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 motifs 25,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 reporters 24 .
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 regulation 28,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 networks 38 ), 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 proteins 71 . 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 complex 38 . Likewise, miR-200 has been identified as a powerful marker of the epithelial phenotype of cancer cells 72 , while miR141-200c targets β-catenin, the downstream effector of the Wnt proliferation pathway 73 , and miR-545 targets complex-forming CCND1 and CDK4 with potentially suppressive effects on the proliferation of lung cancer cells 58 . 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.  Species m T , m C , μ, p T , p C are synthesized (degraded) with the rates b T , b C , β, α T , α C (d T , d C , δ, δ T , δ C ) correspondingly. Protein complex C p undergoes spontaneous degradation with a rate δ p . Finally, c T and c C 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.

Methods
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 = {m T , m C , μ, c T , c C , p T , p T } and by x its steady state value. Assuming that the divergence from the steady state, δ = − x x x, is small and expanding Eq. (3) around x, 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 that 43 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 reactions 42 . 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 rates 42 . 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 (b T ) (variance of the output O) obtained at stationarity by the GA for any given value of the input variable b T in the range b b ( , )

T T min max
Step 2: Compute the optimal input distribution ( ) T is the mean expression level of the output for the input b T . Following 49,50 , in the limit of small noise the channel's capacity coincides with the mutual information obtained when the input variable b T is distributed according to P opt (b T ).
Step 3: Generate an input signal according to P opt (b T ), record corresponding output signal O(b T ) and calculate the joint probability distribution P(b T , O). Step