A combination of transcriptional and microRNA regulation improves the stability of the relative concentrations of target genes

It is well known that, under suitable conditions, microRNAs are able to fine tune the relative concentration of their targets to any desired value. We show that this function is particularly effective when one of the targets is a Transcription Factor (TF) which regulates the other targets. This combination defines a new class of feed-forward loops (FFLs) in which the microRNA plays the role of master regulator. Using both deterministic and stochastic equations we show that these FFLs are indeed able not only to fine-tune the TF/target ratio to any desired value as a function of the miRNA concentration but also, thanks to the peculiar topology of the circuit, to ensures the stability of this ratio against stochastic fluctuations. These two effects are due to the interplay between the direct transcriptional regulation and the indirect TF/Target interaction due to competition of TF and target for miRNA binding (the so called"sponge effect"). We then perform a genome wide search of these FFLs in the human regulatory network and show that they are characterizedby a very peculiar enrichment pattern. In particular they are strongly enriched in all the situations in which the TF and its target have to be precisely kept at the same concentration notwithstanding the environmental noise. As an example we discuss the FFL involving E2F1 as Transcription Factor, RB1 as target and miR-17 family as master regulator. These FFLs ensure a tight control of the E2F/RB ratio which in turns ensures the stability of the transition from the G0/G1 to the S phase in quiescent cells.


Models in detail
Gillespie reactions In the following we briefly introduce the reactions involved in the circuits considered in the main text that we used to simulate the models through the implementation of Gillespie's direct algorithm [43]. As discussed in the text, we assume a titrative interaction of the miRNA with both its targets while the interaction between transcription factor (TF) and target (T) is modelled as a Hill function.
Set of reactions for the NM1: Set of reaction for the NM2: Set of reaction for the maimed NM4: Set of reaction for the micFFL: 2 Master equation In the following we report the master equations for the various circuits analyzed, with {s → n 1 , m 1 → n 2 , c 1 → n 3 , p 1 → n 4 , m 2 → n 5 , c 2 → n 6 , p 2 → n 7 } and the step operator E k j = ∞ l=0 k l l! ∂ l ∂n l j as defined in [32]. The rates are rescales as in the main text.
From each master equation we obtained the first two moments of the various distributions through the linear noise approximation (LNA) [32]. Briefly, given R reactions and N molecular species, a general master equation takes the form where Ω is the volume of the system, v ij is a matrix with dimensions # reactions · # species (in which each element is equal to the concentration change of species j due to reaction i) and E is the step-operator.
With the following change of variables n = Ω x + Ω 1 2 ξ (13) and after its substitution in the master equation, we obtain From the highest order (Ω 1 2 ) we obtain the equations for the mean concentrations per unit of volume. From the order Ω 0 and under the assumption that Ω is large enough (so that smaller orders are negligible) we obtain the equations for the fluctuations. Expanding the functions f i in power of ξ we have to consider only the 0th-order: Then, defining the two matrices and using the multivariate Fokker-Planck equation's solution for the first two moments, we obtain the following relationships: From these equations it is straightforward to calculate the approximated correlations r x,y and variances σ x for each concentration. In supplementary figure 1 we show the comparison between simulations and approximated analytical results.

Linear noise approximation and simulations: steady state discrepancy
The main problem of the above analysis is that the LNA does not properly take into account the correlation between chemical species in the evaluation of steady state concentrations. Indeed, if we solve the generating function for the first moments for mRNA and protein we obtain: The LNA expectation for the mixed moments is < xy >=< x >< y > while from the Master equation we find < xy >=< x >< y > +σ x σ y r xy = (1 + η x η y r xy ) < x >< y > where ηs are measures of noise. Using equation (20)  (1+ηsηm 2 r (s,m 2 ) )<s> < p 2 >= k p2 < m 2 > As it is possible to notice in Figure 2, these corrections can well explain the discrepancy occurring between simulations and approximated first order differential equations. Notice that the highest values of the correlation are obtained in the parameter region in which concentrations of microRNA and messengers are similar.

Comparison between NMand micFFL
To better investigate the increase in correlation due to the transcriptional link, we compare NM4 and micFFL (i.e. the two circuits which show the best results in terms of correlation between TF and T). We fix mRNA, protein, microRNA and complexes amounts and then evaluate the ratio between p 1 and p 2 correlations in NM4 and micFFL. The constraint is thus the following: Thanks to this constraint it is possible to compare directly the gain of correlation due to micFFL with respect to NM4. Our results are shown in Figure 3. The circuit with the transcriptional link between TF and T proteins always reaches higher values of correlation.

Comparison with the explicit promoter dynamics
In all the models used we assumed that the promoter dynamics is much faster than the other reactions. In order to test if such common approximation does not affect our results we substitute the reaction for the target mRNA production with three reactions describing the promoter interaction with the transcription factor (see Table below).

Reaction
Rate P icture From the reactions in Table 23 it is possible to calculate the probability of having the TF bound to the promoter. The activated promoter (p bound ) can then produce mRNA with a fixed rate: Comparing the two models (without and with the explicit promoter), we find the following mapping: (25) Figure 4 shows that for fast promoter dynamics there is a very good agreement between the two models in terms of mean number of molecules ( Figure 4A and B), noise ( Figure 4C) and correlation coefficient ( Figure 4D). Each curve is the mean over 1000 Gillespie's trajectories.

Stability analysis
In this section we accomplish a brief stability analysis of the micFFL circuit. From the definition d dt we expand the vectors for the micFFL case, so that Depending on the sign of ∇· − → F the system could be dissipative (∇· − → F ≤ 0) or conservative (⇔ ∇· − → F = 0). In our particular case the system is dissipative: Since the divergence of − → F is always negative (the parameters have to be always positive to make sense), the phase-space is not conserved and a 7-dimensional attractor exists. Further information on the stability of the system can be obtained if one fixes the values of the parameters in the equations. For example, with the set of parameters we used in our numerical experiments (k m1 = 23.5, k m2 = 41, k s = 17.5, k on c1 = k on c2 = 1.5, k p1 = k p2 = 117, γ m1 = γ m2 = 2, γ s = γ p1 = 1, αγ c1 = αγ c2 = 1.5, (1 − α)γ c1 = (1 − α)γ c2 = 1, h = 200, n = 1), the system has a linearly stable fixed point.

Steady state analysis with the logic approximation
We discuss here the steady state analysis of the Eq.(1) of the main text in the framework of the logic approximation in which the Hill function is approximated with the Heaviside step function f (p 1 ) = H(p 1 − h). Even if this approximation is very crude it may help to have an intuitive picture of the behaviour of the various players of the circuit as a function of the parameters of the circuit. The equations for p 1 and p 2 (see Eq. (1) of the main text) can be solved immediately leading to the steady state values: p 0 1 = k p1 m 1 /γ p1 and p 0 2 = k p2 m 2 /γ p2 . We can rescale the activation coefficient h s ≡ hγ p1 /k p1 and then write the step function as a function of m 1 and eliminate p 1 from the equations. Following [16] we introduce the quantities (i = 1, 2) • In the λ → 0 limit we find for the transcription factor m 1 the same threshold behaviour discussed in [16] as a function of the miRNA concentration. The same effect should be present also in the target concentration m 2 , but is hidden by the fictious step behaviour due to the logic approximation. It is easy to understand the origin of this threshold behaviour: if the number of free miRNA molecules greatly exceeds the number of transcripts mTF and mT, then these will be almost all bound in complexes and the corresponding proteins will not be expressed. On the opposite side, if the number of mTF and mT molecules overcomes miRNA amount, then nearly all miRNAs will be bound in complexes but there will be a sufficient amount of free mTFs and mTs to be translated.
• As the total miRNA concentration decreases the TF concentration increases following the trajectory plotted in fig 5. When the TF concentration reaches the threshold h s for the m 2 activation we observe a sudden enhancement in the TF concentration due to the sponge interaction between m 1 and m 2 . When also m 2 is present then the two mRNAs start to compete for the same miRNAs and as a net effect there is a smaller amount of miRNA available to downregulate m 1 . This non linear behavior of the TF concentration as the miRNA concentration increases is in our opinion one of the most effective ways to detect sponge-like interactions.
• The ratio m 2 /m 1 (and thus p 2 /p 1 ) can only take two possible values: m 2 /m 1 = m 0 2 /m 0 1 for m 1 > h s and m 2 /m 1 = 0 for m 1 < h s . However this is clearly an artifact of the logic approximation.