Translating ceRNA susceptibilities into correlation functions

TL;DR

This paper demonstrates that ceRNA susceptibilities to perturbations can be inferred from correlation functions, providing a practical method to analyze ceRNA crosstalk in biological datasets.

Contribution

It establishes a direct link between ceRNA susceptibilities and correlation functions, enabling easier estimation of ceRNA interactions from data.

Findings

Correlation functions encode ceRNA susceptibilities

Validated approach with TCGA breast cancer data

Numerical simulations support the theoretical framework

Abstract

Competition to bind microRNAs induces an effective positive crosstalk between their targets, therefore known as `competing endogenous RNAs' or ceRNAs. While such an effect is known to play a significant role in specific conditions, estimating its strength from data and, experimentally, in physiological conditions appears to be far from simple. Here we show that the susceptibility of ceRNAs to different types of perturbations affecting their competitors (and hence their tendency to crosstalk) can be encoded in quantities as intuitive and as simple to measure as correlation functions. We confirm this scenario by extensive numerical simulations and validate it by re-analyzing PTEN's crosstalk pattern from TCGA breast cancer dataset. These results clarify the links between different quantities used to estimate the intensity of ceRNA crosstalk and provide new keys to analyze transcriptional datasets and effectively probe ceRNA networks in silico.

Figures

Tables

Table 1. Table 1: Parameters used in Fig 2 . Note that i ∈ { 1 , 2 } 𝑖 1 2 i\in\{1,2\} while a = 1 𝑎 1 a=1 .

Parameters [units]	Figs 2A	Figs 2B	Figs 2C
$b_2$, ${\beta }_1$ [molec. min-1]	10	10	10
$k_{11}^{+}$ [molec.-1 min-1]	$e^{-5}$	$e^{-5}$	$e^{-5}$
$k_{21}^{+}$ [molec.-1 min-1]	$e^{-6}$	$e^{-6}$	$e^{-6}$
$k_{11}^{-}$, $k_{21}^{-}$ [min-1]	0.001	0.001	0.001
$d_1$, $d_2$, ${\delta }_1$ [min-1]	0.1	0.05	0.2
${\sigma }_{11}$, ${\sigma }_{21}$ [min-1]	1	0.05	0.05
${\kappa }_{11}$, ${\kappa }_{21}$ [min-1]	0.001	0.001	0.001

Equations

$${\rho }_{ij}=\frac{⟨m_i{m}_j⟩-⟨m_i⟩⟨m_j⟩}{[(⟨m_i^2⟩-{⟨m_i⟩}^2)(⟨m_j^2⟩-{⟨m_j⟩}^2){]}^{1/2}},$$

$${\chi }_{ij}=\frac{\partial ⟨m_i⟩}{\partial b_j}\ge 0,$$

$${\dot{\mu }}_a={\beta }_a-{\delta }_a{\mu }_a-\sum _i{k}_{ia}^{+}{m}_i{\mu }_a+\sum _i(k_{ia}^{-}+{\kappa }_{ia})c_{ia}$$

$${\dot{m}}_i=b_i-d_i{m}_i-\sum _a{k}_{ia}^{+}{m}_i{\mu }_a+\sum _a{k}_{ia}^{-}{c}_{ia}$$

$${\dot{c}}_{ia}=k_{ia}^{+}{m}_i{\mu }_a-\frac{c_{ia}}{{\tau }_{ia}},$$

$$⟨c_{ia}⟩=k_{ia}^{+}{\tau }_{ia}{m}_i{\mu }_a.$$

$${\dot{m}}_i$$

$${\dot{\mu }}_a$$

$$\begin{gathered} L=-\sum _i(b_i\log m_i-d_i{m}_i) \\ -\sum _a({\beta }_a\log {\mu }_a-{\delta }_a{\mu }_a)+\sum _{i,a}{k}_{ia}^{+}{m}_i{\mu }_a. \end{gathered}$$

$$⟨f⟩=\frac{1}{Z(T)}\sum _{\mathbf{m},\mathbf{\mu }}f(\mathbf{m},\mathbf{\mu })e^{-L/T},$$

$$Z(T)=\sum _{\mathbf{m},\mathbf{\mu }}{e}^{-L/T}$$

$${⟨fg⟩}_c\equiv ⟨fg⟩-⟨f⟩⟨g⟩,$$

$$⟨m_i⟩=-T\frac{\partial }{\partial d_i}\log Z(T),$$

$${⟨m_i{m}_j⟩}_c=-T\frac{\partial ⟨m_i⟩}{\partial d_j}\equiv -T{\omega }_{ij},$$

$${⟨m_i\log m_j⟩}_c=T{\chi }_{ij}.$$

$$X_{ij}={⟨m_i\log m_j⟩}_c$$

$$C_{ij}={⟨m_i{m}_j⟩}_c$$

$${\omega }_{ij}\simeq -\frac{⟨{\mu }_1⟩}{{\mu }_{0,i1}{\mu }_{0,j1}}\frac{b_i{b}_j}{d_i{d}_j}{\left({\delta }_1+\sum _{\ell }{a}_{\ell }\frac{b_{\ell }{k}_{\ell 1}^{+}}{d_{\ell }}\right)}^{-1},$$

$$C_{\text{PTEN,ceRNA}}={⟨m_{\text{PTEN}}{m}_{\text{ceRNA}}⟩}_c$$

$$X_{\text{ceRNA,PTEN}}={⟨m_{\text{ceRNA}}\log m_{\text{PTEN}}⟩}_c$$

$$X_{\text{PTEN,ceRNA}}={⟨m_{\text{PTEN}}\log m_{\text{ceRNA}}⟩}_c$$

$${\dot{\mu }}_a={\beta }_a-{\delta }_a{\mu }_a-\sum _i{k}_{ia}^{+}{m}_i{\mu }_a+\sum _i(k_{ia}^{-}+{\kappa }_{ia})c_{ia},$$

$${\dot{m}}_i=b_i-d_i{m}_i-\sum _a{k}_{ia}^{+}{m}_i{\mu }_a+\sum _a{k}_{ia}^{-}{c}_{ia},$$

$${\dot{c}}_{ia}=k_{ia}^{+}{m}_i{\mu }_a-\frac{c_{ia}}{{\tau }_{ia}},$$

$${\dot{\mu }}_a\simeq {\beta }_a-{\delta }_a{\mu }_a-\sum _i\frac{{\sigma }_{ia}}{{\sigma }_{ia}+{\kappa }_{ia}+k_{ia}^{-}}{k}_{ia}^{+}{m}_i{\mu }_a,$$

$${\dot{m}}_i\simeq b_i-d_i{m}_i-\sum _a\frac{{\sigma }_{ia}+{\kappa }_{ia}}{{\sigma }_{ia}+{\kappa }_{ia}+k_{ia}^{-}}{k}_{ia}^{+}{m}_i{\mu }_a.$$

$${\dot{\mu }}_a\simeq {\beta }_a-{\delta }_a{\mu }_a-\sum _i{k}_{ia}^{+}{m}_i{\mu }_a,$$

$${\dot{m}}_i\simeq b_i-d_i{m}_i-\sum _a{k}_{ia}^{+}{m}_i{\mu }_a,$$

$$L=-\sum _i(b_i\log m_i-d_i{m}_i)-\sum _a({\beta }_a\log {\mu }_a-{\delta }_a{\mu }_a)+\sum _{i,a}{k}_{ia}^{+}{m}_i{\mu }_a.$$

$$\frac{{\dot{m}}_i}{m_i}$$

$$\frac{{\dot{\mu }}_a}{{\mu }_a}$$

$$\dot{L}=\sum _i\frac{\partial L}{\partial m_i}\dot{m_i}+\sum _a\frac{\partial L}{\partial {\mu }_a}\dot{{\mu }_a}=-\sum _i\frac{{\dot{m_i}}^2}{m_i}-\sum _a\frac{{\dot{{\mu }_a}}^2}{{\mu }_a}\le 0.$$

$$m_i=\frac{b_i}{d_i}{F}_i(\mu ),F_i(\mu )=\frac{{\mu }_{0,i}}{{\mu }_{0,i}+\mu },{\mu }_{0,i}=\frac{d_i}{k_i^{+}}.$$

$$F_i(\mu )\simeq \{\begin{array}{ll}1-\frac{\mu }{{\mu }_{0,i}}& \text{for }\mu \ll {\mu }_{0,i}\text{(ceRNA i expressed)}\\ \frac{1}{2}-\frac{\mu -{\mu }_{0,i}}{4{\mu }_{0,i}}& \text{for }\mu \simeq {\mu }_{0,i}\text{(ceRNA i susceptible)}\\ \frac{{\mu }_{0,i}}{\mu }& \text{for }\mu \gg {\mu }_{0,i}\text{(ceRNA i repressed)}\end{array}$$

$$\frac{d{F}_i}{d\mu }\equiv F_i^{\prime }\simeq \{\begin{array}{ll}-\frac{1}{{\mu }_{0,i}}& \text{for }\mu \ll {\mu }_{0,i}\text{(ceRNA i expressed)}\\ -\frac{1}{4{\mu }_{0,i}}& \text{for }\mu \simeq {\mu }_{0,i}\text{(ceRNA i susceptible)}\\ -\frac{{\mu }_{0,i}}{{\mu }^2}& \text{for }\mu \gg {\mu }_{0,i}\text{(ceRNA i repressed)}\end{array}$$

$$\mu \simeq \frac{\beta -{\sum }_{i\in \text{Repr}}{b}_i-\frac{1}{4}{\sum }_{i\in \text{Susc}}{b}_i}{\delta +{\sum }_{i\in \text{Expr}}\frac{b_i{k}_i^{+}}{d_i}+\frac{1}{4}{\sum }_{i\in \text{Susc}}\frac{b_i{k}_i^{+}}{d_i}},$$

$$\frac{\partial \mu }{\partial d_j}\simeq {\chi }_{\mu \mu }\frac{b_j}{d_j}\frac{\mu }{{\mu }_{0,j}}\times \{\begin{array}{ll}1& \text{for }\mu \gg {\mu }_{0,i}\text{(ceRNA i expressed)}\\ \frac{1}{4}& \text{for }\mu \simeq {\mu }_{0,i}\text{(ceRNA i susceptible)}\\ 0& \text{for }\mu \gg {\mu }_{0,i}\text{(ceRNA i repressed)}\end{array}$$

Full text

Introduction

MicroRNAs (miRNAs) are small non coding RNA (ncRNA) molecules that post-transcriptionally regulate a significant portion of the eukaryotic transcriptome via sequence-specific, protein-mediated binding in the cytoplasm [1]. Their primary effects on coding transcripts consist in inhibiting translation and fostering degradation [2]. Long ncRNAs, instead, can transiently sequester miRNAs, thereby altering their availability and overall repressive potential [3]. Following early observations concerning small regulatory RNAs in plants and bacteria [4, 5], competition to bind miRNAs has been hypothesized to cause an effective positive interaction (‘crosstalk’) between their coding and/or non-coding targets that may directly affect protein levels [6] (see Fig. 1A,B). Several experimental and modeling studies have clarified the conditions under which such a scenario may become biologically relevant, highlighting specifically how molecular levels and kinetic heterogeneities may control it [1, 8, 9, 10, 11, 12, 13]. So far, such a ‘ceRNA effect’ (whereby ceRNA stands for ‘competing endogenous RNA’) has been quantitatively validated in cases of differentiation [14], disease [15] or in presence of unphysiologically large transcriptional inputs [16]. Its significance in standard physiological conditions is therefore subject to scrutiny [17].

A major difficulty in detecting the ceRNA effect unambiguously in experiments or data lies in the fact that it should be disentangled from other mechanisms that may bear a similar impact, i.e. an effective positive coupling, on transcripts. Imagine a network of $N$ ceRNA species interacting with $M$ miRNA species. ceRNA levels $m_{i}$ ($i=1,\ldots,N$) fluctuate stochastically in time due to random synthesis and degradation events and to interactions with miRNAs, whose levels are also subject to random fluctuations. Denoting by $\left\langle\cdot\right\rangle$ the time-average in the steady state, an effective ceRNA-ceRNA dependence can be signaled by a statistical correlation coefficient such as Pearson’s [8], i.e.

[TABLE]

with the idea that, if $\rho_{ij}$ is large enough, a perturbation altering the level of ceRNA $j$ will cause part of the miRNA population to move from one target to the other, effectively broadcasting the perturbation from ceRNA $j$ to ceRNA $i$ through miRNA-mediated interactions.

A more direct description of this mechanism is attained instead via susceptibilities like [1]

[TABLE]

where $b_{j}$ stands for the transcription rate of ceRNA $j$. $\chi_{ij}$ quantifies the shift in the mean level of ceRNA $i$ caused by a (small) variation in $b_{j}$, and a large $\chi_{ij}$ (assuming no direct control of ceRNA $i$ by ceRNA $j$) points to miRNA-mediated crosstalk between ceRNAs $i$ and $j$ (see Fig. 1C).

While both $\chi_{ij}$ and $\rho_{ij}$ capture aspects of ceRNA crosstalk seen in experiments, their underlying physical meaning is a priori different. Fluctuating miRNA levels naturally correlate co-regulated targets, so that a large $\rho_{ij}$ is obtained when both ceRNAs respond to the stochastic dynamics of their regulator. This however does not necessarily imply a large $\chi_{ij}$. In fact, $\chi_{ij}$ can be large even in absence of fluctuations in miRNA levels, i.e. as a consequence of competition alone. In such conditions, $\rho_{ij}$ vanishes. $\chi_{ij}$ has indeed been found to be asymmetric under exchange of its indices (i.e. $\chi_{ij}\neq\chi_{ji}$ in general) [1], at odds with $\rho_{ij}$ which is necessarily symmetric. It would therefore be important to clarify how quantities like (1) and (2) are related in miRNA-ceRNA networks, especially to understand whether responses to perturbations (a central quantity of interest for many potential applications of the ceRNA effect) can be encoded in quantities as intuitive and as simple to measure experimentally or from data as a Pearson correlation coefficient.

Here we show that the information conveyed by $\chi_{ij}$ is indeed captured by a correlation function similar to $\rho_{ij}$. On the other hand, $\rho_{ij}$ is linked to a susceptibility, i.e. to the response of a target to a perturbation altering the level of its competitor, but the perturbation concerns the intrinsic decay rate of the competitor rather than its transcription (as is the case for $\chi_{ij}$). In the following, we will derive these results and validate them by computer simulations and gene expression data analysis, and explore their consequences.

Results

Theory

We start from the dynamics of molecular populations in a miRNA-ceRNA network, denoting by $m_{i}$ the level of ceRNA species $i$ (ranging from 1 to $N$), by $\mu_{a}$ the level of miRNA species $a$ (ranging from 1 to $M$), and by $c_{ia}$ the levels of miRNA-ceRNA complexes. In the deterministic limit where stochastic fluctuations are neglected, the time evolution of concentration variables is described by

[TABLE]

with the different parameters denoting intrinsic synthesis ($b_{i}$, $\beta_{a}$) and degradation rates ($d_{i}$, $\delta_{a}$), complex association/dissociation rates ($k_{ia}^{\pm}$) and complex processing rates ($\sigma_{ia}$ and $\kappa_{ia}$ for stoichiometric and catalytic processing, respectively), while $\tau_{ia}=(\sigma_{ia}+\kappa_{ia}+k_{ia}^{-})^{-1}$ represents the mean lifetime of the complex formed by miRNA $a$ and ceRNA $i$. We note that if the mean lifetime of complexes is much shorter than that of free molecular species, i.e. if $\tau_{ia}\ll 1/d_{i}$ and $\tau_{ia}\ll 1/\delta_{a}$ for each $i$ and $a$, miRNA-ceRNA complexes achieve a steady state much faster than miRNA and ceRNA levels. In such conditions, $\dot{c}_{ia}\simeq 0$ and one can eliminate complexes from (3) by replacing $c_{ia}$ with its steady state value

[TABLE]

For $\kappa_{ia}+k_{ia}^{-}\ll\sigma_{ia}$ (i.e. when stoichiometric degradation without miRNA recycling is the dominant channel of complex processing), this allows to re-cast (3) in the form (see Supporting Text)

[TABLE]

where $L$ is a function of all miRNA levels $\bm{\mu}=\{\mu_{a}\}$ and all ceRNA levels $\mathbf{m}=\{m_{i}\}$ given by

[TABLE]

One easily sees (see Supporting Text) that $L$ decreases along trajectories of (23), implying that its minimum describes the physically relevant steady state of (3) with $\mathbf{m\neq 0}$ and $\bm{\mu\neq 0}$.

If intrinsic molecular noise (arising from stochastic transcription and degradation events and from titration due to miRNA-ceRNA interactions) is added to (3), after a transient molecular levels will eventually stabilize and fluctuate over time around the steady state described by the minimum of $L$. We are interested in finding a compact and intuitive mathematical form for the correlations arising between the different components in such conditions. Molecular noise is Poissonian, namely the strength of fluctuations affecting each variable is proportional to the square root of mean molecular levels (see e.g. [12] for an explicit representation in the context of a miRNA-ceRNA network), which makes our goal especially challenging. However we will see that the effects of molecular noise can be remarkably well approximated by a uniform “effective temperature” $T$ representing the strength of fluctuations affecting all molecular species involved. In this case, one can describe fluctuations around the steady state as thermal fluctuations around a Boltzmann-Gibbs equilibrium state. This allows to compute averages of generic functions of $\mathbf{m}$ and $\bm{\mu}$ as “thermal averages”, i.e.

[TABLE]

where

[TABLE]

is a normalization factor, the deterministic limit being obtained for $T\to 0$. In particular, defining

[TABLE]

by straightforward calculations one finds

[TABLE]

Therefore, in this approximation, the susceptibility $\chi_{ij}$ [Eq. (2)] is linked to the correlation function [Eq. (12)]

[TABLE]

which, as $\chi_{ij}$, is not symmetric under the exchange of $i$ and $j$, while the ceRNA-ceRNA covariance

[TABLE]

is tied to the susceptibility $\omega_{ij}$ quantifying the change in $\left\langle m_{i}\right\rangle$ induced by a (small) change of the intrinsic degradation rate $d_{j}$ of ceRNA $j$ [Eq. (11)]. (Note that $\omega_{ij}\leq 0$.) The constant linking these quantities is the temperature $T$ quantifying the strength of the uniform “effective noise”.

Somewhat unexpectedly, the above results suggest that the susceptibility $\omega_{ij}$ must be symmetric under exchange of $i$ and $j$, i.e., for instance, if the level of ceRNA $i$ is altered by changing the intrinsic degradation rate of ceRNA $j$, then the reverse is also true. To check this property, one can calculate $\omega_{ij}$ explicitly for a system formed by $N$ ceRNA species interacting with a single miRNA species at steady state by following a different route, specifically along the lines of [1]. Considering the repression strength to which ceRNAs $i$ and $j$ are subject at a given (mean) level $\left\langle\mu_{1}\right\rangle$ of miRNA species 1 ($M=1$ in this case), one finds, for each ceRNA, a soft “threshold” value of $\left\langle\mu_{1}\right\rangle$, denoted by $\mu_{0,i1}\simeq d_{i}/k_{i1}^{+}$, such that $i$ is unrepressed (resp. repressed or susceptible to changes in $\left\langle\mu_{1}\right\rangle$) if $\left\langle\mu_{1}\right\rangle\ll\mu_{0,i1}$ (resp. $\gg\mu_{0,i1}$ or $\simeq\mu_{0,i1}$). A direct calculation (see Supporting Text) shows that $\omega_{ij}$ can attain large values only if both ceRNAs are susceptible to $\left\langle\mu_{1}\right\rangle$, in which case one has ($i\neq j$)

[TABLE]

where $a_{\ell}=0,1,1/4$ if ceRNA $\ell$ is repressed, unrepressed or susceptible, respectively. Eq. (15) confirms that $\omega_{ij}$ is indeed symmetric under exchange of $i$ and $j$.

Concerning the approximations under which the the above results were obtained, we remark that we started by considering (3) in the limit of (i) fast complex equilibration, and (ii) miRNA-ceRNA complex processing dominated by the stoichiometric channel, with the former playing the key role in deriving the function $L$ (see Supporting Text). We note however that the overall scenario just described also holds for when complexes evolve over time scales much longer than those of free molecular levels, i.e. for $\tau_{ia}\gg 1/d_{i}$ and $\tau_{ia}\gg 1/\delta_{a}$. In particular, (3) can again be re-cast in the form of (23) with $L$ given by (22), albeit with re-scaled transcription rates (see Supporting Text for details).

Therefore we conclude that, as long as molecular noise can be approximated by a uniform effective temperature,

1. (i)

the ceRNA-ceRNA covariance $C_{ij}=\left\langle m_{i}m_{j}\right\rangle_{c}$ is a proxy for the susceptibility $\omega_{ij}$, and 2. (ii)

the correlation function $X_{ij}=\left\langle m_{i}\log m_{j}\right\rangle_{c}$ is a proxy for the susceptibility $\chi_{ij}$.

Validation

We have validated the above scenario by simulating a small network involving 2 ceRNA and a single miRNA species via the Gillespie algorithm [2], where molecular noise is accounted for explicitly (see Supporting Text). Results are summarized in Fig. 2, where we compare $\omega_{12}$, $\omega_{21}$ and $C_{12}\equiv C_{21}$ on on hand, and $\chi_{12}$, $\chi_{21}$, $X_{12}$ and $X_{21}$ on the other, as computed from simulations (i.e. with the actual molecular noise), against the theoretical predictions. We considered three scenarios for the mean lifetime of miRNA-ceRNA complexes, namely those covered by the theory (i.e. complex equilibration much faster and much slower than the equilibration of miRNA and ceRNA levels) as well as the intermediate case where characteristic timescales are comparable for all variables.

One sees that theoretical predictions obtained in the “thermal noise” approximation agree remarkably well with simulations including the actual molecular noise. In particular, the full correspondence between the susceptibilities $\omega_{ij}$ and $\chi_{ij}$ and the (re-scaled) correlation functions $C_{ij}$ and, respectively, $X_{ij}$ is evident. Notice that a single global parameter $T\geq 0$ has been used to fit all data in each of the conditions. This shows how accurately the assumption of a uniform effective temperature can mimic the effects of intrinsic stochasticity. On the other hand, its limits might be reflected, at least in part, in the discrepancies that occur at high transcription rates.

These results confirm that (12) and (14) are indeed good predictors of the response of a ceRNA to a perturbation affecting one of its competitors within a miRNA-ceRNA network. Notably, such correlation functions are easy to estimate from transcription data sets. Our framework therefore has the potential to offer new insight into post-transcriptional regulation, its system-level organization and its impact on cellular functions.

In order to test this idea, we analyzed the ceRNA scenario emerging from 1098 breast cancer samples obtained from The Cancer Genome Atlas [19], focusing on the widely studied oncosuppressor PTEN and its immediate competitors (i.e. the ceRNAs sharing at least one miRNA regulator with PTEN). In particular, we computed

[TABLE]

for a set of candidate PTEN ceRNAs found in [15] by means of Mutually Targeted miRNA-Response Element Enrichment Analysis. Notice that the average appearing in Eqs. (16–18) is over samples and not over time. We expect however that, if the interaction network is conserved across samples, averages over samples should reproduce statistical averages such as (13), as different samples effectively represent different snapshots of the state of the network. Fig. 3A shows that when (16) (whose value is encoded in the color of markers) is large, both (17) and (18) tend to be large. According to (15), a large $C_{ij}$ (or $\omega_{ij}$) signals that both PTEN and its competitor are susceptible to changes in the level of at least one of their shared regulators.

For such pairs, in addition, it has previously been shown that both $\chi_{ij}$ and $\chi_{ji}$ are expected to be large [1]. This implies a fully bi-directional crosstalk, i.e. any perturbation affecting the level of one species should affect the level of the other via miRNA-mediated regulation. Remarkably, this was experimentally shown to be the case in [15] for some of the ceRNAs we tested (e.g. SERINC1, VAPA), all of which are in this regime according to our analysis. Adding to this, we are also able to point to a number of other PTEN competitors, a perturbation of which should trigger a response by PTEN.

On the other hand, smaller values of (16) (orange markers in Fig. 3A) are associated to strongly asymmetric PTEN-ceRNA pairs for which (18) is much larger than (17). This suggests that PTEN will respond to an increase of its competitor’s bare transcription rate (and not vice-versa), while no response of PTEN should be expected upon perturbing the bare decay rate of the same ceRNA as $C_{ij}$ is small. Within the steady state theory of [1], ceRNA pairs with strongly different values of $\chi_{ij}$ and $\chi_{ji}$ pertain to cases where the responding ceRNA (PTEN here) is susceptible to variations in the miRNA levels while the perturbed one (PTEN’s competitor) is fully repressed. Our data analysis fully confirms both this scenario and the theory presented here in linking such cases to low values of the bare covariance (14).

Finally, note (Fig. 3B) that the above information can not be retrieved if $C_{ij}$ is replaced by the Pearson coefficient $\rho_{ij}$, Eq. (1), which just amounts to normalizing $C_{ij}$ by the product of the standard deviations $\sigma_{m_{i}}$ and $\sigma_{m_{j}}$ of $m_{i}$ and $m_{j}$. Indeed, using the value of $\rho_{ij}$ to color-code PTEN’s ceRNAs, one sees that the Pearson coefficient can mislead into expecting (or not expecting) a response to a perturbation when the actual susceptibilities are small (resp. large).

For instance, $\rho_{ij}$ is rather small for the pair formed by PTEN and SLC1A2, which seems to suggest absence of mutual cross-talk between these two transcripts. However, while both $C_{\text{PTEN,SLC1A2}}$ and $X_{\text{SLC1A2,PTEN}}$ are small, $X_{\text{PTEN,SLC1A2}}$ is significant. This suggests that (i) SLC1A2 will not respond to a perturbation affecting the transcription rate of PTEN, and (ii) the pair will be insensitive to changes in each other’s bare decay rate; however, (iii) PTEN will be affected by a change in the bare transcription rate of its competitor despite the small statistical correlation that exists between their levels. Likewise, the large value of the Pearson coefficient between PTEN and DTWD2 can mislead into generically expecting a response when instead the susceptibility is strongly perturbation-dependent. In particular, the level of DTWD2 should not be significantly modified by a change in the level of PTEN (as $X_{{\rm DTWD2,PTEN}}$ is rather small) in spite of the large Pearson coefficient. Notice that, remarkably, for this pair, $C_{ij}$ and $\rho_{ij}$ take on very different values.

Discussion

To sum up, we have identified [Eq.s (11) and (12)] a set of correlation functions that can serve as proxies for ceRNA susceptibilities to perturbations. Specifically, $C_{ij}=\left\langle m_{i}m_{j}\right\rangle_{c}$ is related to the susceptibility $\omega_{ij}$ quantifying ceRNA $i$’s response to a change of the bare decay rate of ceRNA $j$, while $X_{ij}=\left\langle m_{i}\log m_{j}\right\rangle_{c}$ is related to the susceptibility $\chi_{ij}$ quantifying ceRNA $i$’s response to a change of the bare transcription rate of ceRNA $j$. These relations are valid at steady state and within the approximations discussed, are fully confirmed by numerical simulations.

Most importantly, quantities like $C_{ij}$ and $X_{ij}$ can be easily estimated from data and possibly measured in experiments. An analysis of PTEN’s emergent crosstalk pattern from TCGA breast cancer dataset using these functions has indeed shown that a map of ceRNA responses to perturbations affecting competitors can be constructed by combining the information provided by each, while the Pearson coefficient $\rho_{ij}$ can be inaccurate in this respect. This opens the way to probing the structure and function of ceRNA networks in silico by straightforwardly analyzing transcriptional data, and provides a key to obtain testable transcriptome-scale predictions about ceRNA crosstalk.

Notice that our results apply without any modification to ceRNA pairs that don’t share miRNA regulators, i.e. it is capable of identifying long-range crosstalk (i.e. interactions between ceRNAs that are separated by multiple miRNAs along the miRNA-ceRNA network) of the kind discussed in [20].

From the viewpoint of physics, results like (11) and (12) are akin to the “fluctuation-response relations” that constitute a cornerstone of statistical mechanics [21]. Their derivation in our context has relied on an equilibrium framework that presupposes stationarity of molecular levels. Since ceRNA crosstalk can be substantially more complex away from the steady state [22], a more refined mathematical study will be required to extend the theory developed here to off-equilibrium dynamical regimes. Our results on the other hand may also open the way to the application of recently developed inference techniques [23] to estimate miRNA levels or kinetic parameters from ceRNA levels.

Acknowledgments

We gratefully acknowledge Carla Bosia and Andrea Pagnani for useful insight and suggestions.

I Derivation of the function $L$ for fast complex processing

Starting from Eq. (3) of the Main Text, namely (with $\tau_{ia}=(\sigma_{ia}+\kappa_{ia}+k_{ia}^{-})^{-1}$)

[TABLE]

we assume that complexes equilibrate much faster than miRNA and ceRNA levels and substitute $c_{ia}$ with its steady state value $\left\langle c_{ia}\right\rangle=k_{ia}^{+}\tau_{ia}m_{i}\mu_{a}$. One finds

[TABLE]

For $\kappa_{ia}+k_{ia}^{-}\ll\sigma_{ia}$, this reduces to

[TABLE]

which is easily seen to be equivalent to Eq. (5) in the Main Text, with

[TABLE]

II $L$ decreases along the dynamics

By direct differentiation and using the fact that (see Eq. (5) in the Main Text)

[TABLE]

one finds

[TABLE]

In other words, under the approximations discussed above, $L$ decreases along the dynamics of the miRNA-ceRNA network. Therefore the minimum of $L$ (which is unique by virtue of the concavity of $L$) describes a steady state of the dynamics (19).

III Approximate calculation of $\omega_{ij}$ for a system with one miRNA and $N$ ceRNA species

Starting from Eq. (21) taken for $N$ ceRNA species and a single miRNA species (we suppress its index for simplicity) in the limit $\sigma_{i}\gg k_{i}^{-}+\kappa_{i}$, the steady-state level of $m_{i}$ reads

[TABLE]

Now noting that

[TABLE]

and that the steady state miRNA level can be approximated by [1]

[TABLE]

so that

[TABLE]

we can compute the susceptibility $\omega_{ij}$ as

[TABLE]

One finds

[TABLE]

where $\widehat{W}$ is a $3\times 3$ matrix that only depends on the regime $R(i)$ (repressed, susceptible or expressed) to which ceRNA $i$ belongs. By considering the definitions of the different regimes in terms of the value of $\mu$, all elements of $\widehat{W}$ are found to be $\ll 1$ (for instance, $W_{{\rm Expr,Expr}}=\mu/(\mu_{0,i}\mu_{0,j})\ll 1$ as $\mu\ll\mu_{0,i}$ and $\mu\ll\mu_{0,j}$ if ceRNAs $i$ and $j$ are both expressed) except for $W_{{\rm Susc,Susc}}$, which is given by

[TABLE]

leading immediately to Eq. (15) of the Main Text.

IV Case of slow complex processing

Assuming complex levels $c_{ia}$ are roughly stationary over time scales for which $m_{i}$ and $\mu_{a}$ evolve (i.e. $\tau_{ia}\gg 1/d_{i}$ and $\tau_{ia}\gg 1/\delta_{a}$ for each $i$ and $a$), then all terms in (19) that involve the variables $c_{ia}$ can be taken to be roughly constant for short enough characteristic times. In such conditions, miRNAs are effectively transcribed at rates

[TABLE]

while ceRNAs are effectively transcribed at rates

[TABLE]

In this limit, (19) can again be cast as

[TABLE]

with

[TABLE]

The main difference from the previous case lies in the fact that the minimum of $L$ should now be computed self-consistently from the asymptotic value of $c_{ia}$: after the (fast) equilibration of $m_{i}$’s and $\mu_{a}$’s following (36), a new steady state value for complexes is computed as $c_{ia}=k_{ia}^{+}\tau_{ia}m_{i}\mu_{a}$, leading in turn to new values for the effective transcription rates $\beta_{a}^{{\rm eff}}$ and $b_{i}^{{\rm eff}}$ and hence to new values for $m_{i}$’s and $\mu_{a}$’s from (36), and so on until convergence.

V Stochastic dynamics of a miRNA-ceRNA network

The time evolution of our miRNA-ceRNA network with $N$ ceRNA species (labeled $i$), $M$ miRNA species (labeled $a$) and intrinsic (molecular) noise is described by the system

[TABLE]

where $\tau_{ia}=(\sigma_{ia}+\kappa_{ia}+k_{ia}^{-})^{-1}$ while $\eta_{a}$, $\xi_{a}$ and $\zeta_{ia}$ represent stochastic variables. As each noise source contributes independently to the overall noise level, one has

[TABLE]

where $\xi_{m_{i}}$, $\eta_{\mu_{a}}$, $\zeta_{ia}^{\pm}$, $\zeta_{ia}^{\sigma}$, and $\zeta_{ia}^{\kappa}$ and are mutually independent zero-average random variables representing, respectively, the intrinsic noise in ceRNA levels, in miRNA levels, in the binding/unbinding dynamics of complexes, in the stoichiometric complex degradation channel and in the catalytic complex degradation channel. Correlations are, for each component, described by

[TABLE]

where

[TABLE]

denote the mean steady-state molecular levels. To obtain Fig. 2 of the Main Text, we have simulated the above system with $M=1$, $N=2$ using the Gillespie algorithm [2].