Probing the Limits to MicroRNA-Mediated Control of Gene Expression.

According to the 'ceRNA hypothesis', microRNAs (miRNAs) may act as mediators of an effective positive interaction between long coding or non-coding RNA molecules, carrying significant potential implications for a variety of biological processes. Here, inspired by recent work providing a quantitative description of small regulatory elements as information-conveying channels, we characterize the effectiveness of miRNA-mediated regulation in terms of the optimal information flow achievable between modulator (transcription factors) and target nodes (long RNAs). Our findings show that, while a sufficiently large degree of target derepression is needed to activate miRNA-mediated transmission, (a) in case of differential mechanisms of complex processing and/or transcriptional capabilities, regulation by a post-transcriptional miRNA-channel can outperform that achieved through direct transcriptional control; moreover, (b) in the presence of large populations of weakly interacting miRNA molecules the extra noise coming from titration disappears, allowing the miRNA-channel to process information as effectively as the direct channel. These observations establish the limits of miRNA-mediated post-transcriptional cross-talk and suggest that, besides providing a degree of noise buffering, this type of control may be effectively employed in cells both as a failsafe mechanism and as a preferential fine tuner of gene expression, pointing to the specific situations in which each of these functionalities is maximized.

Introduction

The problem of tuning protein expression levels is central for eukaryotic cell functionality. A variety of molecular mechanisms are implemented to guarantee, on one hand, that protein copy numbers stay within a range that is optimal in the given conditions and, on the other, that shifts in expression levels can be achieved efficiently whenever necessary [1-3] (whereby ‘efficiency’ here encompasses both a dynamical characterization, in terms of the times required to shift, and a static one, in terms of moving as precisely as possible from one functional range to another). Quantifying and comparing their effectiveness in different conditions is an important step to both deepen our fundamental understanding of regulatory circuits and to get case-by-case functional insight about why a specific biochemical network has been selected over the others.

As the major direct regulators of gene expression, transcription factors (TFs) are most immediately identified as the key potential modulators of protein levels [4]. In a somewhat simplified picture, one may imagine that a change in amount of a TF can induce a change in the expression level of the corresponding gene, and that the ability to regulate the latter (the output node) via the former (the input node) can be assessed by how strongly the two levels correlate. The effectiveness of a regulatory element is however limited by the stochasticity of intracellular processes, from the TF-DNA binding dynamics to translation [5]. A convenient framework to analyze how noise constrains regulation is provided by information theory [6, 7]. In particular, the simplest situation in which a single TF modulates the expression of a single protein can be characterized analytically under the assumption that the noise affecting the input-output channel is sufficiently small. The mutual information between modulator and target –a convenient quantity through which regulatory effectiveness can be characterized– depends on the distribution of modulator levels and can be maximized over it. Remarkably, in at least one case this maximum has been found to be almost saturated by the actual information flow measured in a living system (for more details see [8, 9]). In other terms, for sufficiently small noise levels in the channel that links TFs to their functional products, one may quantify the optimal regulatory performance achievable in terms of the maximum number of bits of mutual information that can be exchanged between modulator and target.

Several control mechanisms however act at the post-transcriptional level [10-12]. Among these, regulation by small regulatory RNAs like eukaryotic microRNAs (miRNAs) has attracted considerable attention over the past few years [13-15]. In short, miRNAs are small non-coding RNA molecules encoded by nuclear DNA, that can inhibit translation or catalyze degradation of mRNAs when bound to them via protein-mediated base-pairing. miRNAs appear to be crucial in an increasing number of situations ranging from development to disease [16-18]. Their function however can differ significantly from case to case. For instance, they have been well characterized as noise buffering agents in protein expression [19, 20] or as key signaling molecules in stress response [21]. Recently, though, investigations of the miRNA-mediated post-transcriptional regulatory (PTR) network have hinted at a possibly more subtle and complex role. It is indeed now clear [22-27] that the miRNA-RNA network describing the potential couplings stretches across a major fraction of the transcriptome, with a large heterogeneity both in the number of miRNA targets and in the number of miRNA regulators for a given mRNA. The competition effects that may emerge in such conditions suggest that miRNAs may act as channels through which perturbations in the levels of one RNA could be transmitted to other RNA species sharing the same miRNA regulator(s). Such a scenario has been termed the ‘ceRNA effect’, whereby ceRNA stands for ‘competing endogenous RNA’ [28]. In view of its considerable regulatory and therapeutic implications, the ceRNA effect has been extensively analyzed both theoretically and experimentally [29-44].

The apparent ubiquity of potentially cross-talking ceRNAs however raises a number of fundamental questions about the effectiveness of “regulation via competition” per se. Although hundreds of targets are predicted for a single miRNA, observations show that only few of them are sensitive to changes in miRNA expression levels. Most targets are likely to provide a global buffering mechanism through which miRNA levels are overall stabilized [28, 29]. Effective competition between miRNA targets requires that the ratio of miRNA molecules to the number of target sites lies in a specific range, so that the relative abundance of miRNA and RNA species must be tightly regulated for the ceRNA mechanism to operate [29-34]. On the other hand, the magnitude of the ceRNA effect is tunable by the miRNA binding and mRNA loss rates [33, 34, 42]. The performance of a regulatory element, however, does not only depend on kinetic parameters, but also on the range of variability (and possibly on the distribution itself) of modulator levels (e.g. TFs) [8, 9]. The maximal regulatory effectiveness of a given genetic circuit –quantifying how precisely the output level can be determined by the input level– can therefore generically be obtained by solving an optimization problem over the distribution of inputs. This type of approach provides an upper bound to the effectiveness of a regulatory mechanism as well as indications concerning which parameters, noise sources and/or interactions most hamper its performance.

It would be especially important to understand in which conditions the degree of control of the output variable (i.e. the ceRNA/protein level) that can be accomplished through post-transcriptional miRNA-mediated cross-talk may exceed that obtainable by different regulatory mechanisms. In this work we characterize the maximal regulatory power achievable by miRNA-mediated control and compare it with that of a direct, TF-based transcriptional unit [45]. In principle, since fluctuations can be reduced by increasing the number of molecules, an (almost) arbitrary amount of information can be transmitted through a biochemical network. However, cells have to face the burden of macromolecular synthesis [46-48]. Optimality is therefore the result of a trade-off between the benefits of reduced fluctuations and the drawbacks of the associated metabolic costs. For this reason, we start by fixing a maximal rate of transcription (or, alternatively, the maximal number of output molecules) so as to have a simple but reasonable framework to characterize and compare the capacities of the different regulatory channels. Next, we quantify how an input signal is processed by the transcriptional (TF-based) and post-transcriptional (miRNA-mediated) regulatory elements by characterizing the response in the output ceRNA’s expression levels. In such a setting, information flow is hampered by intrinsic noise if the target gene is weakly derepressed by the activation of its competitor. Otherwise, target derepression appears to have a strong impact on a regulatory element’s capacity. Upon varying the magnitude of derepression by tuning the kinetic parameters, we then show that in certain regimes miRNA-mediated regulation can indeed outperform direct control of gene expression. Finally, we argue that the presence of miRNA molecules in large copy numbers notably reduces the level of intrinsic noise on weakly targeted transcripts. In this case, the mutual regulation of ceRNA molecules by miRNA-mediated channels may become a primary mechanism to finely tune gene expression.

Besides providing a quantitative characterization of the maximal regulatory power achievable through miRNA-based post-transcriptional control, these results provide important hints on the circumstances in which regulation by small RNAs may function as the main tuner of gene expression in cells.

Results

Mathematical model of ceRNA competition

We consider (see Fig 1) a system formed by two ceRNA species (whose levels are labeled m 1 and m 2, respectively) and one miRNA species (with level labeled μ), whose transcriptions are activated by a single TF each (with levels labeled, respectively, f 1, f 2 and f ). Both ceRNAs are in turn targeted by the miRNA. miRNA-ceRNA complexes (levels labeled c with i = 1,2) assemble and disassemble at rates , respectively, whereas complexes can be degraded both stoichiometrically (i.e. without miRNA recycling) at rates σ and catalytically (i.e. with miRNA recycling) at rates κ [49]. In addition, ceRNA and miRNA molecules degrade (resp. synthesize) at rates d (resp. b ) and δ (resp. β), respectively. Steps leading to the formation of the RNA-induced silencing complex (RISC), allowing for the miRNA-ceRNA binding, are neglected for simplicity. TF levels are treated as externally controlled parameters.

Fig 1

Schematic representation of the basic ceRNA network model.

The model comprises three transcription factors (TF1, TF2 and TF) controlling the synthesis of two ceRNA species (ceRNA1 and ceRNA2) and one miRNA species, respectively. The fractional occupancies of TF binding sites are denoted by n 1, n 2 and n , respectively. Both ceRNAs are targets for the miRNA, with whom they form complexes denoted respectively as C1 and C2. The amount of molecules for each species is denoted by the variable next to the corresponding node. Reaction rates are reported next to the corresponding arrow.

Denoting the TF-DNA binding/unbinding rates by k in and k out, respectively (for simplicity, these parameters are taken to be the same for all involved TFs), the TF binding sites’ fractional occupancies n ℓ (0 ≤ n ℓ ≤ 1, ℓ ∈ {1, 2, μ}) obey the dynamics according to which transcriptional activation requires the cooperative binding of h TF molecules for each RNA species involved. In general, the occupancies by different TFs equilibrate on different timescales [50]. However it is often assumed that the transcriptional on/off dynamics is much faster than transcription itself [51, 52]. As a consequence, each n ℓ can be fixed at its ‘equilibrium’ value.

Every process in the above scheme contributes to the overall level of noise. We represent the mass-action kinetics of the system through the set of coupled Langevin processes (i = 1,2) [41, 53] where the mutually independent random (Poisson) ‘forces’ ξ , , ξ , and denote, respectively, the intrinsic noise in ceRNA levels (due to random synthesis and degradation events), in the association/dissociation processes of complexes, in the miRNA level, in the catalytic complex decay and in the stoichiometric complex decay. Each of the above noise terms has zero mean, while correlations are given by where we introduced the steady state molecule numbers

We shall be interested in the fluctuations of molecular levels around the steady state that arise due to intrinsic noise sources. The stochastic dynamics of the system can be simulated via the Gillespie algorithm (GA, see Materials and Methods). Fig 2A shows typical GA results for , and , with the corresponding Fano factors (FFs) shown in Fig 2B, as functions of f 1. Numerical results are matched against analytical estimates obtained by the linear noise approximation (see Materials and Methods).

Fig 2

Cross-talk scenario in the ceRNA network.

Steady state values of (A) the ceRNA and miRNA levels, and (B) the corresponding Fano factors. Markers denote GA computations, curves correspond to approximate analytical solutions obtained by the linear noise approximation. Values of the kinetic parameters are reported in Table 1.

Table 1

Parameters used to obtain each figure.

When the parameter was varied over a range of values, the range is shown in brackets. The parameter ω appearing in the last column denotes the re-scaling factor defined in Eq (24), while the channel to which parameters refer is specified in brackets. Values of f 1, f 2 and f are mapped to those of via the steady state conditions, Eq (2).

Parameter	Fig 2A and 2B	Fig 5A-5C	Fig 6A-6C	Fig 7A	Fig 7B	Fig 8	Fig 9
b 1 [molecule min−1]	120	100	100	100	100	[1, 147] (miRNA); 0 (TF)	90 (miRNA); 0 (TF)
b 2 [molecule min−1]	100	100	100	100	100	110.2 (miRNA); [10, 98] (TF)	22.3 (miRNA); 22(TF)
β [molecule min−1]	100	100	100	100	100	100 (miRNA); 0 (TF)	80 (miRNA); 0(TF)
d 1 [min−1]	0.1	0.1	0.1	0.1	0.1	0.1	0.1
d 2 [min−1]	0.1	0.1	0.1	0.1	0.1	0.1	0.1
δ [min−1]	0.1	0.1	0.1	0.1	0.1	0.1	0.1 ω
σ 1 [min−1]	0.5	1	1	1	1	1	1
σ 2 [min−1]	1	1	1	1	1	1	1
k1+ [molecule−1 min−1]	0.002	[e −10, e 2]	e 2	e 2	e 2	e 3 (miRNA); 0 (TF)	e 3.39 ω (miRNA); 0(TF)
k2+ [molecule−1 min−1]	0.001	[e −10, e 2]	e −3	e −3	e −3	e 0(miRNA); 0 (TF)	e −5.77 ω (miRNA); 0(TF)
k1- [min−1]	0.0005	0.001	0.001	0.001	0.001	0.001	0.001
k2- [min−1]	0.001	0.001	0.001	0.001	0.001	0.001	0.001
κ 1 [min−1]	0.0003	0.002	[e −5, e 5]	0.002	e −8	0.002	0.002
κ 2 [min−1]	0.0004	0.002	[e −5, e 5]	0.002	e 5	0.002	0.002
n¯_1 [adim.]	(0,0.99)	(0,0.99) (miRNA); 0.99(TF)	(0, 0.99) (miRNA); 0.21(TF)	(0.1,0.9)	(0.1,0.9)	(0, 0.99) (miRNA)	(0, 0.99) (miRNA)
n¯_2 [adim.]	0.71	0.99 (miRNA); (0,0.99) (TF)	0.21 (miRNA); (0, 0.99) (TF)	(0.1,0.9)	(0.1,0.9)	0.99 (miRNA); (0, 0.99) (TF)	0.99 (miRNA); (0, 0.99) (TF)
n¯_μ [adim.]	0.9	0.99	0.99	0.99	0.5	0.99	0.99
h [adim.]	5	5	5	5	5	5	5
k out/k in [moleculesh]	63300	63300	63300	63300	63300	63300	63300
f max [molecules]	30	30	30	30	30	30	30

In particular, in Fig 2 one sees that an increase of is accompanied by a concomitant increase of and by a decrease of the average number of free miRNA molecules. This is an instance of miRNA-mediated ceRNA cross-talk. Indeed, upon up-regulating m 1 by injecting the corresponding TF (i.e. by increasing f 1), the level of free miRNAs will decrease as more and more molecules will be actively repressing ceRNA1, causing in turn an up-regulation of ceRNA2. f 1 will thus positively correlate with m 2. One easily sees that steady-state ceRNA levels depend on through a sigmoidal function, namely [42] The constant acts as a ‘soft’ threshold for the miRNA level, allowing to distinguish three situations:

if , ceRNA is free or unrepressed: spontaneous ceRNA degradation dominates over miRNA-mediated decay channels, so that, effectively, the ceRNA level is weakly sensitive to small changes in ;

if , ceRNA is bound or repressed: miRNA-mediated ceRNA decay dominates over spontaneous ceRNA degradation but, again, the ceRNA level is weakly sensitive to small changes in as most ceRNAs are bound in complexes with the miRNA;

if , ceRNA is susceptible to μ: spontaneous and miRNA-mediated decay channels have comparable weights and the ceRNA level is very sentitive to small changes in .

The behaviour of the FFs (see Fig 2B) emphasizes how noise patterns change in the different regimes. In the displayed example, ceRNA1 and ceRNA2 become susceptible for ln(f 1) ≃ 2.5 and ln(f 1) ≃ 2.1, respectively. Indeed, one observes that the corresponding FFs peak close to these values, in accordance with the observation that stochastic fluctuations are enhanced when the rates of substrate supply are adequately balanced in a stoichiometrically coupled system [34, 54, 55]. The FF for ceRNA2 appears to approach one for large values of f 1, as expected for the pure Poisson birth/death process that characterizes the free regime [56]. On the other hand, for very small but nonzero mean fractional occupancy (corresponding to small values of f 1), m 1 will with high probability only take on the values 0 or 1, as a transcribed molecule will quickly undergo degradation or sequestration in a complex. In such a situation, the mean and variance of m 1 will be calculated by summing up zeros and ones over time, leading to a FF equal to one.

Measuring the performance of regulatory elements

The amplitude of ceRNA derepression

We shall focus on quantifying the regulatory capacity of the post-transcriptional and transcriptional channels that link, respectively, TF1 to ceRNA2 and TF2 to ceRNA2, as depicted in Fig 3.

Fig 3

Schematic representation of the channels under consideration.

(A) miRNA-channel, (B) TF-channel.

In what follows, we shall take the number of transcription factors as the key input variables (f 1 for the miRNA-mediated channel and f 2 for the direct channel), and the number of target RNAs m 2 as the key output variable. It is reasonable to think that an effective regulator should have a sensible impact on the target molecule. As a first measure of regulatory performance we will therefore consider the amplitude of variation (AOV) of m via f , defined as corresponding to the difference between the maximum and the minimum steady state level of ceRNA obtained by varying f from to . Clearly, larger values of Δ imply stronger degrees of derepression of ceRNA induced by the corresponding TF. For sakes of simplicity, we shall adopt the following notation to distinguish the AOV induced by transciptional and miRNA-mediated control channels: Note that the AOV only involves mean levels. As such, it is unable to account for the intrinsic degree of stochasticity of the system.

Maximal mutual information

In order to evaluate the impact of noise, as a further quantifier of regulatory power we shall employ the mutual information (MI), a standard information-theoretic measure of channel performance given by (i, j = 1, 2) [52] where p(m , f ) stands for the joint probability distribution of m and f , while p(m ) and p(f ) denote its marginals. The MI is a non-negative quantity that vanishes when m and f are independent. Conversely, when one can reliably distinguish (roughly) different values of the output variable upon tuning the input (see [9, 45, 52] for further details about the usefulness of this scheme for the analysis of transcriptional regulatory or biochemical modules). The maximal achievable information flow for a fixed channel (i.e. for given kinetic parameters) is called the channel capacity and can be computed by maximising the MI over the input variable (in this case, over the distribution of the TF level). In the simplest case, one can assume that the channel is Gaussian, meaning that for every given input f the output m 2 is described by the conditional probability density where both the average (f ) and variance depend on the input variable f . In the limit of small variance, the distribution p opt(f ) of TF levels that maximises the MI between m 2 and f reads [45] where is a normalizing factor. For p = p opt one obtains the channel capacity S1 Fig provides a qualitative yet intuitive representation of how the shape of the input/output relation and the size of fluctuations affect information flow. Intuitively, the normalizing factor Z ‘counts’ the number of distinguishable output levels comparing the local slope of the input/output curve with the noise strength: if the slope is too low or the fluctuations are too high no information can be transmitted through the channel. We shall use the following notation to distinguish the capacities of the direct and indirect channels: We shall furthermore denote their difference by ΔI = I TF − I miRNA.

We are henceforth going to adopt the following protocol. After generating the output data () by simulating the circuit’s dynamics via the GA, the optimal input distribution will be constructed according to Eq (19). We shall then let the circuit process an input signal sampled from that distribution and calculate the MI for both the transcriptional and post-transcriptional channels (see the flowchart in Fig 4). Channel capacities will then be compared, with the goal of clarifying the conditions under which the former can be more effective than the latter.

Fig 4

The flowchart of the method.

A C++ implementation is available at https://github.com/araksm/ceRNA.

Comparing indirect miRNA-mediated regulation with direct transcriptional control

Target prediction algorithms suggest that miRNA binding affinities vary significantly among their targets [57]. Such heterogeneities indeed are mapped to the miRNA binding kinetics and are shown to influence the susceptibility of the targets to the miRNA molecules [36]. Moreover, the level of complementarity between the regulator and the target seems to be decisive for the selection of a decay channel (catalytic or stoichiometric) for the miRNA-ceRNA complex [15, 27]. One may therefore expect that the effectiveness of miRNA-mediated post-transcriptional control depends strongly on the kinetic parameters characterizing the network.

Here in particular, we are going to investigate how the capacities of these regulatory elements are affected by changes in (i) miRNA-ceRNA binding kinetics, (ii) miRNA recycling rates, and (iii) effective transcription rates of all RNA species involved.

In order to contrast the performances of miRNA- and TF-channels we shall start by analyzing their respective responses to the same input signal. More precisely, the input variable f will be varied from 0 to a value f max defined by the condition n (f max) = 0.99 (i.e. from a situation in which the promoter is always free to one in which it is essentially always occupied).

The capacity of the miRNA-channel is maximal in a specific range of miRNA-ceRNA binding rates

Available estimates of RNA binding energetics indicate that the target sites of an individual small non coding RNA can span a wide range of affinity both in mammalian cells [57] and in bacteria [58]. To explore the functional consequences of such heterogeneities we analyse the behaviour of the regulatory channels changing the values of the complex association rates and while keeping fixed all others parameters. As displayed in Fig 5, the maximal amount of information transmitted through the post-transcriptional channel (Fig 5A) increases monotonically with (for fixed ) and displays a maximum versus (for fixed ). Indeed, the higher , the higher c 1, and the more information about changes in f 1 can be transmitted through the miRNA-channel. On the other hand, ceRNA2 is unrepressed by the miRNA for small values of while it is strongly repressed for high values of , for any value of the input. In both regimes, ceRNA2 is only weakly sensitive to changes in the level of its competitor. In such conditions, no information can be transmitted through the miRNA-channel. The information flow is therefore maximal between these two extremes.

Fig 5

Dependence of the transcriptional and post-transcriptional channel capacities on miRNA-ceRNA association rates.

(A) I miRNA. (B) I TF. (C) ΔI = I TF − I miRNA. Values of the kinetic parameters are reported in Table 1.

For the transcriptional channel, we expect to see maximum information transmission when the target, ceRNA2, is unrepressed by the miRNA, in which case the channel is effectively decoupled from the rest of the network, while no information transmission is expected when the target is fully repressed by the miRNA. Indeed, in Fig 5B one sees that, for small values of (unrepressed target), I TF is the highest and it decreases to zero as increases together with the degree of target repression. Moreover, since a higher moves m 2 to larger values (i.e. towards the free regime), I TF increases with and decreases with in the susceptible regime. Notice that, when post-transcriptional capacity is highest, miRNA-mediated regulation is as effective as transcriptional control (see Fig 5C).

The intuition that the channel capacity is strongly coupled to the target’s size of derepression (its AOV) is confirmed by the fact that the latter indeed displays a very similar behaviour (see S2 Fig).

miRNA-mediated regulation may represent the sole control mechanism in case of differential complex processing

Direct experimental measurements of miRNA-recycling rates are challenging [59]. However it is well known that the degreee of complementarity of a miRNA to its target lead to different repression pathways [60] and that differential complex processing strongly affects reciprocity of competition [33]. In order to better quantify the degree of regulatory control achievable in case of different processing mechanisms we changed the values of miRNA-recycling rates κ 1 and κ 2 keeping fixed all the other parameters, obtaining the results displayed in Fig 6. Both I TF and I miRNA are higher for small values of the recycling rates, and decrease as miRNA recycling increases. The observed dependence of I miRNA on κ 1 agrees with the results of [42], where the sensitivity of miRNA molecules to changes in the levels of their targets decreases as the miRNA recycling rate increases, eventually vanishing for κ 1 ≫ σ 1. Indeed, for large κ 1 miRNAs become insensitive to changes in f 1, and hence no information can be transmitted.

Fig 6

Dependence of the transcriptional and post-transcriptional channel capacities on miRNA recycling rates.

(A) I miRNA. (B) I TF. (C) ΔI = I TF − I miRNA. Values of the kinetic parameters are reported in Table 1.

Remarkably (see Fig 6C), for large κ 2 and small κ 1 miRNA-mediated regulation may be the sole mechanism able to effectively control the target level. In this case, in fact, I TF becomes especially small as the target is completely repressed by the miRNA while I miRNA stays finite since the activation of ceRNA1 can titrate away the miRNA and therefore derepress m 2.

As for the previous case, the behaviour of the AOV correlates strongly with this scenario, i.e. larger AOV corresponds to higher achievable information flow in the regulatory element (see S2 Fig). This suggests that, although targets that are degraded purely catalytically cannot compete for miRNA at steady state (as also shown for bacterial small RNAs [61]), they can be notably derepressed if one of their competitors decays stoichiometrically.

Effective derepression of a repressed target may be conditional on the activation of its competitor and of miRNA

Transcriptional and post-transcriptional channel capacities are expected to depend on the effective transcription rates of the ceRNAs and the miRNA, described respectively by the quantities and , since these bear a direct impact on the regime (repressed, susceptible or unrepressed) to which each ceRNA belongs. In particular, upon changing the fractional occupancies, situations where miRNA-mediated regulation outperforms TF-based control can occur. Such a possibility might indeed be realized for relatively large values of , where the miRNA-mediated channel works optimally. This regime is explored in Fig 7.

Fig 7

Dependence of ΔI on the fractional occupancy of the TF binding site in the optimality range for the post-transcriptional miRNA-mediated channel.

(A) Case of targets with weak catalytic degradation (small κ 1 and κ 2). (B) Case of a strongly catalytically degraded target (large κ 2) with a weakly catalytically degraded competitor (small κ 1). Note that for (A) while for (B). Values of the kinetic parameters are reported in Table 1.

Remarkably, one finds that post-transcriptional control outperforms transcriptional one when

the effective trascription rate of ceRNA2 exceeds that of ceRNA1 (, see in Table 1) in presence of weak catalytic degradation (Fig 7A);

the effective transcription rate of ceRNA1 is smaller than that of the miRNA (, see in Table 1) in presence of strongly catalytically degraded target. (see Fig 7B).

A more detailed analysis shows that, in these scenarios, whenever I miRNA > I TF the target’s derepression size is larger when the signal is processed post-transcriptionally rather than through the direct transcriptional channel, as shown in S3 Fig. Again, the AOV appears therefore to be a good proxy for the channel capacity.

Previous work had pointed to near-equimolarity of ceRNAs and miRNAs within a network as one of the central preconditions for significant ceRNA regulation [34, 42, 62]. On the other hand, experiments on liver cells have shown that such conditions may not be met in physiological conditions, where the number of miRNA target sites can vastly exceed that of miRNA molecules [29, 30]. Our results suggest that target derepression may be significant even if the competitor is in low copy numbers, provided a certain heterogeneity in kinetic parameters (e.g. for a catalytically degraded target and a stoichiometrically degraded competitor) is present, in line with recent experiments pointing to a key role of parameter diverseness [36].

A sufficiently high degree of target derepression is required for information transmission

We have seen that AOV displays a dependence on kinetic parameters that very closely matches the one found for the channel capacities. In particular, whenever I miRNA > I TF, we found ΔmiRNA > ΔTF. One may be lead to think that the AOV by itself may fully describe the capacity of miRNA-mediated post-transcriptional information processing and that stochasticity could, to a large degree, be neglected. In some situations, however, the size of the output range is poorly informative about the channel performance. This becomes clear by exploring the relationship between the AOV and the channel capacity at a quantitative level. In Fig 8 the dependence of the optimal MI on the AOV is discussed for two cases: (i) direct transcriptional regulation by TF in absence of miRNA molecules (black curve), (ii) miRNA-mediated post-transcriptional regulation (blue curve).

Fig 8

Channel capacities as a function of the target’s degree of derepression (AOV).

The black curve corresponds to direct transcriptional control in absence of miRNAs, while the blue curve describes the behavior of the miRNA-mediated post-transcriptional channel. For both channels . The predicted maximal MI in the Poissonian limit given in S1 Text is shown as a dashed line. Values of the kinetic parameters are reported in Table 1.

In Fig 8 we indeed see that for large enough AOV, where fluctuations on each level are small compared to ΔTF, the TF-channel capacity approaches the theoretical limit of a Poissonian channel worked out in S1 Text. On the other hand, for small values of ΔTF, the capacity of TF-channel becomes vanishingly small.

Similarly in Fig 8 one can see that I miRNA increases (roughly) logarithmically with ΔmiRNA if the latter exceeds a certain threshold, whereas it is essentially zero below this value. In other words, the miRNA-mediated channel is unable to convey any information unless the degree of target derepression is sufficiently high. Hence I miRNA can be close to zero even though ΔmiRNA is finite. Notice that the derepression threshold for miRNA-mediated information transmission is larger than the corresponding threshold for direct transcriptional information flow.

A close analysis reveals that such a feature is due to the random fluctuations on the lowest expression level of the target, . Below threshold, in fact, the degree of derepression of the target is smaller than the fluctuations on caused by intrinsic noise (reported in Fig 8). Therefore the response observed at the level of the output is indistinguishable from the noise. In the miRNA-mediated post-transcriptional channel the threshold AOV is larger since noise on is higher. As soon as the AOV overcomes random fluctuations on , the channel capacities starts to increase with ΔmiRNA.

Noise reduction makes the capacity of the post-transcriptional channel comparable to that of the transcriptional channel in the limit of large populations of weakly interacting miRNAs

We will now focus specifically on the influence of the miRNA-ceRNA binding component of the intrinsic noise, which will be studied upon varying miRNA-ceRNA binding rates and miRNA population size. We shall consider two different scenarios, namely (a) direct, TF-based transcriptional regulation in absence of miRNAs (corresponding to for i = 1,2), and (b) indirect post-transcriptional miRNA-mediated regulation (corresponding to non-zero complex binding rates). Kinetic parameters will however be re-scaled as where ω > 0 is the re-scaling factor. Note that the above choice does not change as long as and (both of which are true in the case we consider). By changing ω one may therefore characterize the dependence of the channel capacities on the kinetic parameters δ and while keeping the output range approximately fixed, such that ΔmiRNA ≃ ΔTF for every ω.

Fig 9A shows how the the channel capacities change with ω. One sees that in general I miRNA < I TF, so that at fixed AOV the miRNA-mediated channel typically has smaller capacity than the direct one. However, in the limit of small ω, i.e. when miRNA populations are sufficiently large and miRNA-ceRNA couplings are weak, the miRNA channel can perform as effectively as a pure transcriptional channel. This result can be understood by examining the noise levels (see Fig 9B). For small ω (or high μ 0,2), in particular, the FF approaches the Poissonian limit since (i) for μ ≳ μ 02 relative fluctuations in miRNA levels are small, extra noise coming from molecular titration disappears and the target degradation constant is rescaled as ; (ii) μ ≪ μ 0,2 target ceRNA is seen in free regime, where it is not sensitive to miRNA molecules. In such conditions, miRNA-mediated noise buffering [20] cannot be observed, nevertheless transcriptional and miRNA-mediated regulation are effectively equivalent.

Fig 9

Comparison of transcriptional and post-transcriptional regulation for a fixed output variation range.

(A) Dependence of I TF and I miRNA on ω. (B) FF of m 2 at the corresponding steady state. AOV and are the same for both channels, namely ΔmiRNA = ΔTF = 209 ± 1 and . Values of the kinetic parameters are reported in Table 1.

Observed heterogeneities in miRNA-target interaction strength may enable regulation by miRNA-mediated cross-talk

The quality of miRNA-target interaction influences the binding kinetics and may be decisive for the activation of the target decay channel [15, 27, 63, 64]. Estimations of the miRNA-ceRNA binding affinities are experimentally challenging. However, computational methods predict a considerable degree of heterogeneity across different miRNA-ceRNA pairs [57]. Remarkably, for the majority of cases reported in the literature, the predicted miRNA-binding energies of the RNAs on the ‘input side’ of the channel are lower than those of the RNAs on the ‘output side’, in line with the optimal conditions highlighted by our model. For example, binding affinities between the long non-coding RNA linc-MD1 and its regulatory miRNAs (playing a central role in skeletal muscle cell differentiation) are significantly lower than those characterizing miRNA interactions with linc-MD1’s competitors, namely the MAML1 and MEF2C mRNAs, as predicted by miRanda algorithm [39, 64]. Likewise, the circular RNA CDR1 has been found to contain around 70 binding sites with high complementarity for miR-7, corresponding to a strong effective coupling through which it can regulate the expression of miR-7’s target genes [37, 65, 66]. Finally, the high sequence homology of pseudogenes (long non-coding RNA genes developed from protein-coding genes but unable to produce proteins) with their parental gene allows them to compete for a large number of shared miRNAs [28, 38, 67].

Our study also points to the potential relevance of miRNA-ceRNA complex decay channels. It is known that, in case of sufficient complementarity, miRNA molecules can function as siRNAs and cleave their targets after binding them [15, 27]. Such targets decay catalytically and are therefore effectively degraded by the miRNAs. Our model predicts that miRNA-mediated control may be the preferred regulatory mechanism in presence of kinetic heterogeneities at the level of miRNA recycling rates. In particular, targets that undergo catalytic degradation may be efficiently derepressed by their competitors. This type of scenario has been observed in experiments concerning the bacterial small RNA Qrr [61]. Qrr represses its targets by distinct mechanisms. For instance, luxR is repressed catalytically, luxM stoichiometrically, while luxO is silenced through translational repression. luxM and luxO are however able to derepress LuxR in the presence of Qrr. In the light of our results, these observations may therefore point to a higher than expected role for the ceRNA effect in vivo, especially in cases in which heterogeneities in kinetic parameters are thought to be strong [34, 59, 62, 68, 69].

Discussion

Non-coding RNA molecules, and miRNAs specifically, are increasingly associated to regulatory functions. Besides making it mandatory to characterize the specific role of ncRNAs on a case-by-case basis, especially for situations like disease or differentiation, this fact also raises the question of what ingredients can make miRNAs a preferred tool to regulate the level of a target RNA over, for instance, the target’s TF. A possible answer lies in the noise-buffering role that miRNAs can play, which is especially evident in genetic circuitries like incoherent feed-forward loops [19, 20]. By reducing relative fluctuations in the output level, miRNAs can confer robustness to gene expression profiles. However the so-called ‘ceRNA hypothesis’ opens the way to the possibility that their regulatory functions are carried out at a broader, though more subtle, level. In short, according to the ceRNA scenario miRNAs can mediate an effective positive interaction between their target RNAs driven by the targets’ competition to bind them. In this sense, miRNAs can be seen as a sort of ‘channel of communication’ between RNAs through which RNA levels can be altered and noise can be processed (both buffered and amplified). Previous work [42] has shown that the ceRNA effect may generate both highly plastic and highly selective ceRNA-ceRNA couplings, thereby representing a potentially powerful mechanism to implement gene regulation at the post-transcriptional level.

Although predicted theoretically, the extent and relevance of ceRNA effect in vivo is poorly understood. On one hand, considerable evidence points to the ceRNA effect playing a major role in certain dis-regulated or transient cellular states. For instance, it has been shown that the expression of the tumor-suppressor gene PTEN can be regulated by its miRNA-mediated competitors VAPA, CNOT6L, SERINC1 or ZNF460 [70]. Furthermore, many pseudogenes have been found to compete with their parental genes for a shared pool of common microRNAs, thus regulating their expression as competitive endogenous RNA [28, 38, 67, 71]. Such mechanisms seem to be of particular relevance in cancer. For instance, murine models engineered to overexpress the pseudogenes of the proto-oncogene BRAF develop an aggressive malignancy resembling human B cell lymphoma since, by functioning as ceRNAs, they elevate BRAF expression both in vitro and in vivo [38]. Likewise, the long noncoding RNA linc-MD1 has been shown to regulate the skeletal muscle cell differentiation clock by sponging miRNAs from its competitors, thereby enacting a ceRNA mechanism. In particular, MAML1 and MEF2C (coding for transcription factors that activate muscle-specific gene expression) compete with linc-MD1 for miR-133 and miR-135 respectively [39]. Taken together, the available evidence indicates that miRNA activity depends on the miRNA:target ratio, on miRNA target site abundance and on miRNA binding affinities. Further analyses of high throughput datasets confirm this observation [30, 36]. One may therefore question how said factors may influence miRNA-mediated post-transcriptional control.

The problem however arises of quantifying the degree of control that can be exerted through miRNAs. Taking the ‘channel’ analogy more strictly (as done before for simpler regulatory elements [9, 45, 52]), one may resort to information theoretic concepts and tools to characterize precisely how well miRNAs can process fluctuations coming from the modulator nodes and transfer them to the target nodes. This issue goes beyond noise buffering, specifically including the ability to respond to large changes in the mean levels as well as to changes in the structure of fluctuations. As the properties of a channel are conveniently encoded in the mutual information between the input and output nodes, asking how well a channel can function amounts to asking what is the channel’s capacity, i.e. the maximum value of the input-output (or modulator-target) mutual information achievable through that channel. This work aimed precisely at quantifying the effectiveness of microRNA-mediated post-transcriptional control of gene expression by computing the capacity of the corresponding regulatory channel and comparing it to that of direct, TF-driven transcriptional regulation. Evidently, multiple factors can influence the flow of information across nodes in a biochemical network, starting from the intrinsic noisiness of each reaction step.

Our basic challenge was therefore understanding in which circumstances miRNA-mediated control can outperform the TF-based one, thereby obtaining insight on why the ceRNA effect appears to be so often employed by cells in situations where accurate tuning and/or shifts of expression levels are required. We have therefore considered, along the lines of [33, 34, 41–44], a mathematical model of the ceRNA effect and characterized its steady state in terms of both mean molecular levels and regulatory capacities of the miRNA-mediated and TF-based channels via stochastic simulations.

We have first considered how the two channels process inputs (the TF levels f 1 and f 2) that vary in the same range. We have shown that, while the capacity of the TF-channel depends monotonously on each miRNA-ceRNA binding rate and is largest when the target is unrepressed by miRNAs (as might have been expected), the capacity of the post-transcriptional channel is maximal in a specific range of values of the miRNA-ceRNA binding rates. In agreement with [33], we found that miRNA-channel’s efficiency is tunable to optimality by the binding kinetics. Furthermore, our model suggests that both capacities decrease as the miRNA recycling rates increase, confirming previous indications obtained by different analytical techniques [42]. Consistently with the scenario observed experimentally for the bacterial small RNA Qrr [61], our model finally suggests that catalytically regulated targets are weakly capable of competing for miRNAs but might be significantly derepressed by their competitors.

In addition, post-transcriptional miRNA-mediated information processing was shown to be characterized by a threshold behaviour as a function of the AOV. In other terms, no information can be transmitted across the channel unless the target’s degree of derepression is sufficiently large. This implies that the regulatory effectiveness of the channel is well encoded by the degree of target derepression when the latter is sufficiently high, in which case it is possible to identify regimes in which post-transcriptional regulation is more accurate than transcriptional control.

To get a deeper insight on the origin of the observed threshold behaviour one must however go beyond the AOV and consider more carefully how the miRNA-ceRNA binding noise affects the overall picture. After showing that miRNA-ceRNA binding noise is indeed at the origin of the threshold behaviour that limits the miRNA-channel capacity, we have uncovered the rather remarkable property that in presence of large but weakly interacting miRNA populations the ceRNA effect can regulate gene expression as effectively as the target’s modulator node itself.

The present work has focused on a small genetic circuit made up of a single miRNA species and two target RNA species at steady state. Previous work has however shown that cross-talk is possible even during transients [44]. Going beyond stationarity is therefore likely to bring to light new scenarios where miRNA-mediated regulation plays possibly a yet more prominent role.

On the other hand in a typical eukaryotic cell there are thousands ceRNAs, hundreds miRNAs and a rich structure of conserved targeting patterns [72]. Moreover, cells might be interested in tightly controlling not only each output individually but also particular combinations of output levels (which might be required e.g. for the efficient operation of metabolic pathways). In such a scenario, miRNA-mediated control could represent a powerful mechanism to increase robustness and flexibility in specific directions of the output space. In view of this, it would be important to consider a more general multi-source network coding problem in which a large number of transcription processes are seen as mutually independent information sources, and each of the information sources is multicast to sets of output nodes through the effective network of miRNA-mediated cross-talk interactions. The information-theoretic scheme employed in this work is easily generalized to deal with more complex networked situations. Novel insight might finally shed light on the partly controversial picture unveiled by recent experiments addressing the relevance of the ceRNA effect in vivo [30, 35–39].

Materials and Methods

Gillespie algorithm

Numerical simulations have been carried out using the Gillespie algorithm (GA), a standard stochastic method to analyze the time evolution of a system of chemical reactions which is exact for spatially homogeneous systems [73]. In short, based on the reaction rates, GA calculates when the next reaction will occur and what reaction it will be, and then modifies the amount of each molecular species in the system according to the process that took place. If we denote the probability of reaction r to occur in the time interval (τ, τ + dτ) by P(r, τ)dτ, the algorithm proceeds through the following steps:

Initiate the number of reactants in the system and the termination time;

Generate a random pair (r, τ) according to P(r, τ);

Using the pair (r, τ) just generated, advance time by τ and change number of species involved in reaction r accordingly;

Read out the molecular population values. If the termination time is reached, stop the simulations, otherwise return to Step 2.

After a long run, independently on the initial setup, the system of chemical reactants will come to the equilibrium state.

Linear noise approximation

The mathematical model of ceRNA competition can be solved numerically in the so-called linear noise approximation. Letting x = (m 1, m 2, μ, c 1, c 2) stand for the vector of molecular levels, the kinetic mass action Eq (3) can be re-cast in compact form as where the vector g encodes for the deterministic part of the dynamics, while the vector η represents the aggregate noise terms. Each element of η has zero mean, and we shall denote its correlations by 〈η (t)η (t′)〉 = Γ δ(t − t′). Denoting by the steady state, small deviations from it (i.e. ) relax, in the linear regime, according to where . In this approximation, the correlation matrix C = 〈δx δx 〉 is given by [74] where λ’s and B’s are, respectively, eigenvalues and eigenvectors of the matrix A.

Model parameters

Table 1 reports the values of the kinetic parameters (or of their range of variability) employed in the different figures. As and d are varied in wide ranges in order to test how channel capacities depend on them, our choice was mainly guided by the need to focus the analysis on regimes where ceRNA cross-talk is established (so that the miRNA-mediated channel can actually convey information).

Information transmission and noise.

Qualitative depiction of how noise affects information flow. Three different situations for the input-output curve m 2 (target) vs f (modulator) are shown, namely, long (A) target expression level independent of the modulator; (B) derepressed target, with large fluctuations; (C) derepressed target, with small fluctuations. In the first two cases, little or no information can be transmitted from modulator to target (hence the regulatory effectiveness is severely limited), as either the target is insensitive to changes in modulator levels, or its response is strongly hindered by noise. On the other hand, in case (C) information will be transmitted, since modulating the input one can clearly distinguish different output levels. The number of distinguishable levels is linked to the local slope of the input/output curve, as shown mathematically in Eqs (20) and (21) and is mainly limited by the noise strength. In the limit of vanishing noise, when the input-output relationship becomes deterministic, the mutual information between m 2 and f diverges.

(TIF)

Click here for additional data file.

Dependence of the AOV on kinetic parameters.

(A) ln ΔmiRNA, (B) ln ΔTF, (C) ln ΔTF—ln ΔmiRNA. Values of the kinetic parameters are as in Fig 5 for the panels in the top row and as in Fig 6 for the panels in the bottom row.

(TIF)

Click here for additional data file.

Dependence of ln ΔTF − ln ΔmiRNA on the fractional occupancy of the TF binding site.

Values of the kinetic parameters are as in Fig 7A for panel (A), and as in Fig 7B for panel (B).

(TIF)

Click here for additional data file.

Capacity of a Poissonian channel.

(PDF)

Click here for additional data file.

Dataset used in the figures.

(ODS)

Click here for additional data file.