On the role of extrinsic noise in microRNA-mediated bimodal gene expression

Several studies highlighted the relevance of extrinsic noise in shaping cell decision making and differentiation in molecular networks. Experimental evidences of phenotypic differentiation are given by the presence of bimodal distributions of gene expression levels, where the modes of the distribution often correspond to different physiological states of the system. We theoretically address the presence of bimodal phenotypes in the context of microRNA (miRNA)-mediated regulation. MiRNAs are small noncoding RNA molecules that downregulate the expression of their target mRNAs. The nature of this interaction is titrative and induces a threshold effect: below a given target transcription rate no mRNAs are free and available for translation. We investigate the effect of extrinsic noise on the system by introducing a fluctuating miRNA-transcription rate. We find that the presence of extrinsic noise favours the presence of bimodal target distributions which can be observed for a wider range of parameters compared to the case with intrinsic noise only and for lower miRNA-target interaction strength. Our results suggest that combining threshold-inducing interactions with extrinsic noise provides a simple and robust mechanism for obtaining bimodal populations not requiring fine tuning. We furthermore characterise the protein distributions dependence on protein half-life.


I. INTRODUCTION
Environmental fluctuations can be source of noise in molecular networks, usually referred to as extrinsic noise [1,2]. Extrinsic noise, together with intrinsic fluctuations due to the probabilistic nature of chemical reactions, shapes gene expression and thus cell decision making and differentiation. Cells fate is normally triggered by non-trivial patterns of chemicals, such as morphogens, which are differentially expressed throughout the tissues. Cell phenotypic variability can be therefore measured by quantifying the underlying distribution of the triggering chemicals or of the downstream proteins in the same pathway. Bimodal distributions turn out to be a common outcome of gene-expression data, ranging from cancer to immune cells [3,4]. The two peaks of the distribution are usually associated with different physiological states of the system, may they refer to different organ primordia or different disease states or cancer subtypes [5][6][7][8].
In the past, some studies highlighted the possibility that microRNAs (miRNAs), similarly to their bacterial counterpart [9], may induce bimodality in the expression of their targets thanks to their specific titrative interactions [10][11][12].
MiRNAs are small molecules of non-coding RNA found in eukaryotes to act as posttranscriptional regulators. Although they were found in several different eukaryotic kingdoms, their role is known to be vital in multicellular organisms. They perform this function by recognising mRNA targets through Watson-Crick base pairing. Once bound to the target, they can both enhance the instability of that mRNA by degrading it and decrease its translation by keeping it bound. Interestingly, different levels of interaction miRNA-target can be achieved by different numbers of miRNAs [13][14][15]. Theoretical predictions together with in vitro single-cell experiments suggested that bimodality in the expression level of miRNA targets can be achieved with a high miRNA-target interaction strength [11,12]. In terms of genetic sequences, this would imply a high specificity between target and miRNA and therefore a high number of complementary binding sites (bs) per target.
As long as a miRNA molecule is bound to the target, such target cannot be translated.
It is then possible to define a threshold for the mRNA transcription rate such that below the threshold all the mRNA target molecules are bound to miRNAs and above the threshold there are molecules of mRNA free for translation [11,16,17]. Close to the threshold the number of both free miRNAs and targets is small, their fluctuations are highly coupled thanks to the non-linear interaction between the two and a small fluctuation in their amount of molecules may lead the system from the "bound" to the "unbound" state [11,12,17]. If the interaction strength between miRNA and target is high, then the transition from the bound to the unbound state is sharp. Close to the threshold, simply because of intrinsic fluctuations in the amount of both miRNA and target, part of the targets will be bound to the miRNA and part unbound for the same transcription rate. Picturing this in terms of target distribution would lead to a bimodal distribution whose two modes are associated to the bound and unbound state. It is worth to underline that this kind of bimodality is due to the presence of noise and not to peculiar molecular mechanisms introducing multiple deterministic stable states in the system.
MiRNAs are predicted to regulate more than 60% of our genome through a combinatorial action: every single miRNA can regulate several targets and one target can be regulated by different miRNAs [15,18]. The variety of targets they regulate is so wide and important for different signalling pathways or developmental stages [19,20] that the alteration of their expression levels may contribute to diseases such as tumour development and metastatisation [21][22][23][24]. It is nowadays also well established that multiple cell-cycle regulators are controlled by miRNAs, whose regulation could be in turn cell-cycle dependent [25][26][27][28]. The expression level of miRNAs may thus change with the cell-cycle progression, and there are indeed miRNAs differentially expressed according to the particular phase of the cell cycle [29]. As a consequence, in a population of cells heterogeneous with respect to the cell cycle, such as non-quiescent cancer cells, the amount of miRNAs can strongly fluctuate from cell to cell.
This introduces an extra source of noise in the system besides the intrinsic stochasticity of chemical reactions involving gene transcription, translation and regulation.
Here we study, both analytically and with simulations, how this extrinsic source of noise influences miRNA-mediated regulation. Although extrinsic noise affects also the expression of the other species in the network, we are here interested in specifically understanding the effects of a noisy source on miRNAs. We show how a distribution of miRNA transcription rates reshapes the threshold between miRNA and target and defines a wide region of bimod-ality. Such a bimodal distribution can be seen at a "population level", being the amount of miRNA heterogenous throughout the different cells. This outcome is deeply different than previous results where differential phenotypic expression is induced by the strong coupling between miRNA and its target at the "single-cell level". Interestingly, we also show that, if the miRNA target is protein coding, the protein half-life can alter the protein distribution.
With respect to the shape of the mRNA distribution, an increased protein half-life leads to a narrowing of the protein distribution around its mean. This may promote or suppress bimodality, suggesting that bimodal distributions at the level of mRNA may still correspond to a specific single phenotype at the protein level. Conversely, repressed heavy tailed mRNA distributions may give rise to bimodal protein distributions.
Finally, given the existence of multiple targets competing for one type of miRNA, we ask whether these properties can be maintained in a more complex circuit made of two competing endogenous RNAs (ceRNAs) and one miRNA [30]. The different target genes act indeed as sponges for the miRNA molecules and may sequester them from the environment. As a result, the overexpression or underexpression of one of the targets can lead respectively to an increase or decrease in the level of expression of the other competitors. The intensity of such cross regulation depends on the distance from the threshold of quasi equimolarity between miRNAs and targets [11,17]. This suggests that, if one target has a bimodal distribution, such bimodality may be influenced by the expression levels of the other miRNA competitors.

A. Stochastic model for miRNA-target interaction with extrinsic noise
Models of microRNA-mediated circuits have been the subject of several studies over the past years [16,34,35]. Hereby, we will refer to one of the simplest ways of accounting for microRNA-driven inhibition, depicted in Figure 1A. The molecular species involved in this circuit are miRNAs (S), target mRNAs (R) and proteins (P), result of the translation of the target mRNA.
In the following, we shall assume miRNAs and mRNAs to be transcribed from independent genes. For simplicity we neglect all the intermediate reactions leading to the synthesis of mRNAs, therefore assuming they are produced at constant rate k R . For the miRNA we consider it to be transcribed with a constant rate k S which we let fluctuate between different cells to probe the effects of extrinsic noise on the system. MiRNAs and mRNAs can also be degraded by the action of specialised enzymes. Hereby, we assume these reactions to be governed by mass-action law with rates g S and g R . The associated molecular reactions read: MiRNAs act as post-transcriptional regulatory elements, by binding the target mRNAs in a complex that can be subsequently degraded. Such interaction between miRNAs and mRNAs is quantified by the effective parameter g, which takes into account the strength of the coupling miRNA-target: from a biochemical point of view, it depends on the number of miRNA binding sites dedicated to a specific target [16]. The formation of miRNA-mRNA complex reads: While the mRNAs are always degraded due to the titrative interaction, the miRNAs can be recycled with probability 1 − α in the following way: Whenever the mRNAs are not bound to miRNAs, they can be translated into proteins with translation rate k P and, as assumed for the other molecular species, proteins can be as well degraded with rate g P , i.e.: From now on, we define as "intrinsic noise" the fluctuations due to the stochasticity of the chemical reactions with constant rates ( Figure 1A) and as "extrinsic noise" those due to the fluctuating miRNA transcription rate (see Figure 1C).
The system can be described by the probability distribution P (n S , n R , n P , t|K) of observing n S molecules of miRNA, n R molecules of mRNA and n P proteins at time t given a set of parameters K = {k R , k S , k P , g R , g S , g P , g, α}. This probability distribution follows the same master equation presented in [11] that can be either solved numerically or at the steady-state with some approximations. If the parameters fluctuate, in order to obtain the full distribution at the steady state P ss (n S , n R , n P ) one has to take into account such fluctuations. This can be achieved by using the law of total probability [32], which states that P (n S , n R , n P ) = P (K)P (n S , n R , n P |K)dK. As our aim is to test the effects of a fluctuating miRNA transcription rate, we shall assume this to be the only parameter drawn from a probability distribution, specifically a Gaussian centred around k S with variance σ k S and defined only for positive values of k S .
To obtain the steady-state distribution P (n S , n R , n P |k S ) conditioned on a specific miRNA transcription rate we could chose different approximation methods. Pivotal examples are the Van Kampen [31] and the gaussian approximations [11]. In the following we focused on the first one, leaving to the Supplementary Information (SI) a comparison between the two methods. We therefore performed a system-size expansion, thus assuming the system distribution at fixed parameters to be gaussian.
The marginal distribution P (n S , n R , n P ) is then found by using the law of total probability, i.e., by performing a weighted average on all possible values of k S .
Finally, we applied the same approach when considering two targets interacting with the same miRNA ( Figure 5A). In this case the conditional distribution is P (n S , n R1 , n R2 , n P 1 , n P 2 |k S ) from which one can obtain the full one by integrating over the values of the miRNA transcription rate.

Analytical approach
The equations governing the dynamics of the system considered in Figure 1 are given by: where R, S and P are the concentrations of the three species involved in the circuit expressed in nanomolars and the parameters are the same as in Figure 1.
Because of the inherent stochastic nature of molecular reactions, intrinsic noise should be taken into account by defining the probability distribution of observing n = (n R , n S , n P ) molecules at time t, namely P (n, t). The number of molecules of the species X, n X , relates to the concentrations, ρ X , as n X = V cell ρ X . The dynamics of this system can be rewritten  in terms of the master equation, that reads: Solving the master equation or even only determining the first and second moments of P (n, t) might be a difficult task due to the non-linear terms appearing from the miRNA-mRNA interactions. It is common practice to refer to approximations, amongst all the van Kampen's system-size expansion [31].
The van Kampen's system-size expansion relies on the assumption that the number of molecules of any species X can be split into two contributions, the former coming from the deterministic system and the latter from the intrinsic noise in the system, namely: where ξ X is a gaussian-distributed variable with zero average. The cell volume, V cell , represents the system size and it is assumed to be sufficiently large. Rewriting the master equation by using Eq. 7, at the leading order in V cell one finds the probability distribution to be a Dirac's delta in terms of the number of molecules, therefore recovering the deterministic equations Eqs. 5. The next-to-leading order in the expansion instead leads to a linear Fokker-Planck like equation in terms of the noisy variables {ξ X }, the probability distribution of the number of molecules then being gaussian centred around the deterministic concentrations φ X and with finite variance. By using this equation, one can compute the variances and the cross-correlations of the noisy variables, namely ξ 2 X and ξ X ξ Y that can be related to the molecules' variances and cross-correlations.
In the specific case of the circuit we analysed, the variances and cross-correlations are steady-state solutions of the following system: Using the expressions obtained from the previous system, specifying all the parameters, one can compute the probability distribution of the system, P (n, t).
To investigate the role of environmental fluctuations (i.e. extrinsic noise) in our system of interest, we consider a fluctuating rate of miRNA production k S . For the sake of simplicity, we assume this parameter to be drawn from a gaussian distribution, P (k S ) = , defined for positive values of k S , where k S and σ k S are respectively the average and variance of k S .
Due to the further stochasticity introduced by the extrinsic noise, the master equation previously derived does not hold anymore. Fluctuations on the parameter k S do not allow to rewrite a similar equation for this system in a simple way. However, as also shown in [32], the probability distribution of the entire system, P (n) can be rewritten in terms of conditional probabilities by using the law of total probability in the following way: where P (k S ) is the gaussian distribution in k S and P (n|k S ) is the conditioned probability of observing a certain configuration of the system, n, given a specific value of k S . This probability distribution is solution of the master equation Eq. 6 for any given k S . One can therefore apply again the van Kampen's expansion on the master equation of P (n|k S ) and obtain all the moments of this distribution (that will be all functions of the fluctuating parameter k S ). The full solution can be obtained by averaging the result over all the values of k S .

Molecular simulations of microRNA mediated circuits via the Gillespie's algorithm
Stochastic simulations have been performed via implementation of the Gillespie's direct algorithm [33]. All the results presented in this paper are obtained at the steady state.

A. Single cell versus population-induced bimodality
As discussed in the Introduction, the understanding of bimodal distributions is usually related to cell fate determination and differentiation. These mechanisms are at the basis of organism development and mis-development. It is therefore important to address the question about which might be the underlying molecular mechanisms allowing cell diversity and variability. Given the strong involvement of miRNAs in developmental stages, we focus here on the miRNA network represented in Figure 1A, at the single-cell level. Previous works [16,30] showed that the binding and unbinding reactions between miRNA and target give rise to non-trivial threshold effects in condition of quasi equimolarity between miRNA and target, see Figure 1E, where the threshold is defined in terms of miRNA and mRNA transcription rates as k * S = αk * R [16]. If the miRNA is transcribed at a rate above the threshold value, k S > k * S , the system is enriched in microRNA, which tends to bind most of the present mRNA therefore preventing its translation. In this regard, the system can be seen as below the threshold with respect to the target and we shall refer to it as in the "repressed state".
Above the threshold, the mean amount of free target increases linearly with its transcription rate. The scenario with free mRNA molecules will be denoted as the "unrepressed state" of the system. Such threshold effect displaying a rather sharp transition between the repressed and the unrepressed state gets more marked as the interaction strength between miRNA and targets increases. Close to the threshold value of the target transcription rate, due to the probabilistic nature of chemical reactions, the system will stochastically switch between the repressed and unrepressed state. This stochastic switching is enough to give rise to bimodal target distributions which appear for a narrow range of the target transcription rate k R [11] ( Figure 1B). Such bimodality characterises the single cell where the miRNA network is defined: every single cell can jump from the repressed to the expressed target state if the coupling constant with the miRNA is high enough.
On the contrary, in presence of extrinsic noise, the miRNA transcription rate is not the same for every cell ( Figure 1C,D). Hereby, we model the extrinsic noise through a gaussiandistributed miRNA transcription rate. To understand intuitively the consequences of this kind of extrinsic noise, let us consider the case of a miRNA transcription-rate distribution with fixed average k S . When the mRNA transcription rate (k R ) is very low and the average miRNA transcription rate is much larger than the threshold value ( k S αk R ), most of the drawn transcription rates k S will be larger than the threshold value. This would place the network in the parameter range where the targets are all bound to the miRNAs (see Figure 2). For larger k R , approaching the threshold, values of k S extracted from the left-tail will correspond to the case with some unbound targets. Below the threshold, as k S < αk R , the majority of the drawn k S will belong to the unrepressed state with the right tail of the distribution sampling from the all-bound region (Figure 2A,B). This scheme however is the same as before at a population level. The presence of rates above and below the threshold across the population can give rise therefore to a bimodal distribution in the number of free targets (Figure 2A,C). In particular, the higher the extrinsic noise, the larger the range of target transcription rates for which bimodality is present and the greater the separation between the two phenotypes (bound and unbound targets), see Figure 3A. This implies that, in contrast to the case without extrinsic noise, it is no longer necessary to fine tune the transcription rates to obtain a bimodal distribution. Even for high values of k R , the fraction of randomly picked k S that brings the system in the bound state is not negligible and forms a visible peak in the distribution (Figure 2A-C). The bimodal distribution is in this case given by the superposition of unimodal distributions obtained for different k S and weighted by the probability P (k S ) (see Figure S2 in SI). Focusing on one particular value of the variance σ k S and varying the target transcription rate, k R , we monitored the appearance of bimodal distributions. We ran Gillespie's simulations from which we sampled the number of targets for the histograms shown in Figure 2. By using system-size expansion and the law of total probability, we obtained the analytical targets distributions, shown in Figure  2. The analytical approximation captures the behaviour of the system for the mean, the Coefficient of Variation and the probability distribution, as testified by the agreement with the simulations (see Figure 2A,C,D and S3 in SI). Since in this case the bimodality is a single-cell effect, only those cells having the target highly interacting with the miRNA have chance to experience the repressed and unrepressed state when k R is close to its threshold value.
Adding some extrinsic noise relaxes the constraint on the interaction strength. Bimodality becomes rather a population effect, with some cells being locked in the repressed state (by having large miRNA transcription rates k S ) and some others (with smaller k S ) displaying free targets. Figure 3B shows

Interplay between different targets increases the stability of bimodal phenotypes
The study so far led us to picture the possible importance of extrinsic noise in cell phenotypic variability. Given the existence of multiple miRNA-targets networks, we now investigate how the results for the one-target case, extend to the multiple-target one. Let us consider a minimal model with two targets, R 1 and R 2 , competing for the same miRNA, S ( Figure 4A). We start by investigating the effect of an increase in the expression of target R 2 on target R 1 with and without extrinsic noise. Upon increasing the transcription rate k R 2 of target R 2 , the threshold of the target R 1 shifts towards lower expression levels: the miRNAs are indeed sponged away by R 2 and a lower amount of R 1 is needed to overcome the threshold.
If R 1 has a high interaction strength g 1 with the miRNA, then the range of bimodality shifts towards lower expression levels as well. The width of the range of R 1 bimodality is determined by the interaction strength g 2 of the target R 2 with the miRNA. If g 2 >> g 1 , then the miRNAs are sponged away by the second target with such a high frequency that the net effect is a reduction in the amount of miRNAs available to target R 1 . This entails a shift not only of the k R 1 threshold value but also of the range of bimodality. If g 2 < g 1 , the second target R 2 interacts with low frequency with the miRNA, R 1 is slightly derepressed and the net effect on its bimodal distribution is a shrink of the range of transcription rates for which it is present. The emerging picture is that, for a given transcription rate k R 1 , it is possible to tune the distribution of target R 1 from monomodal to bimodal and from unrepressed to repressed and viceversa via the expression of target R 2 .
The presence of extrinsic noise makes such cross regulation possible also for cases with lower miRNA-target interaction on both targets. In Figure 4B we show the bimodality phase diagram for R 1 at a fixed interaction strength g 1 between miRNA and target R 1 , and for a fixed level of extrinsic noise (see Figure S4 in SI for a different noise level). The interaction strength is such that in case of pure intrinsic noise R 1 does not show bimodal distribution.
As explanatory example, Figure 4C shows that the peaks of R 1 distribution can be tuned towards the repressed or the unrepressed case decreasing or increasing the expression of a second target R 2 . Here, the two targets R 1 and R 2 are both coupled through the noisy miRNA with small interaction strengths. Also in this case, the same analytical approach as before gives good agreement between theory and simulations.
These observations suggest that even if the miRNA repression is low and diluted over a network of multiple targets, the noisy environment allows cross-regulation between ceRNAs at the population level (see Figure 5D, with the intersection between the bimodality phase diagrams for both targets for a fixed level of noise and miRNA interaction strengths).  ] target 1 bimodality region target 2 bimodality region g 2 =30 nM -1 min -1 g 2 =60 nM -1 min -1 g 2 =120 nM -1 min -1 g 2 =300 nM -1 min -1 no bimodality target 2 transcription rate k R2 [nM min -1 ] Figure 4: Competition between two targets of the same miRNA. (A) Reference circuit including extrinsic noise for the case of two genes competing for the same miRNA. k S is the miRNA transcription rate. k R1 and k R2 are the mRNA transcription rates and g R1 and g R2 are the mRNA degradation rates of target 1 and 2 respectively. k P 1 and k P 2 are the translation rates and g P 1 and g P 2 are the degradation rates of protein 1 and 2 respectively. g1 and g2 are the miRNA interaction rates with target 1 and 2. α is the fraction of miRNAs that are not recycled after binding to the mRNAs. (B)
The other parameters are as follows: k R1 and k R2 range from 0 nM min −1 to 5.1 × 10 −3 nM min −1 , k S = 1.2 × 10 −3 nM min −1 , g S = 1.2 × 10 −2 min −1 , g R1 = g R2 = 2.4 × 10 −2 min −1 , k P 1 = k P 2 = 6.0 min −1 , its unrepressed mRNA, that is those molecules not bound to miRNAs and free for translation. In presence of extrinsic noise the target mRNA can show bimodal distributions even without clearly having a double steady state coming from the deterministic system. We here investigate if this is the case for the protein distribution. A key factor to keep into account is the time scale of protein synthesis and degradation. In general, if the protein dynamics are fast, the protein distribution follows closely the one of the mRNA (see Figures 5A1,5B1,  5C1). Conversely, slower protein dynamics tend to filter out the intrinsic fluctuations of the mRNA and lead to narrower distributions (see the supplementary material for a detailed discussion of the case without extrinsic noise). This is a single cell effect, that is, for a given rate of miRNA transcription k s , the corresponding protein distribution gets more peaked as the protein dynamics get slower. Then, also the protein distribution subject to extrinsic noise tends to concentrate around its mean. This feature bears different consequences on the protein distribution shape according to the specific structure of the mRNA distribution.
If the mRNA distribution is bimodal, slower proteins will have a distribution condensing around their mean which is located close to the unrepressed peak. They will therefore preferentially display the unrepressed phenotype and may completely lose their bimodal structure (see Figure 5A2). Bimodality can persist for strongly bimodal mRNA distribution because the noise reduction mechanism is acting at the single cell level so that it cannot overcome the effects of the extrinsic noise (see Figure 5B1,B2).
For a repressed (unimodal) mRNA distribution the mode is far from the mean, so that narrowing around the mean implies the rise of a second (unrepressed) peak. For moderately slower dynamics (see Figure 5C2) the protein distribution may be bimodal and for even slower ones it will be unimodal close to its mean (see Figure 5C3).
Altogether these results suggest that slow proteins promote the expression closer to the mean of the corresponding mRNA distribution. This may be sufficient to remove the bimodal feature of the protein distribution or not depending on the interplay among the amplitude of the extrinsic noise, the coupling between target and miRNA and the transcription rates.

IV. DISCUSSION AND CONCLUSION
Previous studies pointed out the relevance of extrinsic noise in molecular networks in shaping cell decision making and differentiation [1,2]. Gene expression levels can be normally observed into bimodal distributions, whose modes correspond to different physiological states of the system [3][4][5][6][7][8]. In this work we addressed the question on the role of extrinsic noise in shaping bimodal gene distributions in the context of miRNA-mediated regulation, both with stochastic modelling and simulations.
MiRNAs may induce bimodality in the expression of their target genes simply due to peculiar titrative interactions [11,12]. In a system with pure intrinsic noise, such bimodal distributions can be observed in conditions of high miRNA-target interaction strength and for a small range of target transcription rates. Indeed, the binding and unbinding reactions between miRNA and target in condition of quasi equimolarity let the target jump from the bound to the unbound state, giving rise to bimodal distributions. This bimodality is observed at a single-cell level: every single cell can indeed switch from one state to the other so that at a given time part of the population will be "bound" and part "unbound".
We showed that introducing some extrinsic noise on the miRNA transcription widens the range of the target transcription rate for which one observes target bimodality.
In this case the bimodal distributions arise at the population level, made of several cells that are heterogeneous with respect to the miRNA expression and therefore amount. Hence the bimodality comes out from the superposition of those unimodal distributions describing the single cells, i.e., each one obtained for a different value of miRNA transcription rate.
Interestingly, a high miRNA-target interaction strength is not necessary to obtain a population-induced bimodal distribution. We showed that extrinsic noise and miRNA-target interaction strength act at similar levels with respect to the bimodality. The interaction strength between miRNA and target in our model takes into account possible different number of miRNA binding sites on the mRNA target sequence. Since the miRNA repression on a given target is usually small and possibly diluted over multiple targets, our results suggest that some extrinsic noise can compensate for a low interaction strength in order to obtain differentially expressed phenotypes.
Since every single miRNA may have many different targets that in turn compete for the shared pool of miRNAs, a change in the expression level of one of them may alter the expression of the others. If one of these targets is bimodally distributed, then the bimodality may be influenced by the expression level of the other miRNA competitors. We modelled this scenario considering two targets in competition for the same miRNA and showed that cross regulation is possible even in case of small miRNA-target interaction strength if some extrinsic noise is present. In particular, different targets may cross-regulate each other's bimodal distributions and their interplay is pivotal in stabilising the presence of single phenotypes. This suggests that even if the miRNA repression is low and diluted over a network of multiple targets, the noisy environment makes cross-regulation among them possible.
The aim of regulatory systems is to control protein expression. Concerning the effect of extrinsic noise on their distribution, we showed that depending on the time scales of protein synthesis and degradation, the protein distribution may suppress or amplify the bimodality at the level of mRNA.
Altogether our results suggest that the coupling between extrinsic noise and threshold behaviour represents an efficient mechanism to obtain bimodal phenotypes without the need of fine tuning the rates of reactions which was required for the case of intrinsic noise only.
However, questions on the biological relevance of extrinsic noise may arise, especially about the perks of this high variability of different cells within a same tissue. Through our analysis, we observed that in both one-target and multiple miRNA-targets circuits, even in presence of extrinsic noise the region of unimodality is still quite wide (see Fig. 3A, 4B and D): the system is still able to noise buffer the extrinsic fluctuations and channel one only final phenotype into the protein branch within a wide range of parameters. In the remaining part, where the extrinsic noise is not buffered, other downstream mechanisms here neglected could help taking cell fate decision and hence lead the cell in a repressed or unrepressed state. In this case, the bimodal distribution induced by extrinsic noise would work as "enhancer" of cell-cell variability.
A feasible strategy to test the experimental validity of these theoretical results involves building ad-hoc synthetic circuits made of miRNA targets tagged with fluorescent labels, as previously done in [12,36]  In the Main Text we showed how to obtain the steady state distribution P (n S , n R , n P |k S ) conditioned on a specific miRNA transcription rate by a system-size expansion [1]. Alternatively, one could resort to a gaussian approximation of the moments as in [2] which consists in closing the hierarchy of the equations for the moments by expressing the third one as a combination of means and covariances. The simplest way to actually produce the shape of the distribution P (n S , n R , n P |k S ) is to make the additional assumption that it is gaussian with averages and covariances given by the moment closure scheme. Both approximations can be convoluted with the distribution of k S to produce the complete distribution P (n S , n R , n P ).
In this section we compare the results of the two approximations. For the plots in Figure 2 and 4 in the Main Text the two approximations give rise to basically indistinguishable curves (not shown) in agreement with the simulation results. However, they may fail to predict subtle details in the bimodal distribution, such as the initial emergence of the unrepressed peak addressed in Fig. 5 in the Main Text. For such parameters, the two approximations differ more markedly (see Fig. S1). In general, the gaussian approximation is more accurate in predicting the shape of the protein distribution in case of high protein stability. This observation can be traced back to the fact that the gaussian approximation over-performs the van Kampen one in describing the mean value of mRNA and proteins (see Fig. S1 D). As discussed in the Main Text, proteins with slower dynamics narrow their distribution around their mean so that predicting the mean accurately becomes more relevant. In Figure   S1 we show a comparison between the two approximations in relation to the histograms represented in Fig. 5 in the main text.

LOSS IN CASE OF PURE INTRINSIC NOISE
The comparison of coefficient of variations (CV ) are common measures of noise buffering in molecular networks. Here we focus on the comparison between the coefficient of variations of the targets (CV R ) and the proteins (CV P ) in case of pure intrinsic noise. If there is noise buffering, i.e. the amount of noise is reduced by the production/degradation of proteins, the protein CV will be smaller than the target one. Conversely, if there is no noise buffering, the two CVs are exactly the same. To give a prediction on the values of the parameters for which noise buffering is present, one can use van Kampen's system-size expansion's solution and build up analytical expressions for the coefficient of variations of targets and proteins.
Even by fixing all parameters but miRNA transcription rates and proteins transcription and degradation rates, the expression for the coefficient of variations are quite cumbersome.
To understand qualitatevely the behaviour of the Coefficient of Variation as a function of the degradation rate of the proteins, one can fix all the other parameters and compute the expression of the CV. Imposing CV P − CV R ≤ 0 allows to quantify the range of values of g P to have noise buffering in the final protein channel. In particular we found that g P ≤ 0.9 min −1 for the parameters' choice as in Fig. 5A-A1, i.e., for values of g P smaller than a given threshold and hence for long-living proteins.   3.3 g 2 =30 nM -1 min -1 g 2 =60 nM -1 min -1 g 2 =120 nM -1 min -1 g 2 =300 nM -1 min -1 no bimodality FIG. S4: Bimodality phase diagram. The plot shows the bimodality phase diagram for the mRNA 1 in a system with two targets competing for the same miRNA. The parameters here used are the following: k S = 1.2 × 10 −3 nM min −1 , σ S = 2.4 × 10 −4 nM min −1 , g 1 = 1.2 × 10 2 nM −1 min −1 , k R1 and k R2 range from 0 nM min −1 to 5.1 × 10 −3 nM min −1 , g S = 1.2 × 10 −2 min −1 , g R1 = g R2 = 2.4 × 10 −2 min −1 , k P 1 = k P 2 = 6.0 min −1 , g P 1 = g P 2 = 2.4 × 10 −2 min −1 , α = 0.5.