Modelling Competing Endogenous RNA Networks

Abstract

MicroRNAs (miRNAs) are small RNA molecules, about 22 nucleotide long, which post-transcriptionally regulate their target messenger RNAs (mRNAs). They accomplish key roles in gene regulatory networks, ranging from signaling pathways to tissue morphogenesis, and their aberrant behavior is often associated with the development of various diseases. Recently it has been experimentally shown that the way miRNAs interact with their targets can be described in terms of a titration mechanism. From a theoretical point of view titration mechanisms are characterized by threshold effect at near-equimolarity of the different chemical species, hypersensitivity of the system around the threshold, and cross-talk among targets. The latter characteristic has been lately identified as competing endogenous RNA (ceRNA) effect to mark those indirect interactions among targets of a common pool of miRNAs they are in competition for. Here we propose a stochastic model to analyze the equilibrium and out-of-equilibrium properties of a network of [Formula: see text] miRNAs interacting with [Formula: see text] mRNA targets. In particular we are able to describe in detail the peculiar equilibrium and non-equilibrium phenomena that the system displays in proximity to the threshold: (i) maximal cross-talk and correlation between targets, (ii) robustness of ceRNA effect with respect to the model's parameters and in particular to the catalyticity of the miRNA-mRNA interaction, and (iii) anomalous response-time to external perturbations.

Figure 1. Representation of a generic miRNA-target interaction network.

(A) Simplified picture of a miRNA-ceRNA interaction network. (B) For each miRNA () and ceRNA () present in the network we consider the main steps of transcription (rates and respectively) and degradation (rates and respectively) plus a titrative interaction between miRNA and ceRNA. miRNA and ceRNA can therefore form a complex with effective association rate . The parameter (the catalyticity parameter) tells which is the probability a miRNA is recycled after having interacted with one of its targets.

Figure 2. Threshold, noise and Pearson's coefficients varying ceRNA transcription rate.

(A–C) Steady state value for means, Fano factors and coefficients of variation for each free molecular species in a system of two miRNAs (miRNA1 and miRNA2, green and orange lines respectively) interacting with two ceRNAs (ceRNA1 and ceRNA2, blue and red lines respectively) varying the concentration of ceRNA1. In proximity to the threshold the system shows hypersensitivity to changes in the control parameter (ceRNA1 transcription rate), captured by a maximum in the Fano factors (panel B). For the same values of ceRNA1 transcription rate, the local maximum in the coefficients of variation (panel C) is the fingerprint of bimodal distributions in the number of molecules for each molecular species. (D) Pearson's coefficients between the two miRNAs (orange line) and the two ceRNAs (blue line). The two lines show a maximum in proximity to the ceRNA1 transcriptiom rate threshold value, meaning that there is a region of parameters where the fluctuations in the number of ceRNAs or miRNAs are highly correlated. Lines are the results of Gaussian approximation while symbols are Gillespie's simulations. For panels B,C the line color-code is the same as in panel A.

Figure 3. Threshold, noise and Pearson's coefficients varying miRNA transcription rate.

(A–C) Steady state value for means, Fano factors and coefficients of variation for each free molecular species in a system of two miRNAs (miRNA1 and miRNA2, green and orange lines respectively) interacting with two ceRNAs (ceRNA1 and ceRNA2, blue and red lines respectively) varying the concentration of miRNA1. In proximity to the threshold the system shows hypersensitivity to changes in the control parameter (miRNA1 transcription rate), captured by a maximum in the Fano factors (panel B). For the same values of miRNA1 transcription rate, the local maximum in the coefficients of variation (panel C) is the fingerprint of bimodal distributions in the number of molecules for each molecular species. (D) Pearson's coefficients between the two miRNAs (orange line) and the two ceRNAs (blue line). The two lines show a maximum in proximity to the miRNA1 transcriptiom rate threshold value, meaning that there is a region of parameters where the fluctuations in the number of ceRNAs or miRNAs are highly correlated. Lines are the results of Gaussian approximation while symbols are Gillespie's simulations. For panels B,C the line color-code is the same as in panel A.

Figure 4. Selectivity of miRNA and ceRNA interactions.

(A) Example of a network of ten miRNAs interacting with ten ceRNAs in blocks. The interaction links are such that we can define two main blocks (block 1 and block2) of strongly interacting miRNAs-ceRNAs connected by one common miRNAs (miRNA 5 in block 1, miRNA 6 in block 2) and ceRNAs (ceRNA 5 in block 1 and ceRNA 6 in block 2). Panels (B,C) show an example of dynamics of such network. Varying ceRNA1 (panel B) or miRNA10 (panel C) transcription rate during time (every 35 hours in the example, but the time is arbitrary) has a differentiated effect on the other ceRNAs and miRNAs present in the all network. The color-code for lines in panels B and C follows the color of miRNAs and ceRNAs in panel A.

Figure 5. Threshold effect in a miRNA-target catalytic interaction.

Example of a system of one miRNA interacting with two ceRNAs with cataliticity parameter . The threshold effect is possible only if the system is out of equilibrium (A). Numerical integration of Equation (1) in File S1 leads to time evolution of each molecular species for a given set of parameters. In panels A,C we plot "pictures" of the evolving system at different time (panel A , panel C ) as a function of ceRNA1 transcription rate. When t is smaller than the time complexes need to reach their steady state a threshold effect is observed. In panels B,D we plot the corresponding Pearson's coefficient profiles. Where the threshold effect is present (panel A), a peak in the Pearson's coefficient is also observed.

Figure 6. Response times upon one ceRNA perturbation.

Increasing miRNA transcription rate ceRNA1 shows a maximum and a minimum in its response times upon switching on or off ceRNA2 transcription respectively (panel A and B). The maximum (minimum) is located near the threshold, where ceRNA1 initial value (that is its values before switching on (off) ceRNA2) is near to the steady state it will reach upon switching on (off) ceRNA2 (panel C) but is also more sensitive to ceRNA2 variation (look at the maximum in the Pearson's correlation coefficient between ceRNA1 and ceRNA2 in panel D). Different color lines correspond to different numbers of ceRNAs in competition for the same miRNA. The qualitative trend for response times and Pearson's correlation coefficient is robust with respect to increasing such number.