Competing endogenous RNA crosstalk at system level

microRNAs (miRNAs) regulate gene expression at post-transcriptional level by repressing target RNA molecules. Competition to bind miRNAs tends in turn to correlate their targets, establishing effective RNA-RNA interactions that can influence expression levels, buffer fluctuations and promote signal propagation. Such a potential has been characterized mathematically for small motifs both at steady state and \red{away from stationarity}. Experimental evidence, on the other hand, suggests that competing endogenous RNA (ceRNA) crosstalk is rather weak. Extended miRNA-RNA networks could however favour the integration of many crosstalk interactions, leading to significant large-scale effects in spite of the weakness of individual links. To clarify the extent to which crosstalk is sustained by the miRNA interactome, we have studied its emergent systemic features in silico in large-scale miRNA-RNA network reconstructions. We show that, although generically weak, system-level crosstalk patterns (i) are enhanced by transcriptional heterogeneities, (ii) can achieve high-intensity even for RNAs that are not co-regulated, (iii) are robust to variability in transcription rates, and (iv) are significantly non-local, i.e. correlate weakly with miRNA-RNA interaction parameters. Furthermore, RNA levels are generically more stable when crosstalk is strongest. As some of these features appear to be encoded in the network's topology, crosstalk may functionally be favoured by natural selection. These results suggest that, besides their repressive role, miRNAs mediate a weak but resilient and context-independent network of cross-regulatory interactions that interconnect the transcriptome, stabilize expression levels and support system-level responses.


INTRODUCTION
Competition to bind substrates, enzymes or gene expression machinery is ubiquitous in biological networks and impacts regulatory processes in several ways [1][2][3][4][5][6][7][8][9][10][11][12]. For instance, the initiation and translation rates of different transcripts are effectively coupled by the competition for the ribosome pool, so that modifications of a given RNA species can alter the translational dynamics of other RNAs [13]. Quite generally, competition for limited and shared molecular resources induces effective interactions between the competing species, with signs (positive or negative) that depend on the specifics of the underlying processes [14]. While such interactions constitute in principle an additional layer of indirect regulation, their intensity is strongly context-dependent [9]. The functional role of competition-driven crosstalk therefore has to be evaluated on a case-by-case basis.
Competition for miRNAs (or, more generally, small regulatory RNAs) among long transcripts is undergoing much * andrea.demartino@roma1.infn.it scrutiny in this respect [15]. In silico studies of small motifs, summarized e.g. in [16,17], have characterized how the strength, selectivity and directionality of miRNA-mediated RNA crosstalk are modulated by kinetic and topologic ingredients, leading to highly adjustable output profiles [18][19][20], differential processing of intrinsic and extrinsic heterogeneities [21][22][23][24], stabilization of protein levels [25,26] and long-range effects [27], both at steady state and during transients [28]. Experimental evidence, however, suggests that, in order to fully develop its potential, RNA crosstalk presupposes rather specific conditions, either in terms of the size of the perturbation required to generate a significant response [29][30][31] or in terms of molecular abundances and kinetic parameters [32][33][34] (see e.g. [35,36] for reviews). When such conditions are not met, miRNAs only provide a weak coupling channel for RNAs.
Generally speaking, weak individual crosstalk interactions by themselves do not necessarily imply a reduced physiological role. This is especially true in large networks, where many interactions can aggregate and perturbations can propagate by exploiting topological and kinetic heterogeneities [18,22,37]. On-going explorations of miRNA-RNA networks are indeed uncovering a high degree of hard-wired complexity [38,39].
In the light of these studies, achieving a better understanding of RNA crosstalk from a systemic perspective has become a pressing issue.
Our goal here is to examine RNA crosstalk in silico in extended miRNA-RNA interactomes as a function of various parameters, including global miRNA levels, degrees of parameter heterogeneity, and topological characteristics of the networks. To cope with the lack of knowledge about kinetic parameters, we make use of a maximum-entropy assumption [40]. In short, after obtaining the steady states of the miRNA-RNA network in terms of a small number of kinetic parameters, we focus on the statistics of different quantities induced by a probability distribution of these parameters. This allows to extract context-independent, or typical, features at the cost of weakening our ability to make predictions for individual crosstalk interactions.
In short, our main results can be summarized as follows.
1. Although typically weak, the emergent crosstalk structure is a robust feature of the miRNA-RNA network; for instance, its mean intensity is modulated by miRNA levels but is otherwise weakly affected by transcriptional and/or kinetic heterogeneities (including binding affinities).
2. On the other hand, variability in transcription rates generically enhances the maximal crosstalk intensity achievable as well as non-local effects (i.e. the emergence of long-range crosstalk mediated by chains of miRNA-RNA interactions).
3. The stability of expression profiles is generically higher when crosstalk is strongest. 4. The degrees of RNA and miRNA nodes are the key topological controllers of the above picture.
Overall, these points suggest that miRNA-RNA networks encode for complex and adaptive crosstalk patterns that feed back on the stability of expression profiles despite the fact that the typical crosstalk link is very weak. A relatively small number of stronger couplings drives this scenario, while transcriptional and topologic heterogeneities allow to extend the range of crosstalk up to network scale.

Mathematical model
To model a network comprising M miRNA species and N RNA species we have extended the mathematical framework employed in [18] for the study of small motifs. Conforming to experimental evidence according to which mature miR-NAs are mostly bound to Argonaute [41], the model assumes molecule availability as the only inhibition-limiting factor and describes the interaction between miRNA species a (ranging from 1 to M ) and RNA species i (ranging from 1 to N ) in terms of (see Fig 1a) • synthesis rates (b i for the RNA, β a for the miRNA) • degradation rates (d i for the RNA, δ a for the miRNA) • miRNA-RNA association and dissociation rates (k + ia and k − ia respectively) • miRNA-RNA complex degradation rates (κ ia for the catalytic pathway leading to miRNA recycling [42] and σ ia for the stoichiometric pathway without recycling) Assuming deterministic mass-action kinetics, molecular levels (m i for RNA species i, µ a for miRNA species a) evolve according to where c ia denotes the level of the complexes formed by RNA species i and miRNA species a. Such a system possesses a unique asymptotically stable steady state [43], where molecular levels attain the values (here and in the following, we denote the steady state value of variable x as [x]) , with m i ≡ b i /d i and µ a ≡ β a /δ a the concentrations of RNA and miRNA species in absence of inhibition. The quantities µ 0 ia and m 0 ia are given respectively by and effectively quantify the inverse repression strengths of miRNAs and RNAs. Specifically, see Eq (2), RNA species i is unrepressed (or, respectively, repressed) by miRNA species Hence the smaller is µ 0 ia the stronger the repression that a can exert on i. Similar considerations hold for m 0 ia : the smaller it is, the more miRNA species a will be sequestered by RNA species i.
The strength of miRNA-mediated RNA crosstalk at steady state can be estimated by the change in the steady-state level of RNA species i induced by a (small) variation in the transcription rate of species j, quantified by the susceptibility [18]  RNA species can crosstalk via chains of miRNA-mediated effective interactions, as can do species 1 and 3 in this example. (c1-c4) Classes of miRNA-RNA interactions considered in this work (following [45]): (c1) perfect k-base-pairing in the seed region ('k-mer' mode); (c2) seed base-pairing with up to one mismatched or bulged nucleotide (non-canonical mode or 'seed-nc'); (c3) non-seed base pairing with bulged and/or mismatched nucleotides ('noseed-9nt' mode); (c4) non-seed binding within weak diffuse regions ('noseed' mode). (d) Frequency of each binding mode in the CLASH dataset (from [45]). (e) Distributions of RNA transcription rates used in this work: each rate is assumed to be drawn independently from a lognormal distribution with given mean (same for each RNA species (The prefactor d j in Eq (5) serves the only purpose of making χ ij dimensionless.) Note that (i) χ ij ≥ 0 (i.e. the effective interaction tends to increase or decrease the levels of both RNAs) and (ii) χ ij and χ ji are a priori different (see [44] for a detailed discussion of this aspect). The advantage of the susceptibility over alternative measures of crosstalk, like the Pearson correlation coefficient, lies in the fact that it focuses on the role of competition, disregarding indirect effects due e.g. to fluctuations in miRNA levels. An extended comparison of different crosstalk measures can be found in [17]. Starting from Eq (2), one can derive an analytical expression allowing for the convenient computation of χ ij for each (i, j) pair in any miRNA-RNA network specified by a given set of kinetic parameters (see Supporting Text, Section 1). In compact form, the susceptibility matrix χ = {χ ij } N i,j=1 turns out to be given by where diag m m denotes the N -dimensional diagonal matrix with elements m i /m i (i = 1, . . . , N ) while W is an N × N matrix with elements the sum running over all miRNA species that co-target RNAs i and j. Note that Eq (6) implies that i and j need not be targeted by a common miRNA species in order for χ ij to be non-zero, as crosstalk can propagate through chains of miRNA-mediated interactions [17,27] (see Fig 1b). A toy model explicitly displaying this mechanism is discussed in Supporting Text, Section 2.

Choice of networks and parameters, and simulated scenarios
We shall mainly focus on the human miRNA interactome reconstructed in [45] using the CLASH (Crosslinking, Ligation And Sequencing of Hybrids) protocol. We refer to this as the 'CLASH interactome' for short; see Materials and Methods for details. The 4 types of miRNA-RNA couplings we consider are described in Fig 1c: (c1) perfect pairings of k miRNA seed nucleotides ("k-mer" for brevity, with k ranging from 6 to 9); (c2) sequence-specific pairings with up to one bulge or mismatch in the seed region (non-canonical pairings, or "seed-nc" for short); (c3) a 9 nt stems no-seed interaction allowing for bulged nucleotides in the target ("noseed-9nt"); and (c4) a no-seed interaction with distributed weak pairings ("noseed" for short). Non-canonical pairings are the most abundant in the CLASH interactome, accounting for roughly 77% of all miRNA-RNA interactions [45] (see Fig  1d). They are also weaker than canonical ones and seem to exert a very limited repressive role [46]. Nevertheless, they in principle contribute to miRNA titration and hence to RNA crosstalk. We therefore included them in our analysis. Our results will however turn out to be qualitatively independent of whether non-canonical sites are accounted for. A summary of basic features of the CLASH subnetworks spanned by different classes of interactions is given in Supporting Text, Supplementary Table 1.
For sakes of simplicity, we assume κ ia = κ and σ ia = σ for each (i, a) pair, δ a = δ for all a and d i = d for all i. With this choice, one has, for each miRNA-RNA pair, with λ = σ σ+κ the 'stoichiometricity ratio'. Using values of λ, d and δ compatible with empirical evidence (see Table I  miRNAs and b i for RNAs), and (ii) the values of either µ 0 ia or m 0 ia for each miRNA-RNA pair. We shall consider different scenarios for these quantities (see below). Once parameters are set, emergent crosstalk patterns are obtained by solving Eq (6) numerically.

Transcription rates and transcriptional heterogeneity (TH)
Throughout this study, we assume that both RNA transcription rates b i (i = 1, . . . , N ) and miRNA transcription rates β a (a = 1, ..., M) are log-normal i.i.d. random variables with means b and β, and variances σ 2 b and σ 2 β , respectively. The mean RNA transcription rate b is kept fixed at 8 molecules/h (see Table I), while we use the mean miRNA transcription rate β as a control parameter upon varying which crosstalk patterns are analyzed. To assess the impact of heterogeneity in transcription rates across molecular species (TH for short), we study how crosstalk patterns change as the magnitude of fluctuations increases, assuming the same transcriptional variability for miRNAs and RNAs. Our goal is to understand how the effective interaction network processes different degrees of variability in transcription rates, particularly at the level of RNAs. We hence tune TH by changing the coefficient of variation of individual rates (standard deviation over mean, see Fig 1e), which we denote by CV tr . In particular, we have exploited the log-normality of transcription rates to explore a 20-fold range of values of CV tr , from CV tr = 0.1 to 2.

Binding strengths and binding heterogeneity (BH)
To appraise how heterogeneities in the miRNA-RNA interaction strengths (binding heterogeneity or BH for short) affect the emergent crosstalk landscape, we consider three variants of the structure of binding affinities encoded in the CLASH interactome (see Fig 1f). At the lowest level of diversity, we assume a homogeneous network in which µ 0 ia = 2µ 0 for each miRNA-RNA pair, with µ 0 a constant taken to be equal to 4 molecules (see Table I). (Assuming a negligible miRNA-RNA unbinding rate, this corresponds roughly to an association rate of 0.02/molecules/hour, in agreement with values reported in [37].) At the intermediate level, we discriminate (i, a) pairs interacting via k-mer pairing (for which we take µ 0 ia = µ 0 , i.e. stronger coupling) from the rest (µ 0 ia = 2µ 0 ). Finally, at the highest level we associate different binding strengths to each of the four types of miRNA-RNA pairs, assuming a 2-fold change in µ 0 ia between groups in agreement with esti-  Table I. The yellow shaded area qualitatively marks the region where the mean susceptibility is significantly different from zero, which coincides with the susceptible regime [18]. In each case, the standard error of the mean is equal to or smaller than the size of the markers.

Role of network topology
We furthermore characterize the extent to which crosstalk patterns are induced by the specific wiring of the CLASH interactome by comparing it to the patterns that arise in randomized versions of the same network. In particular, we study ensembles of networks obtained by re-assigning each (i, a) link to a miRNA-RNA pair drawn randomly among all possible pairs with equal probability. This type of re-wiring disregards topological correlations of all orders, including node connectivities [47]. To evaluate the impact of the specific degree sequences encoded in the mapped miRNA interactome on the emergent system-level crosstalk patterns, we also analyze networks generated by a more conservative procedure based on degree-preserving edge-swaps [48]. Details of the latter are given in Supporting Text.
The mean RNA crosstalk intensity is a robust property of the miRNA-RNA network Expectedly, the overall abundance of free RNAs and free miRNAs change in opposite directions as miRNA transcription is globally upregulated and RNAs are increasingly repressed (see Fig 1g). Susceptibilities are bound to be larger when RNAs are more sensitive to changes in miRNA levels, i.e. in the so-called susceptible region at intermediate values of β [17]. Representative susceptibility distributions derived by solving Eq (6) and describing the CLASH network's crosstalk pattern for different degrees of TH and BH are displayed in miRNA availability modulates χ so that it peaks within the susceptible region and is vanishingly small outside of it (see Fig 2b), where molecular levels are practically unaffected by varying miRNA transcription rates. Notably, this picture is substantially unchanged by modifying the degrees of TH and/or BH, save for a modest expansion of the susceptible region. Such a behaviour therefore describes a 'basal level' of crosstalk that occurs in the network in any given condition. To appraise its significance, one can gauge it against the selfsusceptibility χ ii = d i ∂mi ∂bi , which quantifies the change in the level of free transcripts of species i induced by a small modification of its own transcription rate. (Note that, by definition, One sees that χ is about four orders of magnitude smaller than the mean self-susceptibility. In this respect, crosstalk appears to be on average very weak.
On the other hand, the picture just derived strongly suggests that the mean susceptibility profile is determined to a large extent by the topology of the network. We shall validate this hypothesis in the following. This conclusion, as well as the overall qualitative crosstalk characteristics illustrated by Fig  2a, will be seen to remain valid also when the contribution of non-canonical binding sites is disregarded.
The achievable crosstalk strength is enhanced by transcriptional heterogeneities Fig 2c displays the behaviour of the mean maximum susceptibility χ max = max (i,j) χ ij , where the maximum is taken over all pairs of different RNA species (i.e. with i = j). χ max quantifies the maximum achievable intensity of crosstalk interactions in each scenario, therefore providing a proxy for the strength of the most significant miRNA-mediated couplings arising between different RNA species in the network. Like χ , χ max also peaks in the susceptible regime, albeit for smaller values of the mean miRNA transcription rate. The strongest crosstalk hence typically occurs when RNA levels are just weakly sensitive to changes in miRNA availability. Remarkably, χ max is four orders of magnitude larger than χ . The backbone of the RNA crosstalk network formed by the most intense interactions is therefore comparable in intensity to the maximum achievable self-susceptibilities, see Supporting Text, Supplementary Fig S5. At odds with χ , however, χ max is strongly contextdependent, being modulated by both BH and (more significantly) TH. This finding agrees with the proposed role of kinetic heterogeneities in creating favourable paths in the miRNA-RNA network through which perturbations can efficiently propagate, as discussed e.g. in [17,27].
Crosstalk becomes more selective upon increasing heterogeneity Along with a higher potential for propagation, increased TH makes crosstalk more selective by systematically involving a smaller number of targets, both in terms of in-coming regulation (i.e. of the number of different transcripts that can regulate a given RNA) and, more significantly, in terms of outgoing regulation (i.e. of the number of different transcripts that are regulated by a given RNA). To quantify this aspect, we evaluated the quantities Both g i and h j vary between 0 and 1, as do S in and S out . A value g i 0 indicates that a large number of RNA species can almost equally affect the steady state of RNA species i, whereas a value of g i 1 indicates that RNA i is regulated by a small number of other RNA species. Likewise, when h j 0 a perturbation of the transcription rate of RNA j affects the steady state of a large number of other RNA species almost equally, while if h j 1 RNA j only affects a small number of other RNAs. In turn, g −1 i and h −1 i provide an indication of the number of upstream and, respectively, downstream miRNAmediated contacts of a given RNA species. If follows that S in (respectively S out ) represents the average of the inverse number of RNAs, a perturbation of which can considerably affect the level of a given RNA (resp. whose level can be affected by a perturbation of a given mRNA). We will call S in the incoming selectivity and S out the outgoing selectivity, respectively. Fig 3 displays the inverse incoming (Fig 3a) and outgoing (Fig 3b) selectivities as functions of the mean miRNA transcription rate β for the CLASH interactome. One sees that the typical number of crosstalk partners is modulated significantly only in the susceptible region. Higher degrees of transcriptional heterogeneity in particular tend to make crosstalk increasingly more selective (i.e. to lower S −1 in and S −1 out ). On the other hand, different degrees of binding heterogeneity (BH) appear to impact this scenario rather weakly.
Because the selectivity is ultimately a global property tied to susceptibility distributions, a few characteristics of these curves can be understood from features of the latter. For instance, the high selectivity achieved outside the susceptible region is likely due to the existence of strongly crosstalking pairs, enhanced by transcriptional heterogeneities (cf. Fig 2c). Peak inverse selectivity is instead achieved when susceptibility distributions tend to become more homogeneous (cf. Fig  S4). Likewise, transcriptional heterogeneity makes distributions less homogeneous, thereby increasing the selectivity. On the other hand, the divergent behaviour of incoming and outgoing components is harder to understand based on these aspects alone, as it possibly involves topological ingredients.
Keeping in mind that susceptibilities χ ij are not a priori symmetric, one can also quantify the degree of asymme- try in terms of the mean relative difference between χ ij and χ ji . Results (see Supporting Text, Section 3) show that, in the CLASH network, a properly defined asymmetry index is robustly maximized in the susceptible regime, where it can achieve a significantly high value that is weakly modulated by TH.
Summing up, while the behaviour of the mean crosstalk intensity appears to be hard-wired in the topology of the CLASH interactome, other features are tuned by the degree of heterogeneity. Most notably, crosstalk gets stronger and more selective as transcription rates become more diverse, while binding heterogeneities appear to specifically affect the maximum crosstalk intensity achievable. Finally, when crosstalk is strongest, individual crosstalk interactions tend to become more asymmetric, i.e. χ ij and χ ji are typically different. As this feature is observed independently of the degree of TH, the emergence of directional crosstalk appears to be an inherent property of miRNA-RNA networks.
Transcriptional heterogeneities elicit non-local RNA crosstalk Increased crosstalk intensity and selectivity are accompanied by the establishment of non-local effects, represented by strong effective interactions coupling RNAs that are separated by more than one miRNA species in the miRNA-RNA network. This phenomenon has been addressed e.g. in [27,44] in the context of small motifs. To quantify it in a large-scale network, we consider the quantity By definition, K ij is non-zero only if RNAs i and j are cotargeted by at least one miRNA species, while it vanishes for pairs (i, j) that do not share a miRNA regulator. (About 90% of potentially crosstalking RNA pairs involves species that are not co-regulated in the CLASH interactome.) In brief, as µ 0 ia is inversely proportional to the binding affinity between miRNA a and RNA i, larger values of K ij imply a stronger crosstalk potential between RNA species i and j based only on the network's local interaction structure and kinetic parameters. If crosstalk mostly occurs between co-regulated RNAs one should therefore expect the pattern of susceptibilities to match that of K ij s, at least qualitatively. We hence focus on the Pearson correlation coefficient between the K ij s and the susceptibilities χ ij s, i.e. of non-locality in crosstalk patterns (higher ρ implying more local crosstalk). Fig 4 shows the behaviour of ρ, the average being over realizations of TH. While the correlation peaks in the susceptible region, crosstalk patterns generically appear to correlate poorly with local topology in the CLASH interactome, as ρ 0.2. Most notably, ρ decreases significantly as TH is strengthened. This implicates kinetic heterogeneities in the establishment of extended interaction paths that reduce the effective diameter of the interactome by connecting distant RNAs via miRNA-mediated interactions. In this respect, miR-NAs appear to operate on RNAs both as specific repressors of individual transcripts and as a diffuse regulatory layer affecting the transcriptome as a whole.
The most marked effect induced by binding heterogeneities consists in an increase of ρ at small β and a shift of the peak correlation at smaller values of β. Interestingly, changes appear only when the full-fledged variability of binding sites is considered (high BH), while both the homogeneous case (low BH) and the case in which only k-mer and non-k-mer interactions are distinguished (medium BH) return very similar results. The particular structure of non-k-mer interactions reported in the CLASH data therefore only seems to bear a weak impact on the structure of crosstalk patterns. ia (or, equivalently, on the binding affinity k + ia ) to predict crosstalk interactions could be ineffective due to the significant long-range crosstalk that emerges as the network becomes more and more heterogeneous, especially in terms of transcription rates. This conclusion is most relevant in the susceptible regime, where cells presumably op-erate and RNA levels are more sensitive to changes in miRNA levels.
More robust expression profiles are associated to stronger RNA crosstalk By controlling the availability of their targets, miRNAs effectively process the variability induced by RNA transcription rates. In some cases (e.g. in presence of specific patterns of correlation between transcription rates), fluctuations can be reduced leading to more finely tuned expression levels [17,18,22,23]. In general, though, crosstalk tends to amplify target variability, especially when different species are transcribed independently [18]. The exact relationship between crosstalk intensity and transcript variability in extended networks is however bound to depend on the specific features of the crosstalk patterns.
In Fig 5a-b we show the coefficient of variation of RNA levels, averages being taken over many independent realizations of TH, as a function of the mean miRNA transcription rate β in different BH scenarios. Relative fluctuations exhibit a maximum at large values of β within the susceptible region and generically increase with the degree of TH. Variability in transcription rates therefore expectedly promotes variability in the resulting expression profiles. However the increase of fluctuations with respect to the unregulated case (β → 0) is negligible or very modest in a broad range of values of β within the susceptible region. On the other hand, at fixed TH, different BH scenarios do not appear to affect the robustness of expression profiles (see Fig 5b).
Recalling the behaviour of the maximal susceptibility χ max (see Fig 2b), one notices that the strongest maximal crosstalk is associated to more robust expression profiles within the susceptible region and, vice-versa, stronger fluctuations in expression profiles occur when crosstalk gets weaker (see Fig 5c). In other terms, uncorrelated transcriptional heterogeneities tend to be amplified when crosstalk is suppressed (higher miRNA expression levels), while they are more efficiently contained when the strongest crosstalk emerges. This scenario is summarized in Fig 5d: for any given degree of TH, as miRNA availability increases, crosstalk intensity on one hand and fluctuations of the output levels on the other are subject to a tradeoff that gets stronger as transcription rates becomes more homogeneous.
These results clearly implicate transcriptional heterogeneities as a key determinant of the stability of expression profiles even in presence of crosstalk, in line with previous observations on small networks [17,23]. It is however important to remark that this picture was obtained under the assumption of uncorrelated extrinsic fluctuations in RNA transcription rates. The presence of correlations might considerably alter this conclusion, as was first discussed in [18].

Crosstalk patterns are resilient to transcriptional heterogeneity
After analysing systemic properties, we now ask to what degree crosstalk patterns are preserved upon increasing the level of transcriptional heterogeneity. A global analysis shows (see Fig 6) that susceptibilities are remarkably well preserved statistically as the degree of transcriptional heterogeneity increases. Most notably, about 75% of the RNA pairs that are in the top sextile for crosstalk intensity at the lowest CV tr (CV tr = 0.1) persist in the top sextile when TH is 20-fold larger (CV tr = 2). Such a fraction is larger than would be expected by chance (about 58%), implying the existence of a significant backbone of RNA-RNA interactions resilient to transcriptional heterogeneity. A similar picture holds for the other sextiles. It is also instructive to inspect robustness specifically for RNA pairs that do not share any miRNA regulators, which amount to roughly 90% of the total. carry a susceptibility that is two orders of magnitude smaller than the maximum but two orders of magnitude larger than the average.) Weak sensitivity to changes in transcriptional heterogeneity would be expected if crosstalk interactions were functionally significant. Remarkably, this appears to be the case across a broad range of degrees of TH, both for short-range (mediated by a single miRNA species) and long-range (resulting from extended miRNA-mediated chains) crosstalk interactions.

Node degrees are the key topological determinants of the crosstalk scenario in the CLASH interactome
To appraise the role of the specific wiring encoded by the CLASH data in determining the scenario described so far, we compared our results against a null model obtained by randomly re-wiring the CLASH interactome. Specifically, we re-assigned each link to a randomly chosen miRNA-RNA pair, thereby preserving only the overall numbers of links and nodes while altering all other topological features like node degrees, degree-degree correlations, etc (see Methods). Each independent re-wiring process leads to a different final network. These randomized versions diverge from the original miRNA-RNA network in two basic aspects. In first place, they are slightly more compact, as evidenced by the distribution of the shortest miRNA-mediated paths between any two RNAs shown in Fig 7a. In addition, the randomization alters the distribution of node degrees by effectively eliminating the most highly connected RNA and miRNA species that are found in the CLASH data (see Fig 7b and 7c). Results obtained for key crosstalk descriptors in the CLASH and randomized networks (averaged over many realizations of the randomization protocol) are illustrated in Fig 7d- Randomized networks display a much larger (about twofold) mean susceptibility for crosstalk than the CLASH interactome, possibly due to the fact that miRNA targets are generically closer in the randomized versions. However, the maximum achievable crosstalk strength χ max is about 4 times smaller in the random networks compared to CLASH. Moreover, the susceptibility profile is more concentrated in the randomized network than it is for the CLASH network, reflecting a significantly narrower susceptible region (see Supporting Text, Fig S9). Naturally selected miRNA-RNA networks therefore appear to foster the emergence of stronger crosstalk links. In addition, the Pearson coefficient ρ quantifying the linear correlation between susceptibilities and local interaction parameters attains a much larger value in the randomized network with respect to the CLASH interactome throughout most of the susceptible region (see Fig 6f). miRNA-mediated crosstalk in random networks is therefore significantly more local, and thereby easily predictable e.g. by miRNA-RNA affinities, than it is in a network shaped by natural selection. Finally, expression profiles generated in the randomized network are slightly more stable than those found in the CLASH interactome (as quantified by the coefficient of variation, see Fig 6f). This feature is however more marked at higher miRNA expression levels, where RNA crosstalk is generically weaker. The basic traits of the RNA crosstalk emerging in randomized versions of the CLASH data are hence substantially different from those characterizing the interactome. Supporting Text, Section 3 and Fig S10 report the behaviour of the asymmetry and selectivity indices in randomized networks. At odds with the results obtained for the interactome (for which the asymmetry is weakly dependent on parameter heterogeneity), crosstalk in randomized networks becomes drastically more asymmetric and selective with increasing degrees of TH, although the number of interaction partners is generically higher in the randomized topology than it is in the interactome. In other terms, such features appear to be less robust to parameter heterogeneity in random structures than they are in naturally selected networks.
Note that, by applying a more conservative protocol that reshuffles miRNA-RNA links while preserving node degrees, one retrieves a crosstalk scenario that is essentially identical to that found for the original CLASH interactome (see Supporting Text, Section 4). This indicates that degree sequences (i.e. the topology of miRNA-RNA interactions encoded by the different types of couplings), as opposed to e.g. degree-degree correlations or other higher-order topological features, are the key geometric controllers of RNA crosstalk patterns. Enhanced crosstalk and non-locality therefore appear to be encoded by selection within the structure of the miRNA-RNA network interaction.

Canonical and non-canonical binding sites control different aspects of RNA crosstalk
A key question at this point is whether the observed crosstalk scenario is mainly due to the canonical (stronger) kmer pairings or, rather, if non-canonical (weaker) binding sites contribute to its establishment. A breakdown of the topology of the subnetworks induced by the different classes of interactions in the CLASH data shows significant similarities (see Supporting Text, Fig S11). Based on topology alone, then, appraising the role of non-canonical interactions is not simple.
To clarify this point, one can repeat the above analysis by successively adding each type of pairings shown in Fig 1c to the subnetwork induced by k-mer interactions in the CLASH data. After evaluating susceptibilities in each case, one sees (see Fig 8) that the crosstalk scenario underlied by the k-mer layer is qualitatively similar to that retrieved for the complete CLASH interactome. In particular, the k-mer network alone expectedly suffices to explain the maximum achievable crosstalk with quantitative accuracy. However, k-mer interac- The degree of TH is fixed at CVtr = 0.4 (with averages performed over 1000 independent realizations of TH) and high BH is assumed. In each case, the standard error of the mean is equal to or smaller than the size of the markers. The term 'CLASH-noseed' indicates the full CLASH interactome except for 'noseed' type of interactions. tions by themselves would yield stronger mean susceptibility, slightly more local crosstalk patterns and significantly larger variability of output profiles compared to the full network. Perhaps surprisingly, each of these aspects therefore appear to be quantitatively modulated to some degree by the weaker non-canonical interactions.
To further validate this picture, we have analyzed the RNA crosstalk scenario in the cancer-specific interactomes reconstructed in [37]. These networks comprise canonical pairings only and their basic topological characteristics are noticeably distinct from those found in CLASH (see Supporting Text, Fig S12). Results are summarized in Supporting Text, Fig  S13. The low-β behaviour starkly contrasts with that found in the CLASH reconstruction, in that crosstalk carries a stronger local component. In addition, maximal susceptibilities are roughly four times weaker in these networks. As a consequence, crosstalk is generically attenuated compared to CLASH and the potential to process (amplify) transcriptional heterogeneities is limited. Output profiles are consequently more stable against transcriptional variability across the whole range of levels of miRNA expression. These differences aside, the emergent crosstalk pattern robustly shows enhanced maximal intensity and non-locality with respect to their randomized counterpart, in qualitative agreement with the emergent crosstalk picture derived for the CLASH interactome.

Methodological choices
The system-level crosstalk scenarios studied here were derived under a few key methodological choices. First, we focused on the steady states of the mass action kinetics of miRNA-RNA interactions, Eq (2). While reasonable for timescales of the order of 1/d and 1/δ (and compatible with those analyzed in the experimental literature), this choice excludes from our analysis the rich phenomenology observed during transients [28], when crosstalk can be amplified over timescales determined by the details of the interaction kinetics. Likewise, we can't account directly for intrinsically dynamic regulatory mechanism like the dependence of miRNA decay rates on the round of recycling observed in [42]. Second, we opted to set a few parameters to values compatible with experimental evidence (see Table I) while treating miRNA and RNA transcription rates as independent, identically distributed quenched random variables with prescribed probability distributions. All our results were then obtained by averaging over many realizations of such parameters. Such an approach highlights features of the emerging crosstalk patterns that may be expected to be independent of the specific choice of transcription rates. On the flip side, we are un-able to characterize crosstalk for specific, possibly more realistic choices of transcription rates. Finally, we used the mean miRNA transcription rate as the only control parameter. While we explored a very broad range of values thereof, the physiologically relevant regime is likely to lie at intermediate miRNA transcription intensities, where RNA levels are more sensitive to changes in miRNA levels.
By assuming constant mean transcription rates we are effectively discarding the possibility that different RNAs or miRNAs are transcribed at very different rates (at least at low enough values of CV tr ). Significant inhomogeneities in the average biosynthesis rates of individual molecular species may affect our results. Yet, the highly interconnected structure of miRNA-RNA networks and the hierarchical organization of miRNA-RNA binding affinities [38] suggest that RNA crosstalk could be more influenced by global miRNA availability than by the specific structure of the miRNA population, at least in standard physiological conditions (e.e. in absence of strong miRNA induction). This is ultimately the scenario we probe in our study. Likewise, we are ignoring the possibility that heterogeneity parameters are correlated. As mentioned before, correlations between transcription rates would alter the picture regarding the processing of fluctuations [18]. Likewise, correlations between topological and transcriptional parameters, like those observed in [49], would the effects of heterogeneity, thereby significantly affecting crosstalk patterns.
To test the degree to which selection-shaped features of the miRNA interactome influence the emergent crosstalk pattern, we then studied how effective interactions are modulated by the structure of miRNA-RNA binding strengths and by the specific wiring encoded in data. In the former case, we were interested in evaluating the relevance for crosstalk of noncanonical binding sites, whose repressive efficiency is likely very limited [46]. In the latter, we aimed instead at understanding (i) whether crosstalk in naturally selected miRNA-RNA networks is qualitatively different from that arising in random networks and, if so, (ii) which topological features of real networks underlie the observed disparities.

Outlook
In broad terms, our analysis shows that RNA crosstalk in an extended network is modulated by miRNA availability both in terms of its basal level measured by the mean intensity and in terms of its maximal achievable strength. The typical crosstalk interaction is relatively weak. In specific, it is roughly four orders of magnitude smaller than the mean self-susceptibility, i.e. the mean change in the level of free transcript induced by a variation of its own transcription rate. Still, a multitude of strong crosstalk interactions arise, whose intensity is comparable to that of self-interactions. This in turn generates highly selective and directional crosstalk patterns. Notably, while co-regulated RNAs typically bear the strongest crosstalk links, non-co-regulated, or distant, RNAs can still crosstalk at significant intensities (roughly two orders of magnitude above the basal level). In such conditions, the typical RNA-to-RNA distance, in terms of number of links of the miRNA-mRNA interaction networks, above which one species can be considered to be insensitive to perturbations carried out on another species becomes comparable to the diameter of the network. A perturbation in the transcription level of one RNA can then be broadcast (via a chain of miRNA-mediated effective interactions) to distant RNA nodes, potentially propagating over the entire network. Such a feature is intrinsically due to competition and renders local kinetic parameters ineffective priors to predict crosstalk interactions. miRNAs therefore appear to manage a system-level regulatory layer where they operate collectively to mediate a complex, heterogeneous and robust network of RNA crossregulation. More work is however required to fully unravel its functional capabilities, especially concerning the buffering of fluctuations and gene expression noise.
The scenario we describe is qualitatively preserved if crosstalk is assumed to be carried by canonical interactions alone. In particular, the latter are highly effective modulators of crosstalk intensity. Non-canonical binding sites however, while substantially ineffective repression-wise, can enhance non-locality thereby extending the crosstalk range. Topological features of the naturally selected interactome were also found to bear a significant effect on crosstalk. Specifically, they lead to a broader susceptible region, higher maximal susceptibility, and more pronounced non-local effects than their randomized counterparts. In this respect, selection appears to have favoured the emergence of such features at system level.
It is important to stress that the crosstalk interactions on which we focus are quantified by susceptibilities, Eq (5). This implies that (i) they are driven by competition effects exclusively, and (ii) they are generated by small perturbations of RNA levels (as opposed e.g. to the models of [37,38]). Our scenario might therefore be close to a standard, homeostatic physiology in which transcription rates only undergo small variations. In this respect, the emergence of significant non-local effects is a surprising consequence of networking. Large perturbations, like the strong induction of a particular miRNA species, should be expected to amplify the crosstalk picture discussed here. However, responses to perturbations may become non-linear when the applied stimulus exceeds a threshold [28]. In such conditions, susceptibilities or standard correlation coefficients are likely inappropriate to describe crosstalk. More theoretical work on miRNA-RNA networks is required to fully sort out this case.
Unfortunately, probing the crosstalk scenario we describe in experiments could be challenging essentially due to weakness and non-locality. To validate the picture we describe, both in terms of individual interactions and of global features, one may however resort to transcriptomic data. Recent work has indeed identified a specific group of correlation functions that, under certain conditions, yield excellent approximations for the real susceptibilities [44]. Evaluating such quantities on RNA readouts would then provide a direct, data-driven snapshot of RNA crosstalk.

Conclusion
Besides their important role as negative controllers of gene expressiom, miRNAs mediate the establishment of extended networks of RNA cross-regulation. Several features of these networks appear to be hard-wired in the topology of the underlying miRNA-RNA interactome, while others are modulated by transcriptional and/or binding heterogeneities. Whereas the typical crosstalk interaction generated by small changes in RNA availability is weak, non-local effects are significant. Crosstalk-based regulation therefore appears essentially as a system-wide phenomenon, enhanced by variability in kinetic parameters. In physiological conditions, such a regulatory layer can potentially contribute to a variety of functions, such as the processing of transcriptional heterogeneities and the coordination of large-scale rearrangements of RNA levels, similar to the responses observed in [38]. The broader picture we have derived might however apply more generally to networks of molecular species competing for a common resource.

MATERIALS AND METHODS miRNA interactomes
For the CLASH interactome, after parsing the original bipartite network derived in [45] to remove degeneracies and disjoint nodes, we found N = 6, 943 RNA species (implying about 4.8 × 10 7 potential crosstalk interactions) and M = 383 miRNA species connected by 17,411 edges carrying different binding strengths. The same pipeline was applied to the tumor-type specific miRNA-RNA networks obtained in [37] and based on the Cupid protocol for predicting microRNA-target interactions [50], which accounts for canon-ical pairings exclusively. The resulting miRNA-RNA networks are considerably larger than the CLASH interactome, as evidenced by the comparison of degree distributions given in Supporting Text, Fig S12. Computational analysis With parameters set as described, Eq (6) was solved numerically for each of the networks cosidered using Python scripts based on NumPy [51] and SciPy [52]. The code is available from https://github.com/matmi8/RNAnet. In presence of TH, results were averaged over multiple independent realizations of the vectors b = {b i } N i=1 and β = {β a } M i=1 of RNA and miRNA transcription rates (respectively) for each value of CV tr . The number of realizations was chosen in each case to ensure a stable estimation of different quantities. Details are given in figure captions. All other parameters, both kinetic and topologic, were kept fixed. Likewise, in the case of topological heterogeneity, results were averaged over 100 networks obtained by independently randomizing the original miRNA-RNA network while keeping all other parameters, both transcriptional and kinetic, fixed. 100 independent randomizations of the interactome sufficed to ensure stable averages in each condition.

Derivation of RNA susceptibilities in generic miRNA-RNA networks
To derive Eq (6) of the main text, we start recasting the expressions for [m i ] and [µ a ] (see Eq. (2) of main text) as Using these, we immediately obtain where we used the identities Eq (20) can be re-cast in the compact form where χ is the susceptibility matrix (with elements χ ij ), diag m m denotes the diagonal matrix with elements It follows that Recalling that, if all eigenvalues of W are strictly smaller than 1 in absolute values (as is easily verified numerically to the case in this study), one has 1 − Z −1 = n≥0 Z n , one finds that Expressions (20) and (21) clarify an important point. While W ij is different from zero only if RNAs i and j are coregulated by at least one miRNA species, the elements of W n are different from zero if there is at least one chain of n miR-NAs joining RNAs i and j. In practice, this is what allows for crosstalk to occur even between RNAs that are not directly co-regulated, as shown explicitly within a toy model in the following section.

Susceptibility between distant RNA pairs
To show explicitly how non-zero susceptibilities can arise between pairs of RNAs connected by chains of miRNAmediated couplings from Eq 20, we compute here the susceptibility matrix for a toy network formed by three RNA and two miRNA species, Fig S1. The W matrix for this network reads Two elements (w 13 and w 31 ) are nil since RNAs 1 and 3 are not co-targeted by any miRNA. Nevertheless, using (20), the susceptibility matrix turns out to be given by where w ij = w ij − 1. Hence a non-zero susceptibility binds RNA species 1 and 3, which are connected by the chain of interactions passing through RNA 2. This connection is also evidenced by the form of the corresponding elements of χ.
The above equation also shows explicitly that, in general, χ ij and χ ji are different.
FIG. S1. Toy miRNA-RNA network. Two miRNA species and three RNA species interact by direct couplings represented by the continuous blue lines. Dotted lines denote instead the effective crosstalk interactions that are established between RNAs as a consequence of competition to bind miRNAs. These correspond in turn to the non-zero elements of the susceptibility matrix χ.

Crosstalk asymmetry
To measure crosstalk directionality, we define the quantity such that 0 ≤ ∆ ij ≤ 1. In short, the closer ∆ ij is to zero (resp. one) the closer crosstalk between RNAs i and j is to being symmetric (resp. fully asymmetric). A global measure of asymmetry is conveniently obtained by computing the average asymmetry over all pairs of different RNAs in the network, i.e.
Results for this quantity are reported in Fig S2 for both the CLASH network and its randomized variant. Crosstalk asymmetry is generically larger in the susceptible regime, more pronouncedly so in the CLASH network than in its randomized version. Notably, the asymmetry profile is roughly independent of the degree of binding heterogeneity while it is only weakly modulated by transcriptional variability in the CLASH network. As seen for the mean crosstalk intensity, this state of things suggests that the way in which crosstalk asymmetry is tuned by the mean miRNA transcription rate β is an inherent property of miRNA-RNA networks, that is mainly encoded in their topology. The striking difference that can be seen between the behaviour of ∆ in real (CLASH) and random networks (see Fig S2b) supports this intuition.
FIG. S2. Crosstalk asymmetry in the CLASH network and its randomized counterparts. Profile of ∆ obtained (a) in the CLASH network, and (b) in its randomized version in the 3 scenarios considered for binding heterogeneity and for various degrees of transcriptional heterogeneity. Curves are averaged over 100 independent realizations of transcription rate profiles. Results for the random case are additionally averaged over 100 independent realizations of the randomization process. In each case, the standard error of the mean is equal to or smaller than the size of the markers.  Table 1 (Main Text). The yellow shaded area qualitatively marks the region where the mean susceptibility is significantly different from zero, which coincides with the susceptible regime. In each case, the standard error of the mean is equal to or smaller than the size of the markers. The self-susceptibility is maximal when miRNA levels are low, in which case the availability of free RNA molecules increases roughly linearly with the transcription rate. As β increases, miRNA repression gets stronger and selfsusceptibilities decrease until, at large enough β, RNAs are fully repressed and therefore insensitive to small changes in their transcription rates. (c) Comparison between maximum self-susceptibility (averaged over TH realizations), mean self-susceptibility (averaged over TH realizations) and χmax for different degrees of TH in the high BH scenario. The intensity of crosstalk between different RNAs, measured by the latter quantity, is indeed of the same order of magnitude as self-susceptibilities.  1c).