Independent channels for miRNA biosynthesis ensure efficient static and dynamic control in the regulation of the early stages of myogenesis

TL;DR

This study models the miRNA regulatory network in early myogenesis, revealing how independent biosynthesis channels and competition control molecular levels and enable rapid differentiation responses.

Contribution

It provides a quantitative analysis of two independent miRNA biosynthesis pathways and their combined role in controlling early myogenesis, highlighting their necessity for efficient regulation.

Findings

Optimal static control constrains kinetic parameters.

Alternative transcription ensures rapid concentration shifts.

Both channels are essential for circuit functionality.

Abstract

Motivated by recent experimental work, we define and study a deterministic model of the complex miRNA-based regulatory circuit that putatively controls the early stage of myogenesis in human. We aim in particular at a quantitative understanding of (i) the roles played by the separate and independent miRNA biosynthesis channels (one involving a miRNA-decoy system regulated by an exogenous controller, the other given by transcription from a distinct genomic locus) that appear to be crucial for the differentiation program, and of (ii) how competition to bind miRNAs can efficiently control molecular levels in such an interconnected architecture. We show that optimal static control via the miRNA-decoy system constrains kinetic parameters in narrow ranges where the channels are tightly cross-linked. On the other hand, the alternative locus for miRNA transcription can ensure that the fast concentration shifts required by the differentiation program are achieved, specifically via non-linear response of the target to even modest surges in the miRNA transcription rate. While static, competition-mediated regulation can be achieved by the miRNA-decoy system alone, both channels are essential for the circuit's overall functionality, suggesting that that this type of joint control may represent a minimal optimal architecture in different contexts.

Tables

Table 1. Table 1: Processes included in the model and their associated parameters. The transcription rates b h subscript 𝑏 ℎ b_{h} (controller) and b μ subscript 𝑏 𝜇 b_{\mu} (miRNA) are the key control parameters used in this study. Note that α 𝛼 \alpha and k μ ​ ℓ subscript 𝑘 𝜇 ℓ k_{\mu\ell} are functions of [ h ] delimited-[] ℎ [h] , see Eq. ( 1 ), that vary respectively in [ 0 , 1 ] 0 1 [0,1] and [ 0 , k μ ​ ℓ max ] 0 superscript subscript 𝑘 𝜇 ℓ max [0,k_{\mu\ell}^{\rm max}] . The value of k μ ​ ℓ max superscript subscript 𝑘 𝜇 ℓ max k_{\mu\ell}^{\rm max} is as indicated. Numerical values used to obtain Fig. 2 are shown separately, since a specific fold-size perturbation Δ Δ \Delta , describing the increase in μ 𝜇 \mu biosynthesis due to the alternative channel (see ( 10 )), was chosen in order to reproduce empirical time courses.

Processes and rates		Parameters [units]	In all Figs.	In Fig. 2
Ex-novo synthesis and degradation	$\varnothing \underset{d_m}{\overset{b_m}{\rightleftharpoons }}m$	$b_m$ [$\text{nM}\cdot {\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
		$d_m=d$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
	$\varnothing \underset{d_q}{\overset{b_q}{\rightleftharpoons }}q$	$b_q$ [$\text{nM}\cdot {\text{s}}^{-1}$]	${10}^{-4}$	$6\cdot {10}^{-5}$
		$d_q=d$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
	$\varnothing \underset{d_h}{\overset{b_h}{\rightleftharpoons }}h$	$b_h$ [$\text{nM}\cdot {\text{s}}^{-1}$]	variable	$8\cdot {10}^{-5}$
		$d_h=d$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
	$\varnothing \underset{d_{\mu }}{\overset{b_{\mu }}{\rightleftharpoons }}\mu$	$b_{\mu }$ [$\text{nM}\cdot {\text{s}}^{-1}$]	variable	$4\cdot {10}^{-6}$
		$d_{\mu }=d$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
Synthesis from the precursor	$q\stackrel{\left(1-\alpha \right)b}{⟶}\mu$	$\alpha$ [adim.]
	$q\stackrel{\alpha b}{⟶}\ell$	$b$ [${\text{s}}^{-1}$]	${10}^{-4}$	$8\cdot {10}^{-5}$
Degradation of the sponge ($\ell$)	$\ell \stackrel{d_{\ell }}{⟶}\varnothing$	$d_{\ell }=d$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
miRNA-ceRNA complex formation	$\mu +m\stackrel{k_{\mu m}}{\to }{c}_m$	$k_{\mu m}$ [${\text{nM}}^{-1}\cdot {\text{s}}^{-1}$]	${10}^{-3}$	$4\cdot {10}^{-4}$
	$\mu +\ell \stackrel{k_{\mu \ell }}{\to }{c}_{\ell }$	$k_{\mu \ell }^{\max }$ [${\text{nM}}^{-1}\cdot {\text{s}}^{-1}$]	${10}^{-3}$	$1.2\cdot {10}^{-4}$
	$\mu +h\stackrel{k_{\mu h}}{\to }{c}_h$	$k_{\mu h}$ [${\text{nM}}^{-1}\cdot {\text{s}}^{-1}$]	${10}^{-3}$	$4\cdot {10}^{-4}$
Stoichiometric complex decay (without recycling of $\mu$)	$c_m\stackrel{{\sigma }_m}{⟶}\varnothing$	${\sigma }_m=\sigma$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
	$c_{\ell }\stackrel{{\sigma }_{\ell }}{⟶}\varnothing$	${\sigma }_{\ell }=\sigma$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
	$c_h\stackrel{{\sigma }_h}{⟶}\varnothing$	${\sigma }_h=\sigma$ [${\text{s}}^{-1}$]	${10}^{-4}$	$4\cdot {10}^{-5}$
Catalytic complex decay (with recycling of $\mu$)	$c_m\stackrel{{\kappa }_m}{⟶}\mu$	${\kappa }_m=\kappa$ [${\text{s}}^{-1}$]	$5\cdot {10}^{-4}$	$2\cdot {10}^{-4}$
	$c_{\ell }\stackrel{{\kappa }_{\ell }}{⟶}\mu$	${\kappa }_{\ell }=\kappa$ [${\text{s}}^{-1}$]	$5\cdot {10}^{-4}$	$2\cdot {10}^{-4}$
	$c_h\stackrel{{\kappa }_h}{⟶}\mu$	${\kappa }_h=\kappa$ [${\text{s}}^{-1}$]	$5\cdot {10}^{-4}$	$2\cdot {10}^{-4}$
Hill indices		$n$ [adim.]	2	2
		$p$ [adim.]	2	2
Dissociation constants		$h_{\alpha }$ [nM]	$1.0$	$0.5$
		$h_{\mu \ell }$ [nM]	$1.0$	$0.7$
Perturbation fold-size		$\Delta$ [adim.]		5

Equations

$$\alpha =\frac{[h{]}^n}{[h{]}^n+h_{\alpha }^n}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{k}_{\mu \ell }=k_{\mu \ell }^{\max }\frac{[h{]}^p}{[h{]}^p+h_{\mu \ell }^p}\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[m]=b_m-(d+k_{\mu m}[\mu ])[m]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[\mu ]=b_{\mu }+(1-\alpha )b[q]+\kappa ([c_{\ell }]+[c_m]+[c_h])-(d_{\mu }+k_{\mu m}[m]+k_{\mu \ell }[\ell ]+k_{\mu h}[h])[\mu ]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[\ell ]=\alpha b[q]-(d+k_{\mu \ell }[\mu ])[\ell ]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[h]=b_h-(d+k_{\mu h}[\mu ])[h]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[q]=b_q-(d+b)[q]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[c_m]=-(\sigma +\kappa )[c_m]+k_{\mu m}[\mu ][m]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[c_{\ell }]=-(\sigma +\kappa )[c_{\ell }]+k_{\mu \ell }[\mu ][\ell ]\text{leavevmode}\text{leavevmode},$$

$$\frac{d}{dt}[c_h]=-(\sigma +\kappa )[c_h]+k_{\mu h}[\mu ][h]\text{leavevmode}\text{leavevmode}.$$

$$b_{\mu }(t)=b_{\mu }\left[1+\Delta \theta (t-t^{\star })\right]\text{leavevmode}\text{leavevmode},$$

$$d\simeq \frac{1}{\tau }\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}\sigma \simeq \frac{1}{\tau }\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}{b}_j\simeq \frac{\gamma }{\tau }\text{leavevmode}\text{leavevmode},$$

$$k_{\mu i}\in \left[\frac{1}{100\gamma \tau },\frac{100}{\gamma \tau }\right]\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}\kappa \in \left[\frac{1}{10\tau },\frac{10}{\tau }\right]\text{leavevmode}\text{leavevmode},$$

$$[\mu ]=\frac{b_{\mu }+(1-\alpha )b[q]+\kappa ([c_{\ell }]+[c_m]+[c_h])}{d_{\mu }+k_{\mu m}[m]+k_{\mu \ell }[\ell ]+k_{\mu h}[h]}\text{leavevmode}\text{leavevmode},$$

$$[m]=\frac{b_m}{d+k_{\mu m}[\mu ]}\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[\ell ]=\frac{\alpha b[q]}{d+k_{\mu \ell }[\mu ]}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[h]=\frac{b_h}{d+k_{\mu h}[\mu ]}\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[q]=\frac{b_q}{d+b}\text{leavevmode}\text{leavevmode},$$

$$[c_m]=\frac{k_{\mu m}[\mu ][m]}{\sigma +\kappa }\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[c_{\ell }]=\frac{k_{\mu \ell }[\mu ][\ell ]}{\sigma +\kappa }\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[c_h]=\frac{k_{\mu h}[\mu ][h]}{\sigma +\kappa }\text{leavevmode}\text{leavevmode}.$$

$${\chi }_{mh}=\frac{\partial [m]}{\partial [h]}$$

$${\chi }_{m{b}_{\mu }}=\frac{\partial [m]}{\partial b_{\mu }}\text{leavevmode}\text{leavevmode},$$

$$[\mu ]=\frac{b_{\mu }+(1-\alpha )b[q]}{d_{\mu }+w(k_{\mu m}[m]+k_{\mu \ell }[\ell ]+k_{\mu h}[h])}\text{leavevmode}\text{leavevmode},$$

$$w=\frac{\sigma }{\sigma +\kappa }$$

$$[m]=\frac{b_m}{d+k_{\mu m}[\mu ]}\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[\ell ]=\frac{\alpha b[q]}{d+k_{\mu \ell }[\mu ]}\text{leavevmode}\text{leavevmode}\text{leavevmode},\text{leavevmode}\text{leavevmode}\text{leavevmode}[h]=\frac{b_h}{d+k_{\mu h}[\mu ]}\text{leavevmode}\text{leavevmode}.$$

$${\chi }_{m{b}_{\mu }}=\frac{[\mu ][m{]}^2{k}_{\mu m}}{b_m}G(\mathbf{\chi })\text{leavevmode}\text{leavevmode},$$

$${\chi }_{\ell b_{\mu }}=\frac{[\mu ][\ell {]}^2{k}_{\mu \ell }}{\alpha b[q]}G(\mathbf{\chi })+\frac{[\ell ]}{\alpha }\left(\frac{d\alpha }{d[h]}-\frac{[\ell ][\mu ]}{b[q]}\frac{d{k}_{\mu \ell }}{d[h]}\right){\chi }_{h{b}_{\mu }}\text{leavevmode}\text{leavevmode},$$

$${\chi }_{h{b}_{\mu }}=\frac{[\mu ][h{]}^2{k}_{\mu h}}{b_h}G(\mathbf{\chi })\text{leavevmode}\text{leavevmode},$$

$$G(\mathbf{\chi })=\frac{w[\mu ]\left[k_{\mu m}{\chi }_{m{b}_{\mu }}+k_{\mu \ell }{\chi }_{\ell b_{\mu }}+\left(k_{\mu h}+\frac{b[q]}{[\mu ]w}\frac{d\alpha }{d[h]}+[\ell ]\frac{d{k}_{\mu \ell }}{d[h]}\right){\chi }_{h{b}_{\mu }}\right]-1}{b_{\mu }+(1-\alpha )b[q]}\text{leavevmode}\text{leavevmode},$$

$$\frac{d\alpha }{d[h]}=\frac{n{h}_{\alpha }^n[h{]}^{n-1}}{([h{]}^n+h_{\alpha }^n{)}^2}\text{leavevmode}\text{leavevmode},$$

$$\frac{d{k}_{\mu \ell }}{d[h]}=k_{\mu \ell }^{\max }\frac{p{h}_{\mu \ell }^p[h{]}^{p-1}}{([h{]}^p+h_{\mu \ell }^p{)}^2}\text{leavevmode}\text{leavevmode}.$$

$$\displaystyle\chi_{mb_{\mu}}=\frac{k_{\mu m}b_{h}[\mu][m]^{2}}{w[\mu]^{2}\mathcal{A}-b_{m}b_{h}\Big{[}b_{\mu}+(1-\alpha)b[q]\Big{]}}\leavevmode\nobreak\ \leavevmode\nobreak\ ,$$

$$\mathcal{A}=k_{\mu \ell }^2[\ell {]}^2\frac{b_m{b}_h}{\alpha b[q]}+k_{\mu m}^2[m{]}^2{b}_h+k_{\mu h}[h{]}^2{b}_m\left[k_{\mu h}+\left(\frac{k_{\mu \ell }[\ell ]}{\alpha }+\frac{b[q]}{[\mu ]w}\right)\frac{d\alpha }{d[h]}+\left(1-\frac{[\ell ][\mu ]k_{\mu \ell }}{\alpha b[q]}\right)[\ell ]\frac{d{k}_{\mu \ell }}{d[h]}\right]\text{leavevmode}\text{leavevmode}.$$

$$\begin{gathered} {\chi }_{mh}=\frac{[m{]}^2{k}_{\mu m}[\mu ]\alpha b[q]}{\alpha b[q]b_m[b_{\mu }+(1-\alpha )b[q]]-w[\mu {]}^2({k_{\mu m}}^2[m{]}^2\alpha b[q]+{k_{\mu \ell }}^2[\ell {]}^2{b}_m)}\times \\ \times \left\{b[q]\frac{d\alpha }{d[h]}+w[\mu ]\left[k_{\mu h}+\frac{d{k}_{\mu \ell }}{d[h]}\ell +\frac{k_{\mu \ell }[\ell ]}{\alpha }\left(\frac{d\alpha }{d[h]}-\frac{[\ell ][\mu ]}{b[q]}\frac{d{k}_{\mu \ell }}{d[h]}\right)\right]\right\}. \end{gathered}$$

$$\beta \equiv \frac{b_{\mu }}{(1-\alpha )b[q]}=\frac{b_{\mu }([h{]}^n+h_{\alpha }^n)}{b[q]h_{\alpha }^n}\text{leavevmode}\text{leavevmode},$$

$$\mathrm{IR}(\Delta)=\int_{0}^{\infty}\Big{[}[m](t+t^{\star})-[m](t^{\star})\Big{]}dt\leavevmode\nobreak\ \leavevmode\nobreak\ ,$$

Full text

1 Introduction

Following years of experimental work, microRNAs (miRNAs) –small, endogenous, non-coding RNA molecules ubiquitously found in plant and animal cells– have emerged as key agents in post-transcriptional regulation [1]. Their primary mode of action is by protein-mediated base-pairing to target RNA molecules. For coding RNAs, this leads to the repression of gene expression through mRNA cleavage or translational inhibition [2, 3, 4, 5]. Their targets however include both coding and non-coding transcripts like long non-coding RNAs (lncRNAs) [6, 7, 8, 9], which places them at the center of the RNA-based regulatory web. Over time, miRNAs have been found to constitute a remarkably versatile regulatory layer, mediating functions ranging from the buffering of gene expression noise [10] to the timing of genetic circuits [11]. Moreover, they are now known to be profoundly implicated in a variety of developmental and disease processes [12]. Mathematical models have been employed to address the role of miRNAs in different contexts [13], highlighting for instance how miRNA-based regulation may combine with circuit topology [14], kinetic heterogeneities [15, 16] and effects due to competition between miRNA targets [17] to generate diverse functional outcomes. Competition, in particular, has been hypothesized to affect regulation in a broader, yet more subtle, way through the so-called ‘competing endogenous RNA’ (ceRNA) effect [18], whereby co-regulated targets can establish, under specific conditions [19, 20, 21, 22, 17], an effective crosstalk with potentially far-reaching implications. While experimental validations are currently putting under scrutiny the question of how effective this mechanism can be in standard conditions [23, 24, 25], competition has been shown to be central in a number of situations, perhaps most notably in cancer development [26, 27, 28] and muscle cell differentiation [29, 30].

Following recent work that has shed new light on its intricate genetic circuitry [31], we focus here on miRNA-based control in early myogenesis. Key actors include two miRNA species (miR-133 and miR-135), two transcription factors (MAML1 and MEF2C), a skeletal muscle-specific lncRNA (linc-MD1) and the RNA-binding human antigen R (HuR) protein (see Fig. 1A). miR-133 can be produced from a precursor RNA (pre-miR-133b) as well as from two independent genomic loci. However, pre-miR-133b also provides the substrate for the synthesis of linc-MD1 through a pathway alternative (and mutually exclusive) to that leading to miR-133. In addition, linc-MD1 possesses two target sites for miR-135 and one for miR-133 and can therefore act as a ‘decoy’ for both miRNAs. The transcription factors MAML1 and MEF2C, both essential in the expression of muscle-specific genes [32], are instead targets of miR-133 and miR-135, respectively. As a consequence, linc-MD1 is a ceRNA of MAML1 (resp. MEF2C) and competes with it to bind miR-133 (resp. miR-135). In specific, miRNA sponging activity by linc-MD1 de-represses MAML1 and MEF2C leading to muscle-cell differentiation via the expression of the specific genes controlled by the latter. The HuR mRNA plays a subtle role in controlling the alternative processing of pre-miR-133b into linc-MD1 or miR-133. Most importantly, it competes with linc-MD1 for miR-133, thereby favoring the former’s sponging activity and giving rise to a positive feedforward loop that ultimately affects the levels of both species.

As discussed in [31], the trigger that possibly causes the system to exit the feedforward loop, repress muscle-specific gene expression and access later stages of differentiation is an endogenous upregulation of miR-133 transcription from the independent genomic loci. This suggests that the complex regulatory circuitry just described, whereby miRNAs can be synthesized both via a protein controlled switch and from an independent locus, can provide effective control of both timing and molecular levels.

In order to analyze this issue in a quantitative framework, here we analyze a schematic version of the above circuitry through a deterministic mathematical model based on mass-action kinetics, focusing specifically on the roles of HuR and of the alternative loci for miRNA transcription. In a nutshell, by characterizing the magnitude of the ceRNA effect and the response to a sudden increase of the transcriptional activity of miR-133, we show, among other things, that, while HuR-controlled regulation of pri-miR-133 processing is crucial to tune molecular levels, miRNA transcription from the alternative loci allows to achieve fast down-regulation of the target transcription factors MAML1 and MEF2C. In particular, fast enough miRNA-ceRNA binding kinetics causes non-linear response of the target level to upshifts in the miRNA biosynthesis rate as moderate as 20%, providing a highly efficient route to amplifying the effect of the differentiation trigger. In order for both mechanisms to be active, though, kinetic parameters need to be coordinated within specific ranges of values, which depend strongly on how sensitive pre-miR-133b processing is to HuR levels. In other words, the space of interaction constants and transcription rates is significantly constrained by crosstalk requirements.

The fact that crosstalk presupposes some degree of parameter tuning is not surprising per se, as strong ceRNA-ceRNA effective interactions at stationarity are known to be mainly achieved through competition when the concentrations of the involved molecular species are nearly equimolar [19, 20]. Remarkably, though, we find that the system’s dynamic behaviour in such ranges is completely compatible with that observed in time-resolved experiments. This supports the conclusion that the two regulatory elements of the myogenesis clock, namely the HuR-controlled miRNA-decoy system and the alternative locus for miRNA transcription, play different yet coordinated functional roles.

2 Results

2.1 Definition of the model

Our model is closely based on the mechanism controlling skeletal muscle cell differentiation identified and discussed in [29, 31] (see Fig. 1A and B). We consider a precursor RNA species (labeled $q$) which can be processed alternatively into a regulatory microRNA (labeled $\mu$) or a lncRNA (labeled $\ell$).

The relative weight of the two processing pathways is controlled by another RNA species (labeled $h$), such that larger values of $h$ increasingly favor synthesis of $\ell$ over $\mu$, thereby effectively repressing the latter. $\mu$, however, can directly repress $h$. In addition, $\ell$ can ‘sponge’ (i.e., transiently sequester) $\mu$ and lift its repressive effect from $h$. Finally, $\mu$ negatively controls a target mRNA (labeled $m$). Ultimately, $\ell$, $m$ and $h$ are ceRNAs, in that they compete to bind $\mu$. Note that, following [31], we allow $\mu$ to be synthesized independently of $q$ from a separate genomic locus. We shall refer to the two distinct miRNA biosynthetic pathways as the ‘HuR-controlled miRNA channel’ and the ‘alternative miRNA channel’, respectively (see Fig. 1A). By comparing against the HuR-controlled miRNA-decoy system of Fig. 1A, one sees that $q$, $\mu$, $\ell$, $h$ and $m$ play the roles of pri-miR-133b, miR-133, linc-MD1, HuR and MAML1, respectively (we neglected the remaining nodes for sakes of simplicity).

Our model includes the processes and rates reported in Table 1 together with the values used in this study.

Note that, since the exact mechanism of target repression by the miRNA is still unclear, we included both a stoichiometric and a catalytic decay mode for miRNA-ceRNA complexes. Indeed, while miRNAs incorporated into the RISC could return to the cytoplasm following target degradation, complexes may be stored in P-bodies, leading to their stoichiometric degradation [33]. As said above, $h$ regulates the values of the probability that $q$ is processed into $\ell$ (denoted by $\alpha$) and of the $(\mu,\ell)$ binding rate (denoted by $k_{\mu\ell}$). In [31], it is shown that the HuR protein downregulates the biogenesis of $\mu$ by repressing pri-miRNA cleavage by the enzyme Drosha, fostering accumulation of $\ell$. We therefore assume that $\alpha$ and $k_{\mu\ell}$ are Hill functions of $[h]$ (the level of $h$), with indices $n$ and $p$ respectively, i.e.

[TABLE]

with $h_{\alpha}$ and $h_{\mu\ell}$ the corresponding dissociation constants. For simplicity, we furthermore neglect miRNA-ceRNA complex dissociation rates, as partially justified by the observation that the ratio between the binding and unbinding rates of complexes is typically large [20]. Finally, in order to limit parameter proliferation and focus on the action of the HuR-controlled switch and of the alternative miRNA channel, we assume homogeneous degradation and complex processing constants for RNA molecules. (The point where biological realism is most sacrificed by this choice lies in the assumption that the intrinsic decay rates for miRNAs and ceRNAs are comparable whereas such rates are most likely linked to RNA length. This is however only going to affect our results in a quantitative way.)

With these choices, the mass-action dynamics of the regulatory element shown in Fig. 1B is described by

[TABLE]

2.2 Comparison with experimental time courses and choice of parameters

Time courses for the levels of linc-MD1 ($\ell$ in our model), miR-133 ($\mu$) and HuR ($h$) for a mouse myoblast cell culture capable of differentiation have been characterized in [31]. A Western blot showing the behaviour of MAML1 ($m$) during in vitro differentiation is instead displayed in [29]. These results indicate that 24 hours after the induction of differentiation the levels of $\ell$ and $\mu$ are low, while the controller $h$ is abundant. At this stage, $\mu$ is only synthesized from the precursor $q$ (i.e. from the genomic locus miR-133b). On the other hand, $[m]$ is at its maximum. In this phase, myoblasts differentiate and consistently express the target. 48 hours after induction, the level of $h$ peaks while that of the lncRNA $\ell$ appears to be increasing. This continues until about 72 hours after induction, during which time frame $[\mu]$ also increases thanks to a rapid increase in miRNA synthesis from the independent genomic locus miR-133a [31]. In this scenario, $\ell$ can sponge miRNAs away from the target $m$. When myoblast differentiation is accomplished (at about 72 hours after induction), though, $m$ must be repressed. Because $[h]$, and therefore $[\ell]$, appear to be decreasing between 72 and 96 hours, the sponge $\ell$ is no longer able to efficiently remove miRNAs from the target, which therefore gets silenced. This brings the cell to the next stage of differentiation, namely the formation of myocytes, where the expression of a different set of genes is paramount [8].

One sees that (i) $\mu$ is the agent that represses the target to complete the differentiation stage, (ii) $\ell$ plays a key role in controlling the timing of differentiation when the miRNA population increases following the activation of the alternative locus, while (iii) the controller ($h$) ensures that $\ell$ is synthesized at sufficiently high rates at the expense of $\mu$ via the $q$-processing switch.

We have modeled the increase in miR-133 production from the independent genomic locus through a step-wise perturbation on $b_{\mu}$ of the form

[TABLE]

where $\Delta$ measures the fold change in transcription rate after time $t^{\star}$ and $\theta(x)$ is the Heavyside step function. Using (10) together with the values of the kinetic parameters listed in Table 1, we obtain the scenario displayed in Figure 2, in which the concentrations of $h$, $m$, $\ell$ and $\mu$ are shown as functions of time.

Prior to the perturbation, which takes place at $t^{\star}\simeq 50$ h, the miRNA level is small while $[h]$ and $[\ell]$ rise together as the target level $[m]$ decreases. When the alternative channel for miRNA transcription sets in, $\mu$ levels shoot up, thereby repressing $m$, $\ell$ and $h$, in reasonable qualitative (for concentrations) and quantitative (for time scales) agreement with observations [31]. (Note that results have been obtained by fixing $\Delta=5$.)

For the remaining analysis, the choice of kinetic parameters was based on loose ranges extrapolated from [34, 35, 36, 37], namely

[TABLE]

where $j\in\{h,m,\mu,q\}$, $i\in\{h,\ell,m\}$, $\tau\simeq 10^{4}\;\mbox{s}$ is a typical mRNA half life in human cells and $\gamma\simeq 1\;\mbox{nM}$ a reference RNA concentration. The values of the kinetic parameters used in the numerical study of the steady state are shown in Table 1.

2.3 ceRNA crosstalk at steady state

We shall now quantitatively assess the ability of the miRNA-decoy system to control the expression of the target by analyzing both steady state and dynamical responses to perturbations. We focus specifically on the miRNA-mediated ceRNA mechanism, whose regulatory effectiveness has been quantified in [15]. Molecular levels at steady state are easily seen to satisfy the conditions

[TABLE]

Solutions (see Fig. 3A) show that a change in the transcription rate of HuR ($h$, the controller) with all other parameters fixed and no miRNA transcription from the alternative site affects the steady state concentrations of all other ceRNAs. For small enough $b_{h}$, the repressive effect of $\mu$ on the other RNAs is dominant, so the target ($m$) attains a low expression level. As $b_{h}$ increases, $[h]$ and, in turn, $\alpha$ and $k_{\mu\ell}$ increase according to (1), enhancing the synthesis of the lncRNA $\ell$ (the sponge), which gradually de-represses the target by sequestering miRNAs. As a consequence, the target level increases. For large enough $b_{h}$, $m$ and $\ell$ attain finite values while $h$ increases linearly (as expected) and the level of free miRNAs gets more and more suppressed.

Overall, Fig. 3A describes how competition for $\mu$ controls the expression levels of the competitors ($m$, $\ell$ and $h$) through modulation of the rate of synthesis of the controller and via the action of the $q$-processing switch alone (as $b_{\mu}$ is taken to be zero). Note that the regulatory element appears to be more sensitive to changes in $b_{h}$ in a relatively narrow range of values. In this range, molecular populations achieve near equimolarity, a condition under which the magnitude of miRNA-mediated crosstalk is maximal [20, 19].

Figures 3B–D display the dependence of $[m]$, $[\mu]$ and $\alpha$ on $b_{h}$ and on the rate $b_{\mu}$ of exogenous miRNA transcription, i.e. from the independent genomic locus miR-133a.

As $b_{\mu}$ increases, the window of values of $b_{h}$ where the ceRNA effect is stronger shifts towards higher values. Expectedly, though, the magnitude of target de-repression decreases significantly as more miRNAs become available (Fig. 3B), implying that effective ceRNA crosstalk requires a fine tuning of transcription levels from the alternative miRNA locus. Likewise, small values of $b_{h}$ lock the precursor processing pathway into synthesizing $\mu$ preferentially over the sponge $\ell$ (i.e. $\alpha\simeq 0$), see Fig. 3D. Note that the transition between the two states of the HuR-controlled switch is sharper when the transcriptional activity from the independent miRNA channel is low and becomes more graded for high values of $b_{\mu}$. This suggests that, in order to avoid step-wise transitions (which would impose a high extra metabolic burden on cells and could be affected by different types of noise, including cell-to-cell variability in transcription rates) as well as slow cross-overs (which would require large changes in transcriptional activity to robustly switch from one state to the other), both $b_{\mu}$ and $b_{h}$ should lie in an intermediate range of values.

2.4 Optimal control is achieved in specific ranges of kinetic parameters

In order to quantify the magnitude of ceRNA crosstalk and to understand how crosstalk requirements constrain kinetic parameters, we follow [19, 22] and focus on the ‘susceptibilities’

[TABLE]

which represent, respectively, the response of the target $m$ to a (sufficiently small) variation in the level of the controller $h$ or the rate of miRNA transcription from the alternative channel. Such quantities, similar to those employed for the study of magnetic spin systems in statistical physics, are especially useful for the analysis of ceRNA networks as they allow to focus specifically on the effects induced by competition to bind miRNAs. In principle, the existence of ceRNA crosstalk could also be signaled by a positive (Pearson) correlation between ceRNA levels. However, positive correlations can simply result from random fluctuations of miRNA levels: as both targets respond to changes in miRNA availability, their fluctuations are expected to be correlated when miRNA levels fluctuate stochastically. This however does not necessarily imply that the level of a ceRNA will change when the level of its competitor is altered. Quantities like (16), on the other hand, can be non-zero even for a deterministic system and therefore capture precisely the competition-induced aspect highlighted by the experimental literature. (See [38] for a detailed study of the relationship between susceptibilities and correlations in ceRNA networks.)

In order to calculate $\chi_{mb_{\mu}}$, we notice that (13) can be re-written as

[TABLE]

where

[TABLE]

is the stoichiometricity ratio quantifying the relative weight of the stoichiometric pathway of miRNA-ceRNA complex degradation. In turn, $[m]$, $[\ell]$ and $[h]$ take the form

[TABLE]

From (18) and (20) one derives a system of three equations for the susceptibilities $\chi_{mb_{\mu}}$, $\chi_{\ell b_{\mu}}$, $\chi_{hb_{\mu}}$, namely

[TABLE]

with $\bm{\chi}=\{\chi_{mb_{\mu}},\chi_{\ell b_{\mu}},\chi_{hb_{\mu}}\}$ and

[TABLE]

where

[TABLE]

Upon solving the above system one obtains explicit expressions for the susceptibility vector $\bm{\chi}$. Specifically, for $\chi_{mb_{\mu}}$, Eq. (17), one gets

[TABLE]

By similar steps, one can compute susceptibilities such as $\chi_{mh}$, Eq. (16), for which we find in particular

[TABLE]

The heat maps displayed in Fig. 4A and B represent $\chi_{mh}$ and $\chi_{mb_{\mu}}$ as functions of $b_{h}$ and $b_{\mu}$, with grey lines showing the curve of maximum $\chi_{mh}$ and minimum $\chi_{mb_{\mu}}$, respectively. (Note that while $\chi_{mh}$ is positive or zero since $h$ represses $\mu$ that in turn represses $m$, $\chi_{mb_{\mu}}$ is negative or zero, since a positive shift of $b_{\mu}$ causes the level of $m$ to decrease.) The key observation is that the miRNA-decoy system appears to achieve optimal target control in a restricted range of kinetic parameters.

Figure 4C and D detail these effects with a higher resolution. From the former one sees that for $b_{h}=10^{-5}\;\mbox{nM}\cdot\mbox{s}^{-1}$ (red curve) $\chi_{mh}$ is maximum for $b_{\mu}=0$ and it is monotone. In other terms, for a sufficiently small $b_{h}$, i.e. when the expression level of the controller is low enough, the target level appears to be maximally sensitive to the controller level when the alternative miRNA channel is inactive. As the controller’s level increases, instead, optimal control requires a non-zero transcriptional activity from the alternative miRNA locus. A similar scenario is obtained by analyzing the sensitivity of the target to $b_{\mu}$ upon changing $b_{h}$. Overall, optimal control is achieved in a restricted range of parameters characterized by $b_{h}\simeq 2.5\cdot 10^{-4}\;\mbox{nM}\cdot\mbox{s}^{-1}$ and small $b_{\mu}\gtrsim 10^{-5}\;\mbox{nM}\cdot\mbox{s}^{-1}$.

To sum up, we see that effective target control via $h$ can be achieved at steady state in absence of the alternative miRNA channel. However, if both the HuR-controlled channel and the alternative miRNA channel are active, then optimal regulation requires that their rates are coordinated (the faster the former, the faster the latter). On the other hand, effective target control by the alternative miRNA channel requires a controller. This provides quantitative support to the idea that the controller ($h$) ensures fine tuning at stationarity, when the alternative miRNA channel is off, while the latter appears to serve the mainly dynamical purpose of guaranteeing fast target downshifts, in full agreement with results obtained in [22], where a dynamical analysis of the ceRNA effect has shown that the optimal strategy to achieve fast target downregulation consists in rapidly upregulating miRNA levels.

2.5 The HuR-controlled channel dominates miRNA biosynthesis when ceRNA crosstalk is optimal

In order to further clarify how the interplay between the two modes of miRNA biosynthesis impacts ceRNA crosstalk, we analyze here how the adimensional quantity

[TABLE]

representing the ratio between the miRNA transcription rates from the independent genomic locus ($b_{\mu}$, corresponding to miR-133a) and from the precursor (corresponding to miR-133b), changes upon changing $b_{\mu}$ and the transcription rate of the controller, $b_{h}$. Figure 5 displays results, with the grey and green curves outlining, respectively, the $\beta=1$ and $\beta=2$ contours as a reference.

We see that the HuR-controlled channel is consistently dominant for small enough $b_{\mu}$ and $b_{h}$ (where ceRNA control is most efficient), while the alternative channel becomes increasingly more relevant as transcription rates increase. This again suggests that activation of transcription from the independent locus mainly plays a role in determining fast target downregulation. Notice that a two-fold increase of transcription rates with respect to optimal values for ceRNA crosstalk would be required in order to substantially alter this scenario.

2.6 The Hill index characterizing precursor processing sharpens the region of maximum ceRNA crosstalk

In Eq. (1), we have assumed that the parameters $\alpha$ and $k_{\mu\ell}$ are Hill functions of the HuR level $[h]$, and the corresponding Hill indices $n$ and $p$ have been so far kept fixed (see Table 1). Interestingly, though, they appear to have different effects on the emerging crosstalk scenario (see Fig. 6).

An increase of the Hill index $p$ controlling the sharpness of $k_{\mu\ell}$’s response to $[h]$ does not significantly influence the behaviour of the susceptibility $\chi_{mh}$. On the other hand, by increasing $n$ (which modulates $\alpha$’s response to $[h]$), the region of maximal susceptibility gets much sharper while preserving its overall qualitative behaviour. Notice that a larger $n$ would allow to stretch optimal control to much higher values of $b_{\mu}$, albeit at increased metabolic costs associated to the fine-tuning of $b_{\mu}$ and $b_{h}$ and without significant gains in terms of $\chi_{mh}$. Whether an increased sensitivity of the miRNA-decoy switch is advantageous therefore appears to depend crucially on the transcriptional activity of the alternative miRNA channel, with lower (resp. higher) levels requiring lower (resp. higher) sensitivity for $\alpha$.

2.7 Fast binding kinetics leads to non-linear response

The effectiveness of the mechanism allowing to exit the myoblast stage of differentiation can be analyzed by studying the change in the steady-state value of $m$ as a function of the perturbation size $\Delta$ (see Eq. (10)), quantified by the Integrated Response [22]

[TABLE]

where $t^{\star}$ stands for the time at which the perturbation is switched on. $\mathrm{IR}$ clearly depends on $\Delta$, as well as on the duration of the perturbation, and it is sensitive to the subsequent relaxational dynamics. In the context of miRNA-ceRNA networks, it has the extra advantage that it is directly linked to the total amount of protein synthesized from $m$ in response to the perturbation, which represents the key output variable in transcriptional regulatory systems. Note that, upon activating the alternative miRNA transcriptional channel with rate $b_{\mu}$, $[m]$ decreases monotonically, so that $\mathrm{IR}$ is always negative. Numerical results obtained for different values of the perturbation size $\Delta$ and for several values of the binding rates $k_{\mu m}$, $k_{\mu\ell}^{\max}$ and $k_{\mu h}$ are shown in Fig. 7.

For small values of the binding rates (see e.g. orange and green curves) $\mathrm{IR}$ is roughly linear in $\Delta$. As the binding rates increase, though, a regime characterized by non-linear response sets in, as seen for instance from the red and purple curves. Notice that a 20% increase of the miRNA transcription rate from the alternative locus can suffice to deviate from a linear behaviour.

Therefore, a sufficiently fast miRNA-ceRNA binding kinetics may lead to a strong variation of the target level even for a modest change in the perturbation size, suggesting an energetically cheap mechanism for the regulation of the expression of a gene involved in a differentiation process, as $m$, that has to be turned on at a definite time in the differentiation program.

3 Discussion

We have defined and studied a minimal, deterministic mathematical model of the regulatory circuit that has been recently found to control the timing and expression levels in the early phase of myogenesis [31]. We aimed at understanding the roles played by the various mechanisms that appear to coordinate the expression of a target mRNA, including a HuR-controlled channel for the mutually exclusive biosynthesis of the target repressor miR-133 and of its sponge linc-MD1, and an externally-controlled alternative channel for the synthesis of miR-133 from independent genomic loci. In summary, our results indicate that HuR, by controlling the linc-MD1/miR-133 switch, provides essential fine tuning of target levels, while the rapid increase of miRNA transcription from the alternative locus contributes to the fast downregulation of target levels when myoblast differentiation is accomplished. The former mechanism exploits the ceRNA effect and competition to bind miRNAs to indirectly regulate target levels. This strategy has been shown to deliver optimal results at stationarity [19, 15]. miRNA upregulation has instead been found to be the optimal mechanism to rapidly decrease the level of a target [22].

Note that the selection of a miRNA-decoy system over a simpler miRNA-sponge pair, in which the sponge’s biosynthesis is unrelated to that of the miRNA, guarantees an optimal level of complementarity (hence a strong coupling) between the miRNA and its lncRNA interaction partner. As we have seen, this condition is crucial for the emergence of non-linear response to a rapid increase of transcriptional activity from the alternative locus. Therefore the presence of a miRNA-decoy system in this context ultimately provides a major advantage both in terms of performance and in terms of metabolic costs.

Revealingly, while static control could be achieved in absence of the alternative channel for miRNA transcription, HuR-based regulation alone would not be as efficient in downregulating the target, as that would require rapid degradation of the controller (i.e. of HuR) to shift miRNA levels up. Both control routes are therefore essential. In addition, the activation of the alternative transcriptional loci can lead to non-linear (negative) response in target levels already for a modest (ca. 20 %) increase of miRNA biosynthetic rates, providing further support to the idea that miRNAs are optimal dynamical down-regulators.

Our analysis of the complex circuitry controlling skeletal muscle-cell differentiation has been highly simplified on several fronts. In first place, we have focused on a deterministic model, ignoring intrinsic noise sources altogether. The central reason for this lies in the fact that the ceRNA effect as quantified by susceptibilities like (16) occurs even in absence of fluctuations. Most importantly, in such conditions competition is the only possible source of ceRNA crosstalk (if direct transcriptional dependencies are excluded, as in the present case), as correlations due to the fact that both ceRNAs may respond to stochastic fluctuations in miRNA levels (which do not necessarily signal a positive effective coupling between ceRNAs) are forced to be absent from the picture. In this sense, deterministic models like the one discussed here are ideal to evaluate the role of competition in miRNA-based regulatory elements. Extending the present analysis to include intrinsic (molecular) noise will however allow to quantify the effectiveness of control more precisely by employing different measures, including correlation coefficients [20] and information-theoretical capacities [15].

Likewise, it would be important to perform an exploration of parameter space (especially for miRNA-ceRNA interaction constants and complex processing rates) so as to shed light on how robust our conclusions are against kinetic changes. miRNA-based regulatory elements have indeed been shown to be able to exploit kinetic heterogeneities in order to further optimize the efficiency of control with respect to kinetically homogeneous systems [15]. While we have focused on the more conservative scenario here, further work in this direction would highlight which functional aspect would most benefit from an accurate tuning of parameters and how the circuit’s overall functionality (which encodes for a key developmental program) can be affected by specific perturbations.

Finally, our model has focused for simplicity on a reduced version of the experimentally characterized system, namely one in which a miRNA (miR-135) and a ceRNA (MEF2C) have been suppressed. By influencing the level of linc-MD1, miR-135 however indirectly affects the feedforward loop controlled by HuR. In a full-fledged model accounting for all molecular species shown in Fig. 1A, then, one should observe a modulation of the levels of linc-MD1 due to its increased sponging activity. As MEF2C is effectively controlled only by miR-135 and, indirectly, by linc-MD1, a down-regulation of the latter will suffice to repress both MEF2C and MAML1, which indeed have similar time courses [31]. Therefore, it is reasonable to expect that the overall scenario will be qualitatively identical to that discussed here, albeit with a quantitative shift of the relevant parameter ranges. In specific, changes in $\alpha$ (the probability with which the HuR-controlled switch leads to linc-MD1) and/or $b$ (the intrinsic synthesis rate of linc-MD1) might suffice account effectively for the suppressed molecular species. A more detailed study (involving more parameters) will provide the complete picture.

Acknowledgments

Work supported by the European Union’s Horizon 2020 research and innovation programme MSCA-RISE-2016 under grant agreement No 734439 INFERNET.