Mesenchymal Stem Cells for the Treatment of Skin Diseases

1.   Introduction

Multistability, the simultaneous existence of two or more stable steady states, is a feature dramatically relevant for the functioning of living systems ^[1]. Deterministic multistability is realized in gene regulatory networks by the inclusion of specific topological features. Feedback loops for instance are known determinant of multistable behaviours ^[2,3,4] and underlie a number of processes involved in cellular decision-making, ranging from the LAC operon ^[5], to quorum sensing ^[6,7], to cell differentiation ^[8].

Exploration of multiple deterministic steady states is realized through stochastic fluctuations at the molecular level, or noise. Noise is ubiquitous in molecular biology, and occurs in many aspects of gene regulation and other cellular processes in bacteria ^[9], yeast ^[10], and mammalian cells ^[11]. By allowing transitions among alternative steady states, it accounts for much if not all of the variability that we see in biological systems ^[12,13,14].

In general, noise is classified in the two broad categories of intrinsic and extrinsic noise. Intrinsic noise originates in the natural randomness of molecular events, and affects directly the dynamical variables of the system, such as for instance protein or mRNA concentrations of a gene regulatory network. Its modelling is based on a wealth of results, which aim to solve the Chemical Master Equation ^[15], analytically for instance in terms of the Linear Noise Approximation ^[16], or by direct stochastic simulations, using the celebrated Gillespie algorithms ^[17,18].

Extrinsic noise in contrast represents the effect of a fluctuating environment on the system. It translates mathematically into stochastic fluctuations acting on the parameters that enter the system, such as pH levels, temperature, or rate and degradation constants. Extrinsic noise too has been largely investigated in biology in the last two decades ^[19,20,21].

Interestingly, extrinsic noise has been shown to produce highly non-trivial effects, responsible for the emergence of multistability or oscillatory behaviours. New mechanisms, purely stochastic, have now been identified as determinants of multistable behaviour, in addition and possibly in alternative to the selection of deterministic tolopogical features. These effects are referred to as 'purely' noise-induced transitions, where the term 'purely' is used to emphasize the difference with respect to transitions driven by noise among deterministically pre-exisiting steady states, and is often omitted in the literature once the context is clear. We will also omit it in the rest of this paper.

Two fundamentally different mechanisms are known to cause noise-induced transitions. The first mechanism refers to the fact that in presence of fluctuating environments the stochastic differential equation (SDE) describing the system contains multiplicative noise terms. If Gaussian white noise is chosen to describe the environmental fluctuations, the presence of these terms makes the problem mathematically ill-defined, and the solution of the corresponding SDE is not unique. To make the problem mathematically well-defined and tractable, the SDE requires a specific discretization prescription that defines its stochastic integral ^[15,22,23]. The Itô and the Stratonovich prescriptions have been put forward as sensible choices to compute a unique solution of the SDE describing the system, but which prescription is the 'right' one has been debated for a long time. It is now widely recognized that the requirement that the noise be charaterized by a small but finite correlation time implies that the Stratonovich prescription needs to be adopted in the corresponding white noise limit. We call this type of noise 'effectively white' or 'physical' noise, in that its frequency spectrum is characterized by a high frequency cut-off, as indeed happens for physical noises. The adoption of the Stratonovich prescription then leads to the emergence of a spurious drift, the so-called Stratonovich drift, that modifies the deterministic dynamics of the system, and gives rise to non-trivial properties and noise-induced transitions.

The emergence of noise-induced transitions relying on the Stratonovich drift has been the subject of intense theoretical and experimental research in a number of different physical and chemical systems, including homogeneous ^[22] and spatially extended systems ^[24], but also in ecological modelling ^[25], and in socio-economic systems ^[26].

In biology, several hypotheses that rely on Gaussian white noise and on the Stratonovich drift have been put forward to explain the emergence of bimodal behaviours in sytems deterministically monostable. For instance, Samoilov et al. ^[27,28] have shown that futile enzymatic cycles, deterministically monostable, can transition to a bimodal regime when noise affects the enzymes’ activities in an appropriate manner. Similar considerations have led to the formulation of the concept of stochastic control in metabolic pathways ^[29]. In development, Pujadas and Feinberg ^[30] have speculated that noise participates to the developmental program, by contributing to the process of cellular decision making through modifications of the morphology of the epigenetic landscape by noise-induced transitions. In bacteria, quorum sensing ^[31] and simple repressive circuits ^[32] have been investigated stochastically, and have been proven to show non-trivial emerging properties related to noise-induced transitions. In all these studies, Gaussian white noise has been assumed.

However, a second mechanism accounting for noise-induced transitions has been recently discovered, formalized theoretically, and verified experimentally. Ochab-Marcinek et al. ^[33,34] have proposed a neat geometric construction based on nonlinear noise processing. In its simplest version, this mechanism relies on the existence of a formal relation between two stochastic variables that mimics a nonlinear input-output response between them. Simple conservation of probabilities then establishes the output probability distribution in terms of the input probability distribution. Because of the assumed nonlinearity of the response function, it may happen that the output distribution becomes bimodal, even though the input distribution was unimodal.

This second mechanism can be formulated for any input probability distribution, not necessarily Gaussian, and for a variety of biochemical systems. For instance, Birtwistle et al. ^[35] have shown that bimodality of extracellular signal-regulated kinase (ERK) emerges as a response to epidermal growth factor (EGF) in presence of a fluctuating activation threshold that can be analyzed in terms of the nonlinear noise filtering methodology ^[33], as highlighted by Kim and Sauro ^[36]. Similarly, Dobrzyński et al. ^[37] have studied a stochastic model of the hypoxia-inducible factor (HIF) pathway, and found that it shows emergence of bimodal dynamics of HIF depending on response threshold variability to dimethyloxalylglycine, which mimics oxygen limited conditions. What is noticeable in these studies is that they provide experimental validation not only of the occurrence of bimodal dynamics, but also that these emerge as purely induced by noise. In a different context, the differential and biphasic tolerance of clonal bacterial cells to drugs, so-called bacterial persistence ^[38,39,40], has also been described in terms of nonlinear noise filtering ^[41,42].

It is interesting to note that the geometric construction presented in ^[33] can be derived in dynamical terms by extending to nonlinear noise a well established formalism, first developed in ^[43] and applied later in a variety of chemical and physical systems ^[44]. The formalism is based on relaxing the assumption of white noise, and on focusing instead on fluctuations characterized by a correlation time much longer than all other relaxational timescales in the system. We present here this derivation, which offers further insights into the nature of the related noise-induced transitions, and allows for a direct comparison with the Stratonovich approach.

Despite its success over the years, and independently of the mechanism responsible for the occurence of noise-induced transitions, the Gaussian white noise assumption is not free of criticism in gene regulatory networks. First, the Gaussian approximation is of limited applicability for fluctuations affecting parameters which are strictly positive. Second, the white noise assumption is also problematic in GRNs. From the biological perspective, the white character of extrinsic noise is largely questionable for gene regulatory networks, where the typical correlation times of extrinsic noise can extend up to and well above typical cell cycle times ^[45,46]. Third, regulatory dynamics described mathematically in terms of Hill functions include nonlinear dependencies on some of the parameters. This makes the mathematical treatement of the corresponding stochastic dynamics for white noise very hard, if not impossible ^[22].

These problems have been recognized in some recent literature, and relevant progress has been made in characterizing the response of different systems under more biologically grounded noise. Non-white noise, so-called coloured noise, characterized by a finite correlation time, can be addressed by the powerful approach of the Unified Colored Noise Approximation (UCNA) ^[47], which captures very well the dynamics of the system under the Gaussian approximation. Recently, Holehouse et al. ^[49] have applied UCNA to a genetic feedback loop, and a thorough analysis of the effect of Gaussian coloured noise has identified noise-induced transitions from uni- to bimodal behavior and vice versa. We review the UCNA approach in the next section, as a first step in addressing the effect of finite noise correlations under the Gaussian approximation.

In order to relax the Gaussian approximation itself, so-called bounded noise has been introduced recently. Bounded noise is generally obtained by starting from a a Gaussian process, and imposing hard bounds on it, so as to constrain it in a bounded domain. Examples are the Sine-Wiener process, the Tsallis noise, or the Borland noise ^[50,51,52]. These different types of noise have been studied in proof of concepts examples, such as the genetic model ^[53], and in applied set-ups ^[54,55]. Gamma-distributed noise and log-normal noise are also of marked interest in that they are particularly successful in fitting single cell experiments ^[56,57], and can be derived under very broad assumptions from the statistical properties of gene expression data ^[58]. They preserve by construction the positivity of the relevant variable, as their domain of definition is semi-bounded, extending from zero to infinity. All these noises have been found to produce noise-induced transitions, and in general a strong divergence from the expected Gaussian behaviour.

In this paper we review the two main mechanisms by which noise-induced transitions can emerge. First we present the standard derivation of the Stratonovich drift for Gaussian white noise and the UCNA approach for small but finite noise correlation time. Then we introduce the noise filtering approach, and show how this method can be derived dynamically by considering slow fluctuations. We consider both Gaussian and log-normal noise, and review the implications of both approaches for a simple but fundamental system, the repressed gene, analysing both linear and nonlinear dependencies. We finally discuss the main biological implications of these results, and highlight a series of open problems and challenging questions.

2.   The Stratonovich interpretation for Gaussian white noise and the UCNA approach

Let us consider the following dynamical system:

where $\lambda$ is a control parameter, related to environmental fluctuations, and which appears linearly in the dynamics specified by the function $f$. Let us assume that $R$ is affected by fluctuations

Here $\bar{\lambda}$ is the average of $\lambda$, $\xi(t)$ is Gaussian, zero-average, white noise, with intensity $D$ and correlator given by

Because of the linear dependency of $f$ on $\lambda$, we can decompose the fluctuations acting in (2.1) into two terms, and write:

where $g(x) = \partial f(x, \lambda)/\partial \lambda$.

Equation (2.4) is a multiplicative noise equation, and it is ill-defined as it stands. Its definition relies on the time discretization adopted to evaluate its stochastic integral, and it only acquires meaning if supplemented with an adequate discretization prescription ^[23]. Historically two prescriptions have been formulated and debated, the Itô prescription and the Stratonovich prescription, which are implemented by considering the two following versions of Eq (2.4):

Setting $\alpha = 2$ corresponds to assuming the Itô prescription, while setting $\alpha = 1$ corresponds to assuming the Stratonovich prescription. The term $D g(x) g'(x)$ in Eqs (2.5) and (2.6) is the so-called 'Stratonovich drift', and is responsible for the emergence of non-trivial dynamics, as we shall discuss. The associated Fokker-Planck Equation reads

Equation (2.6) can be solved for the stationary probability distribution ^[23]:

where $\cal{N}$ is a normalization factor. The modes of the stationary probability distribution correspond to the stochastic extension of the steady states of the system, and are given by

Equation (2.8) clearly shows that depending on the interplay between the functions $f(x)$ and $g(x)$, the modes can be dramatically different from the deterministic steady states of the system, identified by the equation $f(x) = 0$. This is in principle true for both the Itô and the Stratonovich prescription, as the Stratonovich drift appears in Eq (2.8) in both cases, albeit with a different weight. Its effect is to possibly change the number and the stability properties of the deterministic steady states, with the noise intensity $D$ acting as a bifurcation parameter, driving transitions between unimodal to bimodal or multimodal dynamics. Because these transitions are absent in the deterministic system, they are called noise-induced transitions. The noise is not only driving transitions among pre-existing deterministic steady states, but it possibly creates new ones, and drives transitions among them.

Whether the Itô or the Stratonovich prescription is to be assumed, has been extensively debated in the past. While the two prescriptions are mathematically legitimate and equivalent, of course the implications can be different on noise-induced transitions, as the $\alpha$-dependency in Eq (2.8) shows. For physical and chemical systems, the choice depends on the way noise is derived and on the way if affects the system itself ^[15]. For biological systems, similar considerations apply, with different conclusions being drawn according to the type of noise that is examined. For instance, intrinsic noise can be modelled effectively in terms of the Chemical Langevin Equation ^[59,60], which needs to be interpreted in the sense of Itô, and can be shown to lead to noise induced transitions in small gene regulatory circuits ^[32]. It should be emphasized that the appropriate definition of Langevin dynamics has to rely on the definition of the corresponding master equation, as is done in ^[59,60], whereas arbitrary insertion of noise to mimic intrinsic fluctuations leads to unphysical and inconsistent results ^[15]. In contrast, when extrinsic noise is considered, the assumption that white noise is the limit of Gaussian noise with small but finite correlation time, and therefore 'physical', requires the adoption of the Stratonovich prescription ^[61].

In order to explore further the behaviour of the system when driven by Gaussian noise with small but finite correlation time, let us consider the following dynamical system:

In Eq (2.9) the relevant timescale for $x$, $\tau_x$, is now explicity written for clarity. In contrast to the white noise case, with Eq (2.10) we now assume that the parameter $R$ is described by a Ornstein-Uhlenbeck process ^[23], characterized by a correlation time $\tau$, driven by the Gaussian, zero average, white noise $\xi = \xi(t)$, with correlator given by Eq (2.3). This system extends the system defined by Eqs (2.1)–(2.3) for finite correlation time of the noise. In fact, the correlator of the $\lambda$ process can be easily calculated, and results in

For $\tau \rightarrow 0$, the correlator (2.11) tends to the white noise limit, namely Eq (2.3), and the system defined by Eqs (2.9) and (2.10) becomes formally equivalent to Eqs (2.1)–(2.3).

For $\tau$ finite, the Fokker-Planck Equation has been derived within the so-called Unified Coloured Noise Approximation (UCNA) ^[47],

where

and whose solution is

This solution is exact in both the $\tau = 0$ and $\tau = \infty$ cases. An extensive discussion of this result, including its limitations, is presented in ^[47,48]. In the present context, we notice that it is immediate to verify that for $\tau \rightarrow 0$ Eq (2.12) reduces to Eq (2.6) in the Stratonovich interpretation, confirming thereby the general result ^[61].

This latter argument lends itself well to considering also the neighbourhood of white noise, or so-called 'virtually white' noise. As long as $\tau \ll \tau_x$, we can effectively consider the noise as white, and use the white noise Fokker-Planck Equation (2.6) instead of Eq (2.12), as long as we interpret it in the Stratonovich sense. In this respect, the white noise assumption, generically valid when the correlation time of the noise is much smaller than any other timescale of the system, simplifies greatly the treatment, and allows an effective approximate description.

3.   The repressed gene: Gaussian fast noise on degradation

The effect of Gaussian extrinsic noise, both white and coloured, on simple gene circuits has been object of recent studies. Holehouse et al. ^[49], for instance, investigate the effect of noise on degradation rates by using the UCNA approach in a feedback loop in the bistable regime and in concomitance with intrinsic fluctuations on protein concentrations. Previously, Shahrezaei et al. ^[62] had addressed similar issues, but in a single activated gene, under the effect of log-normal noise.

Here, we revisit these results in a simplified but fundamental model of gene expression, the repressed gene. In order to ascertain the effect of extrinsic noise only, we assume that intrinsic noise, related to protein low copy numbers, is absent, and compare the cases of white noise and coloured noise with small but finite correlation time.

The repressed gene is a gene which is constitutively expressed at its maximal expression rate, but it is regulated by a transcription factor (TF), the repressor, that represses protein synthesis (Figure 1).

We describe this system by decomposing the rate equation for the protein $x$ into a production term, described phenomenologically in terms of a Hill function, and a degradation term, straightforwardly written according to the Law of Mass Action:

Here $\rho$ is the association constant of $R$ to DNA, which describes the strength of the repressor binding to its binding site, $g$ is the maximal expression rate for gene $G_x$, and $k$ is the degradation rate of the protein $x$. Due to the much shorter half-life of mRNA with respect to that of proteins ^[63], we here assume that transcription and translation are lumped together in one single step, so that when the gene is in the active state it synthesises directly the protein $x$. Also, when the concentration of $R$ is very high, we can consider it as unaffected by the $x$ dynamics, and therefore as a constant control parameter. Hence, the dynamics specified by the function $f(x, R, k)$ is characterized by a timescale $\tau_x$, which in the specific case considered here is simply given by the protein inverse degradation rate, $\tau_x = 1/k$.

These dynamics are supplemented with the following biologically grounded parameter values, which broadly apply to typical bacterial gene expression. We estimate the effective maximal gene expression parameter in the range $g \approx g_P g_M/k_M = 10^{-2} \div 10^{-1}$ s$^{-1}$, where $g_P$ and $g_M$ are transcriptional and translational rates respectively, and $k_M$ is the degradation rate of mRNA. This choice is compatible with a moderate to strong transcriptional and translational activity and typical mRNA degradation rates in bacteria ^[64]. Also, we set $k = 10^{-3} \div 10^{-4}$ s$^{-1}$, which corresponds to a protein half-life ranging between minutes and a few hours ^[65], and a typical dissociation constant of trancription factors to DNA $\rho = 0.1 \div 10$ nM$^{-1}$ ^[66]. We also choose $R$ in the range $R = 1 \div 100$ nM, which together with the chosen $\rho$ gives on average a weak ($\rho R < 1$) to strong ($\rho R > 1$) repression activity of $R$ on gene $G_x$. Simulations in this paper are carried out with parameters with these orders of magnitude, and specific parameter values will be indicated in the captions of the figures.

Throughout this paper we will use the repressed gene as the preferred model to revisit the effect of different types of noises on different parameters, namely on the degradation rate $k$ and on the repressor concentration $R$. We do not consider fluctuations on $g$ as these would merely correspond to insert additive noise into the system.

Following ^[49,62], let us consider first white noise on the degradation rate:

where we assume $\xi(t)$ to be Gaussian and white, with correlator given by Eq (2.3). Since the degradation rate appears linearly in Eq (3.1), we can apply straigthforwardly the approach introduced in the previous section, and write the corresponding SDE as

which needs to be supplemented with the Stratonovich prescription as discussed, and where we have set $\beta = g/(1+\rho R)$ for convenience.

We can directly compute the normalized solution of the Stratonovich Fokker-Planck Equation (2.7) as:

whose mode is located at

as implied by Eq (2.8), and to be compared with the deterministic steady state located at

We show direct stochastic simulations of the system (3.3) in Figure 2. The agreement between the simulations and the analytical solution (3.4) is excellent, and it captures the predicted dynamics. Not only a broadening of the mode is observed for higher noise intensities, but a shift in the position of the mode itself takes place, which obeys the theoretical prediction (3.5).

We also explore the effect of coloured noise on the system, by adopting the UCNA approach presented in the previous section, and adopted in ^[49,62]. The UCNA solution (2.14) reads in this case

which correctly tends to Eq (3.4) for $\tau \rightarrow 0$. In Figure 3 we show the $\tau$ dependency of the UCNA solution for this system. The agreement between the direct stochastic simulations and the UCNA solution is excellent. It should be noted that the UCNA solution appears to be virtually indistinguishable from the white noise solution for $\tau \ll \tau_x \approx 1/k = 10^3$ s$^{-1}$, while for increasing values of $\tau$, getting closer to $\tau_x$, the UCNA solution and the white noise solution diverge substantially from each other. This implies that as long as $\tau \ll \tau_x$, the white noise approximation works well even if the noise is characterized by a finite timescale, what we have named the 'virtually white' regime.

We notice that despite the fact that the changes seen in probability distributions are limited to shifts of the position of the mode, and to its width, no change in the number of stationary states in the system occurs. This reproduces qualitatively the results ^[62], where a similar shift was observed when the system is driven by log-normal degradation noise, and no noise-induced transitions were reported. We will discuss this further in Section 6. Similarly, the bimodality emerging in ^[49] for a self-activating gene is related to the pre-existent bistability of the system, with the noise changing the way the different (deterministic) attractor states are reached, through excursions in the bistable parameter regime. These conclusions are in agreement with the observation that extrinsic noise acting on parameters appearing linearly in the system is not expected to modify its dynamics substantially, and to cause transitions from unimodal to bimodal behaviours ^[37].

4.   Non-Gaussian noises: log-normal noise as a special case

While the Central Limit Theorem would suggest the modelling of extrinsic fluctuations as a Gaussian process, this choice is questionable for fluctuations affecting a strictly positive parameter. The need of preserving the positivity of noisy parameters has brought along a number of approaches to model non-Gaussian noise sources. For instance, so-called bounded noises have been introduced recently and have attracted attention in abstract and more applied settings. A first example is offered by the so-called sine-Wiener noise, initially introduced in ^[51], and applied for instance in ^[54,55] in tumor biology. A second example is offered by the Tsallis-Borland noise, which with an adequate choice of parameters corresponds to a Tsallis statistics ^[50,53]. Other types of bounded noises, currently being investigated, include the Cai-Lin noise ^[67], and the Kessler-Sørensen noise ^[68].

A particular mention goes to both Gamma-distributed noise and log-normal noise. Gamma-noise emerges as intrinsic noise in simple two-stage models of gene expression ^[69,70]. Remarkably, Gamma distributions are robust to slow heterogeneities in parameter values ^[71], and therefore are a good fit to experimental data, integrating over both intrinsic and extrinsic noise sources, as shown in single cell experiments ^[56,71]. Similarly, log-normal noise is also a good candidate to represent fluctuations on parameter values phenomenologically ^[45,62]. It is supported by both strong experimental evidence ^[57] and theoretical arguments ^[58], and it has been recently investigated to interpret universal features in bacteria and yeast ^[72,73].

In this paper we will focus our analysis on log-normal noise. We define

where $\eta(t)$ is the standard Ornstein-Uhlenbeck noise, with the Langevin representation

with $\xi(t)$ zero average Gaussian white noise. The definition (4.1) guarantees that $\langle \lambda \rangle = \bar{\lambda}$, since $\langle e^{\eta(t)} \rangle = e^{D/2}$. The fluctuations on $\lambda$ can then be described by

and are characterized by the stationary log-normal distribution:

It should be noted that the log-normal distribution (4.4) is always unimodal. In particular $w(0) = 0$ for all parameter values.

5.   Nonlinear noise filtering and its dynamical derivation for slow noise

Ochab-Marcinek et al. ^[33] have recently introduced nonlinear noise filtering as a useful methodology to study the response of a gene under the effect of an external signal. In general terms, let us consider two variables, $x$ and $y$, such that $y = f(x)$. We can think of $x$ and $y$ respectively as input and output of a dynamical system, and the function $f$ as the response function. If the input variable $x$ is stochastic, with a probability distribution given by $p(x)$, this will induce a probability distribution $w(y)$ for the output variable $y$, thanks to the relation $y = f(x)$. The distribution for $y$ is readily computed by invoking conservation of probability:

which then implies:

where $x = f^{-1}(y)$.

Equation (5.2) is at the basis of the noise filtering approach. If the relation $f$ between input and output is linear, the two proabability distributions are trivially linked, but in case of a nonlinear response function $f$, the transformation (5.2) introduces qualitative changes in their functional forms. For instance, depending on the input distribution, highly nontrivial effects can be expected for the family of Hill functions

where $\alpha$ and $\beta$ are shape parameters. This is the form adopted in ^[33] to describe the response of a repressed gene to a regulating transcription factor. Similar considerations hold for signalling networks ^[37], and for recently proposed models of bacterial persistence ^[41].

In this section we will propose a dynamical derivation of the noise filtering approach when the correlation time of the extrinsic noise is much longer than any other timescale in the system. Extrinsic noise is indeed characterized by slow fluctuations, showing memory and correlation over time scales comparable to the entire cell cycle, or multiples of it, in both mammalian and bacterial cells ^{[45,46,71,74]}. We will also address the treatment of nonlinear noise, highly problematic in the white noise case ^[22], but possible and well-defined for finite and large correlation times. The dynamical derivation presented here recovers in the appropriate limit the fundamental result of ^[33] on noise filtering as a purely geometric construction, and will allow us to identifies the limit of validity of the approach.

Let us then consider the stochastic system:

As before, the dynamics for $x$ are specified by the function $f(x, \lambda)$, and are characterized by the timescale $\tau_x$, while the control parameter $\lambda$ exhibits fluctuations described by Eq (5.5). The variable $\xi = \xi(t)$ is again a Gaussian, zero average, white noise, with correlator given by Eq (2.3). The functions $\mu(\lambda)$ and $\nu(\lambda)$ are kept generic for the moment, but can be chosen so as to reproduce different types of fluctuations, characterized by different distributions, and different correlation times. For instance, they can match the Langevin Eq (4.3) to reproduce log-normal noise.

Instead of studying the neighbourhood of white noise, for $\tau \rightarrow 0$, in this section we are rather interested in developing a formalism to analyze the limit $\tau \gg \tau_x$. The rescaling by the factors $1/\tau$ in Eq (5.5) is introduced so as to slow-down the fluctuations and allows us to carry out a systematic expansion in powers of $1/\sqrt{\tau}$, for $\tau \gg \tau_x$.

The Master Equation of the system specified by Eqs (5.4) and (5.5) results in

where $w_t(x, \lambda)$ is the joint probability density for the $x$ and $\lambda$ variables. The term proportional to $\nu(\lambda)\nu^\prime(\lambda)$ is the so-called Stratonovich drift, corresponding to having assumed the Stratonovich interpretation for the multiplicative noise Eq (5.5).

As it stands, Eq (5.6) is difficult to solve, but it can be solved in stationary conditions when the stochastic fluctuations are slow. By following ^[43,44], we apply the so-called ''switching-curve approximation'', and expand the stationary solution $w_S(x, \lambda)$ of Eq (5.6) in inverse powers of the correlation time $\tau$ of the $\lambda$ process:

We then obtain to the zeroth order in $1/\tau$

where $u(\lambda)$ is defined so that $f(u(\lambda), \lambda) = 0$. By assuming that $w_0(x, \lambda)$ is normalized, $w(\lambda)$ in Eq (5.8) can be identified by solving to order $\tau^{-1}$ the stationary Fokker-Planck Equation associated to Eq (5.5) and supplemented with the Stratonovich interpretation:

The marginalized probability density $p(x)$ for the $x$ process can then be obtained as

By using $\delta(\varphi(x)) = \delta(x-x_0)/|\varphi^\prime(x_0)|$, where $x_0$ is the root of $\varphi(x)$, $\varphi(x_0) = 0$, we readily obtain:

The probability density $w(\lambda)$ can be computed from Eq (5.9), and we obtain

where ${\cal N}$ is a normalization constant.

The mode(s) of the probability density $p(x)$ represent the natural extension of the stable steady state of the corresponding deterministic system, and can be computed explicitly by differentiating $p(x)$. This leads readily to

where $\mu = \mu(u^{-1}(x))$, $\nu = \nu(u^{-1}(x))$, and $u^{\prime \prime} = u^{\prime \prime}(u^{-1}(x))$.

Equation (5.13) is an extension of the standard equation identifying the modes of a multiplicative stochastic process supplemented with the Stratonovich prescription, Eq (2.8). The extra term $-D \nu^2 \left(d u^{-1}(x)/dx \right) u^{\prime \prime}$ in Eq (5.13) accounts for dynamical modifications of the deterministic system which are induced by the slow fluctuations. This term is identically zero only in the case when the the relationship between the variable $\lambda$ and the variable $x$ is linear (so that $u^{\prime \prime} = 0$), while in all other cases nontrivial modifications of the dynamics can be expected, further and different from those arising in the white noise limit because of the Stratonovich drift. Again, whether these modifications will correspond to a shift of parameter values, a change of stability properties of the deterministic attractor states, or a change in their number (namely a purely noise-induced transition) will need to be assessed case by case, once the dynamics of the $x$ and $\lambda$ processes are specified.

As expected, Eq (5.11) is equivalent Eq (5.2). Our derivation highlights that the noise filtering approach is in fact valid in the case when no finite-time internal dynamics is present for the $x$ variable, namely in the case of infinite timescale separation $\tau/\tau_x \rightarrow 0$, and $x$ adapts istantaneously to the value dictated upon it by the dynamics of $\lambda$. This is the essence of the geometric construction in ^[33].

It should be noted that the approach here presented has the advantage that it can be applied also in those cases when the driving distribution $w(\lambda)$ is not known. Equation (5.13) allows us to identify the modes of the probability distribution $p(x)$ starting directly from the Langevin description of the process $\lambda$, Eq (5.5), without requiring the analytical solution of its stationary probability density, Eq (5.12).

6.   The repressed gene: slow log-normal fluctuations on degradation

Shahrezaei et al. ^[62] consider log-normal fluctuations on the degradation rate in a simple model of gene expression, and show that no noise-induced transition takes place in this case. Let us examine this case again by using our model system of the repressed gene.

For $\tau\gg\tau_x$ we can apply the results derived in the previous section and explore the effect of slow log-normal noise on the degradation rate $k$. By formally setting $\lambda \equiv k$, the funtions $\mu$ and $\nu$ entering Eq (5.5) are thus identified as

with the function $u(k)$ resulting trivially in $u(k) = \beta/k$. Equation (5.12) can be easily integrated, and results in:

The resulting $p(x)$ is nontrivial in this case, as the original log-normal noise is transformed nonlinearly into the non-Gaussian probability density (6.3). Despite this, the extrema are given by

which is single-valued for all parameter values.

Simulations confirm nicely these analytical predictions, as shown in Figure 4. In particular the single solution given by Eq (6.4) implies that no transition is possible in this system, despite our choice of log-normal noise to model fluctuations on $k$. It should be noted that the same result reported in ^[62] is derived by adopting the UCNA approach, and therefore valid for the correlation time of the noise being either very small or very large.

7.   The repressed gene: fluctuations on the repressor

Let us now come to the core of the nonlinear noise filtering approach ^[33]. By using again our model system, the repressed gene, in this section we want to highlight with concrete examples the two main factors of the nonlinear noise filtering approach responsible for the occurrence of noise-induced transitions: The nonlinear dependency of the gene dynamics on fluctuating parameters, and the input distribution of the fluctuations.

It is worth noticing that already in this simple system, different choices can be made for the identification of the nonlinear relevant parameter of interest, which correspond to different biological interpretations. For instance, Ochab-Marcinek et al.^[33] focus on fluctuations in the relevant transcription factor, in our case the repressor $R$, while Birtwistle et al ^[35] and Dobrzyǹski et al. ^[37] focus instead on fluctuations on activation thresholds. Both choices are legitimate, but correspond clearly to different biochemical processes. An important question here is how to justify the assumed distribution of fluctuations on the relevant parameters. We will discuss this specific point in our Concluding Remarks, Section 8.

By following ^[33] we make the choice of considering slow fluctuations acting on $R$. First we will consider Gaussian noise, and then we will address log-normal noise.

In the case of Gaussian fluctuations acting on the parameter $R$, we assume $R \rightarrow \bar{R} + \eta(t)$, where $\eta(t)$ is the standard Ornstein-Uhlenbeck process, whose dynamics is given by Eq (4.2). Equation (5.12) can be readily integrated once the function $u(R)$ is defined as

By direct integration of Eq (5.12), we readily obtain

again showing a nonlinear tranformation from the original Gaussian noise.

Simulations agree again very well with the analytical predictions, as it it is shown in Figure 5. The timescale of the fluctuations $\tau$ is crucial for matching our theoretical predictions with the numerical results. Given $k = 10^{-4}$ s$^{-1}$, we set $\tau = 10^5$ s, to capture the dynamics of slow fluctuations of $R$.

However, only one positive mode is present, which extends the deterministic solution. In fact, for $x \neq 0$, Eq (5.13) becomes

which for $D \neq 0$ accounts for the two real solutions:

The solution of Eq (7.3) for $D = 0$ coincides trivially with the solution of the deterministic case, with $x_1$ being the stochastic continuation of it (Figure 6). The $x_2$ branch emerges instead for $D > 0$, signalling that no transition is taking place in the system at finite $D$ (the only critical point being the trivial one, $D = D_c = 0$). Despite the appearance of a second mode, this is to be biologically discarded because at negative $x$ values.

Let us consider now log-normal noise. By direct integration of Eq (5.12), we readily obtain

and, for $x\neq 0$ and $x \neq g/k$, from Eq (5.13)

Simulations show an excellent agreement with the theoretical predictions, as shown in Figure 7. Further to the excellent matching of the full distribution for $R$, the predictions for the modes as given by Eq (7.6) is also verified. In Figure 8 we show the extrema of the stationary probability Eq (7.5) as function of $\bar{R}$ for different noise intensities, obtained by the numerical solution of Eq (7.6) with parameter values as in Figure 7. The curve associated to $D = 0$ (deterministic case) is single valued for all values of $\bar{R}$, with further increase of $D$ provoking the appearance of a second stable mode for intermediate $\bar{R}$ values.

These results confirm that the emergence of noise-induced transitions in simple regulatory dynamics depends strongly on the type of noise and on the nonlinearity involved. From a biological perspective it is then essential to assess carefully noise distributions and network reconstructions to be able to make clear predictions. Certainly non-Gaussian noises, in the case explored here log-normal, long fluctuations correlations, and consideration of nonlinear parameters, seem to be essential ingredients for the emergence of noise-induced transitions.

8.   Concluding remarks

In this paper we have summarized how the concept of noise-induced transitions has evolved over the last few decades, and has now become an active area of research in gene systems and molecular cellular processes in general.

With the help of a simple gene expression system, the repressed gene, we have reviewed different approaches to model extrinsic fluctuations, ranging from assuming Gaussian white and coloured extrinsic noise, to the adoption of bounded noises ^[52] and nonlinear noise filtering methodologies ^[33]. The repressed gene model offers scope for interesting analysis as it is defined in terms of parameters that appear in the dynamical equations both linearly and nonlinearly, a fact that is now known to contribute differently to transitions to bimodality ^[37]. In this paper we have reviewed the behaviour of this systems with respect to its linear dependency on the degradation rate $k$, and the nonlinear dependency of the repressor concentration, $R$.

When the noise appears linearly in the repressed gene, namely on the degradation rate, standard approaches can be adopted to carry out analytical predictions in the case of Gaussian white noise. These standard tools, extended to coloured noise by the powerful UCNA approach ^[47], and adopted in a variety of different gene regulatory systems ^[49,62], do not predict noise-induced transitions. The same result holds true for log-normal noise characterized by long correlation time ^[62].

However, when fluctuations appear in a nonlinear fashion, noise-induced transitions become possible. From a mathematical perspective, treating nonlinear noise dynamically is usually a hard problem, but recent progress has made it possible by adopting the noise filtering approach introduced in ^[33,34]. We have here reviewed this approach, and proposed a dynamical derivation of it in the limit of large correlation time of the noise. The key equation describing the approach (Eq (5.2) or Eq (5.11)) highlights how nontrivial behaviours, including transitions to bimodality, depend on the combined effect of two dynamical features. The first one is represented by the derivative term in Eqs (5.2) and (5.11), which represents the response function of the system, whose nonlinearity is responsible for the distortion of the input probability distribution. For linear filtering, this derivative is just a numerical factor, and hence ininfluent in modifying dynamics between input and output. The second important ingredient is of course the nature of the input distribution itself. We show here with a concrete example based on the repressed gene that if the input distribution is Gaussian no biologically meaningful bimodality is obtained (a second mode indeed appears, but it is at negative concentrations), while if it is log-normal a physical transition occurs, with two positive mode emerging. The same analysis holds for Gamma distributed noise, as presented in ^[33,35,36].

Whether the Gaussian noise approximation for extrinsic noise (either white or coloured) is biologically feasible is a delicate question. While it is indeed true that the Gaussian distribution produces negative values for fluctuating parameters, it is unquestionable that it has got its own merits, which also rely on the development of poweful and thoroughly tested analytical results. In fact in recent papers, such as ^[62] and ^[49] the Gaussian approximation has been carried out, leading to clear and reliable predictions on the behaviour of the system. It is fair to say that according to the the system itself, and the type of noise one wants to model, the Gaussian approximation may still be an attractive and valid alternative. Current experimental results however, support the relevance of bounded noises to fit models to data. For instance, the already mentioned studies ^[35,36,37] show an excellent matching of models to experiments by adopting Gamma noise and the nonlinear filtering approach, which allows for the identification of the related noise-induced transitions.

A number of challenging problems remains open. One important one concerns the interplay between intrinsic and extrinsic noise. Intrinsic and extrinsic noises are invitably present in the system, and ideally should be both considered in the models proposed. Several attempt have been made to tackle this issue, with promising results. For instance, in ^[75] the chemical Lengevin equation with the Itô prescription has been adopted to model intrinsic fluctuations, while extrinsic noise is treated with the Stratonovich prescription. A similar approach has been put forward in ^[49] for the self-regulating gene, and in ^[76] in the context of reducing model dimensionality of feedback loops, by adopting the so-called loop opening approach defined in ^[77,78]. It is worth noticing that in all these cases the SDEs obtained are 'doubly' multiplicative, in that both noise sources provoke the apperance of multiplicative noise terms. Since intrinsic fluctuations (white) and dynamical variables are uncorrelated, for these the Itô prescription follows, while the fact that extrinsic fluctuations (coloured) develop time correlations with the dynamical variables is reflected in the adoption of the Stratonovich prescription. Because of these rich dynamics, noise-induced transitions are likely, and indeed emerge in both studies ^[49,75].

A second open problem concerns the origin of the different input noise distributions chosen to model experimental data and noise-induced transitions. A range of different assumptions has been made. For instance, fully bounded distributions have been postulated to describe tumor biology ^[54,55], while semibounded noises (Gamma and log-normal distributed) have been postulated in the studies ^[33,35,37] either on activation thresholds, or on regulating transcription factors. In the case of a truly external parameter, like an activation threshold, this choice is made on a phenomenological basis, and it may be hard to justify it in terms of molecular processes. However, fluctuations affecting a transcription factor deserve more attention. These may in fact have an intrinsic origin, as they come about through the dynamical processes involved in gene expression, as shown in ^[69,70], where it is predicted that fluctuations in gene expression follow a Gamma distribution, a result confirmed experimentally in ^[56]. In the approach adopted in ^[33], and followed by ourselves in this review, it is crucial to treat the transcription factor as an external parameter in the dynamics of the target gene, so that adopting a stationary Gamma (or log-normal) distribution to describe its fluctuations is fully justified. However, a full dynamical treatment would require modelling in terms of a 3-dimensional system, accounting for transcription factor production, target gene transcription factor-mediated regulation, and noise dynamics. To reduce this 3-dimensional system to a 2-dimensional system (target gene dynamics with the transcription factor appearing as an external parameter, plus noise), as ^[33] and we do, may be a delicate step. While achieving this reduction is mathematically manageable for deterministic systems, how to realize it consistently for a stochastic system is more challenging. This requires eliminating fast variables by stochastic adiabatic elimination ^[79,80] or projection techniques, a possibility currently under investigation ^[81,82]. In this sense, intrinsic noise could in fact behave as extrinsic and be ultimately responsible itself for nontrivial behaviours ^[83], with the boundary between intrinsic and extrinsic dynamics becoming fuzzy and dependent on the definition of networks and subnetworks.

A third problem concerns the nature and the implications of the observed bimodality. Even though the population can be distributed in two or more subpopulations in terms of gene profiles, the internal dynamics of switching between them will be relevant for drawing biological conclusions. The same bimodal distribution can be realized in quasi-static terms, with limited transitions between states, or very fast transitions between them, or anything in between. The timescale of the 'hopping' process will be relevant to assign observable phenotypes to those gene configurations corresponding to specific states. Only slow hopping corresponds to what is commonly referred to as a clear switching behaviour, or, mathematically speaking to (possibly weak) ergodicity breaking ^[41]. Hysteresis, as a signature of true bimodal dynamics, can be seen in the same way. P{á}jaro et al. ^[84] point out that indeed hysteresis may only have a transient nature under the presence of intrinsic noise, as eventually the system converges to equilibrium, and restores ergodic behaviour. This phenomenon is particularly relevant for slow fluctuations, when the rate of convergence to equilibrium is itself slow ^[85], and needs to be studied case by case.

All these results indicate that particular care should be taken when reconstructing gene regulatory networks under the evidence of bimodal distributions of gene expression levels. The usual implication that these require a wiring diagram including feedback loops may be false, as the same bimodality may be produced by nonlinear extrinsic noise. It is intriguing to hypothesize that management of noise could provide an alternative way evolution adopts to realize bimodal physiology of cellular systems ^[14,86]. Certainly the mechanisms for noise-induced transitions here reviewed are promising and worthy of further analysis and experimental validation. These will contribute to our understanding of the fundamental relevance of noise in biological systems.

Acknowledgements

This work was supported by the UK Biotechnology and Biological Sciences Research Council (BBSRC) under grant BB/L007789/1. AR wishes to thank the Editors of this Special Issue for the invitation to submit this paper.

Conflict of interest

The authors declare no conflict of interest in this paper.