The Interplay of Intrinsic and Extrinsic Bounded Noises in Biomolecular Networks

After being considered as a nuisance to be filtered out, it became recently clear that biochemical noise plays a complex role, often fully functional, for a biomolecular network. The influence of intrinsic and extrinsic noises on biomolecular networks has intensively been investigated in last ten years, though contributions on the co-presence of both are sparse. Extrinsic noise is usually modeled as an unbounded white or colored gaussian stochastic process, even though realistic stochastic perturbations are clearly bounded. In this paper we consider Gillespie-like stochastic models of nonlinear networks, i.e. the intrinsic noise, where the model jump rates are affected by colored bounded extrinsic noises synthesized by a suitable biochemical state-dependent Langevin system. These systems are described by a master equation, and a simulation algorithm to analyze them is derived. This new modeling paradigm should enlarge the class of systems amenable at modeling. We investigated the influence of both amplitude and autocorrelation time of a extrinsic Sine-Wiener noise on: the Michaelis-Menten approximation of noisy enzymatic reactions, which we show to be applicable also in co-presence of both intrinsic and extrinsic noise, a model of enzymatic futile cycle and a genetic toggle switch. In and we show that the presence of a bounded extrinsic noise induces qualitative modifications in the probability densities of the involved chemicals, where new modes emerge, thus suggesting the possible functional role of bounded noises.


Author Summary
Realistic modeling the dynamics of genetic networks is a non-trivial task that requires choosing a suitable level of abstraction.For example, within cells the molecules interacting in a network can be present in either small or large numbers.In the former case their time course is typically characterized by wild random oscillations closely mimicking the randomness of chemical reactions.In the latter, instead, these fluctuations are invisible, due to the "law of large numbers".So in this case the molecular concentrations should theoretically stabilize to nice smooth steady states.
However, the presence of perturbing external factors may induce noisy fluctuations more or less disrupting the theoretical "nice behavior".Surprisingly, this disruption may be constructive: a singlefunction network may gain additional biological functionalities in presence of perturbations.In real world, these two kinds of randomness are not separated: proteins of a specific genetic network can be present in small number and be perturbed by external noises.This is actually a topic only sparsely explored in systems biology.
A factor that makes complex the study of these phenomena is that external disturbances are classically represented through unbounded Gaussian noises, which are actually unrealistic and may induce nonbiological artifacts.Thus we focus on studying systems with both small number of molecules and external bounded noises.
A second key concept is that the dynamic behavior of a network is never totally deterministic, but it exhibits more or less strong stochastic fluctuations due to its interplay with many, and mainly unknown, other networks, as well as with various random signals coming from the extracellular world.For long time the stochastic effects due these two classes of interactions were interpreted as a disturbance inducing undesired jumps between states or, with marginally functional role, as an external initial input directing towards one of the possible final states of the network in study.In any case, in the important scenario of deterministically monostable networks the stochastic behavior under the action of extrinsic noises was seen as monomodal.In other words, external stochastic effects were seen similarly as in radiophysics, namely as a disturbance more or less obfuscating the real signal, to be controlled by those pathways working as a low-pass analog filter [15,16].For these reasons, a number of theoretical and experimental investigations focused on the existence of noise-reducing sub-networks [15,17,18].However, it has been recently shown the existence of fundamental limits on filtering noise [19].
Moreover, if noises were only pure nuisances, there would be an interesting consequence.Indeed, in such a case a monostable network in presence of noise should exhibit more or less large fluctuations around the unique deterministic equilibrium.In probabilistic languages this means that the probability distribution of the total signal (noise plus deterministic signal) should be a sort of "bell" centered more or less at the deterministic equilibrium, i.e. the probability distribution should be "unimodal".However, at the end of seventies it became clear in statistical physics that the real stochastic scenario is far more complex, and the above-outlined correspondence between deterministic monostability and stochastic monomodality in presence of external noise was seriously challenged [20].Indeed, it was shown that many systems that are monostable in absence of external stochastic noises have, in presence of random Gaussian disturbances, multimodal equilibrium probability densities.This counter-intuitive phenomenon was termed noise-induced transition [20], and it has been shown relevant also in genetic networks [21,22].
We above focused mainly on external random perturbations acting on genetic networks.In the meantime, experimental studies revealed the other and equally important role of stochastic effects in biochemical networks by showing that many important transcription factors, as well as other proteins and mRNA, are present in cells with very low concentrations, i.e. with a small number of molecules [23][24][25].Moreover, it was shown that RNA production is not continuous, but instead it has the characteristics of stochastic bursts [26].Thus, a number of investigations has focused on this internal stochastic effect, the "intrinsic noise" as some authors term it [28,29].In particular, it was shown -both theoretically and experimentally -that also the intrinsic noise may induce multimodality in the discrete probability distribution of proteins [22,30].However, the fact that intrinsically stochastic systems may exhibit behaviors similar to systems affected by extrinsic gaussian noises was very well known in statistical and chemical phsyics, where this was theoretically demonstrated by approximating the exact Chemical Master Equations with an appropriate Fokker-Planck equation [31][32][33], an approach leading to the Chemical Langevin Equation [34].
Thus, after that for some time noise was mostly seen as a nuisance, more recently it has finally been appreciated that the above-mentioned and other noise-related phenomena may in many cases have a constructive, functional role (e.g.see [35,36] and references therein).Indeed, for example, noise-induced multimodality allows a transcription network for reaching states that would not be accessible if the noise was absent [22,35,36].Phenotype variability in cellular populations is probably the most important macroscopic effect of intracellular noise-induced multimodality [35].
In Systems Biology, from the modeling point of view Swain and coworkers [24] were among the first to study the co-presence of both intrinsic and extrinsic randomness, by stressing the synergic role in modifying the velocity and average in the context of the basic network for the production and consumption of a single protein, in absence of feedbacks.These and other important effects were shown, although nonlinear phenomena such as multimodality were absent.The above study is also remarkable since: (i) it has stressed the role of the autocorrelation time of the external noise and, differently from other investigations, (ii) it has stressed that modeling the external noise by means of a Gaussian noise, either white or colored, may induce artifacts.In fact, since the perturbed parameters may become negative, the authors employed a lognormal positive noise to model the extrinsic perturbations.In particular, in [24] a noise obtained by exponentiating the classical Orenstin-Uhlenbeck noise was used [20].
From the data analysis point of view, You and collaborators [37] and Hilfinger and Paulsson [38] recently proposed interesting methodologies to infer by convolution the contributions of extrinsic noise also in some nonlinear networks, including a synthetic toggle switch [37].
Our general aim here is manifold.Indeed, we want to start by investigating the co-presence of both extrinsic and intrinsic randomness in nonlinear genetic networks in the important case where the external perturbations are not only non-Gaussian, but also bounded.Indeed, by imposing a bounded extrinsic noise we increase the degree of realism of a model, since the external perturbations must not only preserve the positiveness of reaction rates, but must also be bounded.Moreover, it has also been shown in other contexts such as oncology and statistical physics that: (i) bounded noises deeply impact on the transitions from unimodal to multimodal probability distribution of state variables [39][40][41][42][43] and (ii) the dynamics of a system under bounded noise may be substantially different from the one of systems perturbed by other kinds of noises, for example there is dependence of the behavior on the initial conditions [40,42].
This approach opens a number of questions, two of which we shall try to assess here.The first question is to identify a suitable mathematical framework to represent mass-action biochemical networks perturbed by bounded noises (or simply left-bounded), which in turn can depend on the state of the system.To this extent, in the first part of this work we derive a master equation for these kinds of systems in terms of the differential Chapman-Kolgomorov equation (DCKE) [31,44] and propose a combination of the Gillespie's Stochastic Simulation Algorithm (SSA) [27,28] with a state-dependent Langevin system, affecting the model jump rates, to simulate these systems.
The second question relates to the possibility of extending, in this "doubly stochastic" context, the Michaelis-Menten Quasi Steady State approximation (QSSA) for enzymatic reactions [45].We face the validity of the QSSA in copresence of both types of noise in the second part of this work, where we numerically investigate the classical Enzyme-Substrate-Product network.The application of QSSA in this network has been recently investigated by Gillespie and coworkers when extrinsic noise is absent [46].Based on our results, we propose the extension of the above structure also to more general networks than those ruled by the rigourous mass-action law via a stochastic QSSA.
Finally, in the third part we investigate the interplay between intrinsic randomness and extrinsic bounded noises in two cases of interest in biology: (i) a futile cycle [22] and (ii) a genetic toggle switch [7], which is a fundamental motif for cellular differentiation and for other switching functions.As expected, the co-presence of both intrinsic stochasticity and bounded extrinsic random perturbations suggests the presence of possibly unknown functional roles for noise in both networks.The described noise-induced phenomena are shown to be strongly related to physical characteristics of the extrinsic noise such as the noise amplitude and its autocorrelation time.

Noise-free stochastic chemically reacting systems
We start by recalling the Chemical Master Equation and the Stochastic Simulation Algorithm (SSA) by Doob and Gillespie [27,28].Systems where the jump rates are time-constant are hereby referred to as stochastic noise-free systems.We consider a well stirred system of molecules belonging to N chemical species {S 1 , . . ., S N } interacting through M chemical reactions R 1 , . . ., R M .We represent the (discrete) state of the target system with a N -dimensional integer-valued vector X(t) = (X 1 (t), . . ., X N (t)) where X i (t) is the number of molecules of species S i at time t.To each reaction R j is associated its stoichiometric vector ν j = (ν 1,j , . . ., ν N,j ), where ν i,j is the change in the X i due to one R j reaction.The stoichiometric vectors form the N × M stoichiometry matrix D = ν 1 ν 2 . . .ν M .Thus, given X(t) = x the firing of reaction R j yields the new state x + ν j .A propensity function a j (x) [27,28] is associated to each R j so that a j (x)dt, given X(t) = x, is the probability of reaction R j to fire in state x in the infinitesimal interval [t, t + dt).Table 1 summarizes the analytical form of such functions [27].For more generic form of the propensity functions (e.g.Michaelis-Menten, Hill kinetics) we refer to [49].
We recall the definition of the Chemical Master Equation (CME) [27,28,47,48] describing the timeevolution of the probability of a system to occupy each one of a set of states.We study the time-evolution of X(t), assuming that the system was initially in some state x 0 at time t 0 , i.e.X(t 0 ) = x 0 .We denote with P(x, t | x 0 , t 0 ) ≡ P(x, t | ω) the probability that, given X(t 0 ) = x 0 , at time t it is X(t) = x.From the usual hypothesis that at most one reaction fires in the infinitesimal interval [t, t + dt), it follows that the time-evolution of P(x, t | ω) is given by the following partial differential equation termed "master equation" The CME is a special case of the more general Kolmogorov Equations [50], i.e. the differential equations corresponding to the time-evolution of stochastic Markov jump processes.As it is well known, the CME can be solved analytcally only for a very few simple systems, and normalization techniques are sometimes adopted to provide approximate solutions [51].However, exact algorithmic realization of the process associated to the CME are possible by using the Doob-Gillespie Stochastic Simulation Algorithm (SSA) [27,28,47,48], summarized as Algorithm 1.Although equivalent formulations exist [27,28,52], as well as some approximations [49,53,54], here we consider its Direct Method formulation without loss of generality.
The SSA is an exact dynamic Monte-Carlo method describing a statistically correct trajectory of a discrete non-linear Markov process, whose probability density function is the solution of equation ( 1) [55].The SSA computes a single realization of the process X(t), starting from state x 0 at time t 0 and up to time T .Given X(t) = x the putative time τ for the next reaction to fire is chosen by sampling an exponentially distributed random variable, i.e. τ ∼ Exp(a 0 (x)) where a 0 (x) = M j=1 a j (x) and ∼ denotes Algorithm 1 SSA (t 0 , x 0 , T ) 1: set x ← x 0 and t ← t 0 ; 2: while t < T do 3: τ ← a 0 (x) −1 ln(r 1 −1 ); set x ← x + ν j and t ← t + τ ; 8: end while the equality in law between random variables.The reaction to fire R j is chosen with weighted probability a j (x)/a 0 (x), and the system state is updated accordingly.
The correctness of the SSA comes from the relation between the jump process and the CME [27,55].In fact, the probability, given X(t) = x, that the next reaction in the system occurs in the infinitesimal time interval [t + τ, t + τ + dτ ), denoted P(τ | x, t), follows (2) since P(τ, j | x, t) = a j (x) exp (−a 0 (x)τ ) is the probability distribution of the putative time for the next firing of R j , and the formula follows by the independency of the reaction firings.Notice that in equation (2) a 0 (t ′ ) represents the propensity functions evaluated in the system state at time t ′ > t, i.e. as if they were time-dependent functions.In the case of noise-free systems that term evaluates as a 0 (x) for any t ∈ [t, t + τ ], i.e. it is indeed time-homogenous whereas in more general cases it may not, as we shall discuss later.Finally, the probability of the reaction to fire at t + τ to be R j follows by conditioning on j, that is

Noisy stochastic chemically reacting systems
We now introduce a theory of stochastic chemically reacting systems with (un)bounded noise in the jump rates by combining Stochastic Differential Equations and the SSA.Here we consider a system where each propensity function may be affected by a extrinsic noise term.In general, such a term can be either a time or state-dependent function, and the propensity function for reaction R j reads now as where a j (x) is a propensity function of a type listed in Table 1.The noisy perturbation term L * j (t) is positive and bounded by some so we are actually considering both bounded and right-unbounded noises, i.e.C j = +∞.In the former case we say that the j-th extrinsic noise is bounded, in the latter that it is left-bounded.Note that in applications we shall mainly consider unitary mean perturbations, that is We consider here that the extrinsic noisy disturbance L * j (t) is a function of a more generic Σ-dimensional noise ξ(t) with 1 ≤ Σ ≤ M so we write L * j (t) = L j (ξ(t)) and equation ( 4) reads as a j (x, t) = a j (x)L j (ξ(t)) .
Notice that the use of a vector in equation ( 6) provides the important case of multiple reactions sharing the same noise term, i.e. the reactions may be affected in the same way by a unique noise source.In equation ( 6) L j is a continuous functions L j : R Σ → R and ξ(t) ∈ R Σ is a colored and, in general, non-gaussian noise that may depend on the state X(t) of the chemical system.The dynamics of ξ(t) is described by a Σ-dimensional Langevin system Here, η is a Σ-dimensional vector of uncorrelated white noises of unitary intensities, g is a Σ × Σ matrix which we shall mainly consider the be diagonal and f, g h,k : R Σ × R N → R Σ .When ξ(t) does not directly depend on X(t), i.e. the extrinsic noise depends on an external source, which is the kind of noise we mainly consider, equation ( 7) reduces to We stress that the "complete" Langevin system in equation ( 7) is not a mere analytical exercise, but it has the aim of phenomenologically modeling extrinsic noises that are not totally independent of the process in study.

The Chapman-Kolmogorov Forward Equation
When a discrete-state jump process as one of those described in previous section is linked with a continuous noise the state of the stochastic process is the vector and the state space of the process is now N N × R Σ .Our total process can be considered as a particular case of the general Markov process where diffusion, drift and discrete finite jumps are all co-present for all state variables [31,44].For this very general family of stochastic processes the dynamics of the probability of being in some state z at time t, given an initial state z 0 at time t 0 shortly denoted as ω, is described by the differential Chapman-Kolgomorov equation (DCKE) [31,44], whose generic form is Here A j forms the drift vector for z, B i,j (z, t) the diffusion matrix and W the jump probability.For an elegant derivation of the DCKE from the integral Chapman-Kolgomorov equation [50] we refer to [44].This equation describes various systems, in fact we remind that (i) the Fokker-Planck equation is a particular case of the DCKE without jumps (i.e.W (z | h, t) = 0), (ii) the CME in equation ( 1) is the DCKE without brownian motion and drift (i.e.A(z, t) = 0 and B(z, t) = 0), (iii) the Liouville equation is the DCKE without brownian motion and jumps (i.e.A(z, t) = 0 and W (z | h, t) = 0) and (iv) the ODE with jumps correspond to the case where only diffusion is absent (i.e.B(z, t) = 0).We stress that, at the best of our knowledge, this is the first time where a master equation for stochastic chemically reacting systems combined with bounded noises is considered.Let be the probability that at time t it is X(t) = x and ξ(t) = ξ, given X(t 0 ) = x 0 and ξ(t 0 ) = ξ 0 .The time-evolution of P(z, t | ω) is equation (10) where drift and diffusion are given by the Langevin equation (7), that is with × the standard vector multiplication and g T the transpose of g.Moreover, since only finite jumps are possible, then the jump functions and diffusion satisfy for any i, j = 1, . . ., N , and noise ξ * ∈ R Σ .Summarizing, for the systems we consider the DCKE in equation ( 10) reads as This equation is the natural generalization of the CME in equation ( 1), and completely characterize noisy systems.As such, however, its realization can be prohibitively difficult and is hence convenient to define algorithms to perform the simulation of noisy systems.

The SSA with Bounded Noise
We now define the Stochastic Simulation Algorithm with Langevin Noise (SSAn).The algorithm performs a realization of the stochastic process underlying the system where a (generic) realization of the noise is assumed.As for the CME and the SSA, this corresponds to computing a realization of a process satisfying equation ( 14).The SSAn takes inspiration from the (generic) SSA with time-dependent propensity functions [56] as well as the SSA for hybrid deterministic/stochastic systems [57][58][59][60], thus generalizing the jump equation ( 2) to a time inhomogeneous distribution, which we discuss in the following.
For a system with M reactions the time evolution equation for X(t) is where {N j (t) | t ≥ t 0 } is the stochastic process counting the number of times that R j occurs in [t 0 , t] with initial condition N j (t 0 ) = 0.For Markov processes N j (t) is an inhomogeneous Poisson process satisfying when X(t) = x.In hybrid systems this is is a doubly stochastic Poisson process with time-dependent intensity, in our case this is a Cox process [61,62] since the intensity itself is a stochastic process, i.e. it depends on the stochastic noise.More simply, in noise-free systems, this equation evaluates as a j (x)dt, thus denoting a time homogeneous Poisson process.As in [57,59,60,63,64] such a process ca be transformed in a time homogenous Poisson process with parameter 1, and a simulation algorithm can be exploited.Let us denote with T j (t) the time at next occurrence of reaction R j after time t, then follows by equation ( 16) and higher order terms vanish by the usual hypothesis that the reaction firings are locally independent, as in the derivation of equation ( 1).Given the system to be in state x at time t, the transformation Algorithm 2 SSAn (t 0 , x 0 , T ) find τ by solving set x ← x + ν j and t ← t + τ ; 7: end while which is a monotonic (increasing) function of τ is used to determine the putative time for R j to fire.Given a sequence r j,k of independent exponential random variables with mean 1 for j = 1, . . ., M and k ∈ N, equation ( 16) implies that This provides that, if the systems is in state X(t) = x, then the next time for the next reaction firing of R j is the smallest time τ > 0 such that with r ∼ Exp(1), and thus the next jump of the overall system is taken as the minimum among all possible times, that is by solving equality with r ∼ Exp(1).This holds because min{T j | j = 1, . . ., M } is still exponential with parameter M j=1 A j (t + τ, t) and the jumps are independent.We remark that for a noise-free reaction A j (t + τ, t) = τ a j (x), thus suggesting that the combination of noisy and noise-free reactions is straightforward.The index of the reaction to fire is instead a random variable following The SSAn is Algorithm 2.
Step 4 is the (parallel) solution of both equation (21) and Langevin system (7), step 5 samples values for j according to equation (22).As far as step 4 is concerned, it is worth nothing that given X(t) = x for any τ the Langevin equation ( 7) depends only on ξ(t) and the constant x.To this extent, a single trajectory of the vectorial noise process in [t, t + τ ] is This is a discretization of a continuous noise, thus inducing an approximation, but is in general the only possible approach.To reduce approximation errors the maximum size of the jump in the noise realization, i.e. the noise granularity ∆ s = t s+1 − t s , should be much smaller than the minimum autocorrelation time of the perturbing stochastic processes L j (ξ(t)).Finally, the integral in equation (21), evaluated in step 4, is a conventional Lebesgue integral since the perturbation L j (ξ(t)) is a colored stochastic process [65].As an example, given ξ t,τ by a linear interpolation scheme it holds where L s j = L j (ξ(t s )) for t s ∈ ξ t,τ and ∆ s the noise granularity.

Extension to non mass-action nonlinear kinetic laws
Large networks with large chemical concentrations, i.e. characterized by deterministic behaviors, are amenable to significant simplifications by means of the well known Quasi Steady State Approximation (QSSA) [5,45,46,66].The validity conditions underlying these assumptions are very well-known in the context of deterministic models [45] , despite not much being known for the corresponding stochastic models.Recently, Gillespie and coworkers [46] showed that, in the classical Michaelis-Menten Enzyme-Substrate-Product network, a kind of Stochastic QSSA (SQSSA) may be applied as well, and that in such its limitations are identical to the deterministic QSSA.Thus, it is of interest to consider SQSSAs also in our "doubly stochastic" setting, even though possible pitfalls may arise due to the presence of the extrinsic noises.As an example, in Section we will present the results of the numerical experiments similar to those of [46], with the purpose of validating the SQSSA for noisy Michaelis-Menten enzymatic reactions.
Of course, in a SQSSA not only the propensities may be nonlinear function of state variables, but they may depend nonlinearly also on the perturbations, so that instead of the elementary perturbed propensities we shall have generalized perturbed propensities of the form where ψ is a vector with elements ψ j = L j (ξ) for j = 1, . . ., M .This makes possible, within the above outlined limitation for the applicability of the SQSSA, to write a DCKE for these systems as As far as the simulation algorithm is concerned, it remains quite close to Algorithm 2 provided that the jump times are sampled according to the following distribution

Results
We performed SSAn-based analysis of some simple biological networks, actually present in most complex realistic networks.We start by studying the legitimacy of the stochastic Michaelis-Menten approximation of when noise affects enzyme kinetics [46].Then we study the role of the copresence of intrinsic and etrinsic bounded noises in a in a model of enzymatic futile cycle [22] and, finally, in a bistable "toggle switch" model of gene expression [13,73].All the simulations have been performed by a Java implementation of the SSAn (freely downloadable at http://sites.google.com/site/giuliocaravagna/)running on a cluster of 15 dual-core nodes with 2.0 Ghz processor and 1 GB of memory.
The Sine-Wiener noise [39].The bounded noise µ(t) that we use in our simulations is obtained by applying a bounded continuous function h : R → R to a random walk W , i.e.W ′ (t) = η(t) with η(t) a white noise1 .We have The effect of the truncation of the tails induced by the approach here illustrated is that, due to this "compression", the stationary probability densities of this class of processes satisfy Probably the best studied bounded stochastic process obtained by using this approach is the so-called Sine-Wiener noise [39], that is where β is the noise intensity and τ is the autocorrelation time.The average and the variance of this noise are: and its autocorrelation is such that [39] Note that, since we mean to use noises of the form 1 + µ(t), i.e. the unitary-mean perturbations in equation ( 6), then the noise amplitude must be such that 0 ≤ β ≤ 1.
For this noise, the probability density is the following [77]: By these properties, this noise can be considered a realistic extension of the well-known symmetric dichotomous Markov noise [67] a(t), whose stationary density is 1 2 (δ(a − β) + δ(a + β)), for a ∈ R and δ the Dirac delta function.

Enzyme kinetics
Enzyme-catalyzed reactions are fundamental for life, and in deterministic chemical kinetics theories are often conveniently represented in an approximated non mass-action form, the well-known Michaelis-Menten kinetics [5,45,46].Such approximation of the exact mass-action model is based on a Quasi Steady-State Assumption (QSSA) [45,66], valid under some well known conditions.In [46] it is studied the legitimacy of the Michaelis-Menten approximation of the Enzyme-Substrate-Product stochastic reaction kinetics.Most important, it is shown that such a stochastic approximation, i.e. the SQSSA previously discussed, obeys the same validity conditions for the deterministic regime.This suggests the legitimacy of using -in case of low number of molecules -the Gillespie algorithm not only for simulationg massaction law kinetics, but more in general to simulate more complex rate laws, once a simple conversion of Table 2. Exact model (left) and Michaelis-Menten approximation (right) of enzymatic reactions: the stoichiometry matrixes (rows in order E, S, ES, P ) and the propensity functions.Figure 1.Noise-free Enzyme-Substrate-Product system.Averages of 1000 simulations for P , plotted with its standard deviation (dotted), for both exact and approximated Michaelis-Menten models.We have set c 1 = c 3 = 1 and c 2 = 10 and the initial configuration is (i) in left (E, S, ES, P ) = (1, 10, 0, 0) and (ii) in right (E, S, ES, P ) = (100, 10, 0, 0).deterministic Michaelis-Menten models is performed and provided -of course -that the SQSSA validity conditions are fulfilled.
In this section we investigate numerically whether the Michaelis-Menten approximations and the stochastic results obtained in [46] still hold true in case that a bounded stochastic noise perturb the kinetic constants of the propensities of the exact mass-action law system Enzyme-Substrate-Product.Let E be an enzyme, S a substrate and P a product, the exact mass-action model of enzymatic reactions comprises the following three reactions where c 1 , c 2 and c 3 are the kinetic constants.The network describes the transformation of substrate S into product P , as driven by the formation of the enzyme-substrate complex ES, which is reversible.The deterministic version of such reactions is where we write S • E to distinguish the multiplication of E and S from complex ES.By the relations a QSSA reduces to one the number of involved equations.Indeed, since ES is in quasi-steady-state, i.e.ES ′ = 0, then Here K M is termed the Michaelis-Menten constant.In practice, the QSSA permits to reduce the threereactions model to the single-reaction model S − → P Figure 2. Stochastically perturbed Enzyme-Substrate-Product system.Averages of 1000 simulations for P , plotted with its standard deviation (dotted), for both exact and approximated Michaelis-Menten models.Parameters are as in Figure 1 (i) in the top four panels, and as in Figure 1 (ii) in the bottom four.Independent Sine-Wiener noises are present in all the reactions.For i = 1, 2, 3, in the left column panels τ i = 1, in the right τ i = 5, in top β i = 0.5 and in bottom β i = 1.
with non mass-action non linear rate (V max S)/(K M + S).In [46] the condition is used to determine a region of the parameters space guaranteeing the legitimacy of the Michaelis-Menten approximation.When condition (31) holds, a separation exists between the fast pre-steady-state and the slower steady-state timescales [66] and the solution of the Michaelis-Menten approximation closely tracks the solution of the exact model on the slow timescale.
Here we show that the same condition is sufficient to legitimate the Michaelis-Menten approximation with bounded noises arbitrarily applied to any of the involved reactions.We start by recalling the result in [46] about the noise-free models given in Table 2.We considered two initial conditions: (i) one with 10 copies of substrate, 1 enzyme and 0 complexes and products, and (ii) one with 10 copies of substrate, 100 enzyme and 0 complexes and products.As in [46] we set c 1 = c 3 = 1 and c 2 = 10; notice that the parameters are dimensionless and, more important, in (i) they satisfy condition ( 31) since E T = 1 and S 0 + K M = 21, in (ii) no.In Figure 1 we reproduced the results in [46] for (i) in right panel and (ii) in left.As expected, in (i) the approximation is valid on the slow time-scale, and not valid in the fast, i.e. for t < 3, in (ii) it is not valid also in the slow time-scale.
If noises are considered the models in Table 2 change accordingly.So, for instance when independent Sine-Wiener noises are applied to each reaction, the exact model becomes and the Michaelis-Menten constant becomes the time-dependent function .
Notice that the nonlinear approximated propensity a 1 (t) is now time-dependent, and, moreover, it depends nonlinearly on the noises affecting the system.  .Stochastically perturbed Enzyme-Substrate-Product system.Averages of 1000 simulations for P , plotted with its standard deviation (dotted) in the fast time-scale 0 ≤ t ≤ 5 for both exact and approximated Michaelis-Menten models with parameters as in Figure 1 (i).Here we use independent Sine-Wiener noises with τ 1 = 1 in left and τ 1 = 5 in right panels.For every parameters configuration both low, i.e. 0.5, and high, i.e. 1, noise level intensities are used.
Thus condition (31) becomes time-dependent and we rephrase it to be Note that if β 1 > 0 then K * M (t) = K M , whereas if β 1 = 0 then K * M (t) = K M .Each of the shown figures is the result of 1000 simulations for model configuration where the simulation times, which span from few seconds to few minutes, depend on the noise correlation.When the same system of Figure 1 (i) is extended with these noises the approximation is still valid, as shown in the top panels of Figure 2. In addition, the approximation is not valid when condition (32) does not hold, as shown in the bottom panels of Figure 2, as it was in Figure 1 (ii).Notice that in there we use two different noise correlations, i.e. τ i = 1 in the left and τ i = 5 for i = 1, 2, 3 in the right column panels, thus mimicking noise sources with quite different charateristic kinetics.Also, we set two different noise intensities, i.e. β i = 0.5 in top panels and β i = 1 (maximum intensity) in bottom panels, whereas all the other parameters are as in Figure 1.Summarizing, we get a complete agreement between enzymatic reactions with/without noise, independently on the noise characteristics when it affects all of the reactions.
To strengthen this conclusion it becomes important to investigate whether it still holds when noises affects only a portion of the network and, also, whether it holds on the fast time-scale.
As far as the number of noises is concerned, we investigated various single-noise configurations in Figure 3.In there we used a single noise, i.e. two out of the three noises have 0 intensity, with both low and high intensities, i.e. 0.5 and 1.Also, in that figure we vary the noise correlation time as τ ∈ [1,5].As hoped, the simulations show that the approximation is legitimate in the slow time-scale for all the various parameter configurations, thus independently on the presence of single or multiple noises.
Finally, as far as the legitimacy of the approximation in the fast time-scale is concerned, i.e. t ∈ [0, 5], our simulations show a result of interest: if the noise correlation is small compared to the reference fast time-scale and if single noises are considered the noisy Michaelis-Menten approximation performs well also on the fast time-scale.We remark that this was not the case for the analogous noise-free scenario in Figure 1 (i).In support of this we plot in Figure 4 the fast time-scale for τ i = 1 and τ i = 5 for the single noise model with a noise in the enzyme-substrate complex formation, i.e. β 2 = β 3 = 0. Similar evidences were found in the configurations plotted in Figure 3 (not shown).
Table 3. Futile cycle model.The noise-free enzymatic futile cycle [22]: the stoichiometry matrix (rows in order S 0 , E, ES 0 , S 1 , F , F S 1 ) and the propensity functions.

Futile cycles
In this section we consider a model of futile cycle, as the one computationally studied in [22].The model consists of the following mass-action reactions where E and F are enzymes, S 0 and S 1 substrate molecules, and ES 0 and F S 1 the complexes enzymesubstrate.Futile cycles are an unbiquitous class of biochemical reactions, acing as a motif in many signal transduction pathways [68].
Experimental evidences related the presence of enzymatic cycles with bimodalities in stochastic chemical activities [69].As already seen in the previous section, Michaelis-Menten kinetics is not sufficient to describe such complex behaviors, and further enzymatic processes are often introduced to induce more complex behaviors.For instance, in deterministic models of enzymatic reactions feedbacks are necessary to induce bifurcations and oscillations.Instead, in [22] it is shown that, although the determinstic version of the model has a unique and attractive equilibrium state, stochastic fluctuations in the total number of E molecules may induce a transition from a unimodal to a bimodal behaviour of the chemicals.This phenomenon was shown both by the analytical study of a continuous SDE model where the random fluctuations in the total number of enzyme E (both free and as a complex with S) is modeled by means of a white gaussian noise on the one hand, and in a totally stochastic setting on the other hand.In the latter case it was assumed the presence of a third molecule N interacting with enzyme E according to the following reactions By using N the stochastic model results to be both quantitatively and qualitatively different from the deterministic equivalent.These differences serve to confer additional functional modalities on the enzymatic futile cycle mechanism that include stochastic amplification and signaling, the characteristics of which depend on the noise.Our aim here is to investigate whether bounded noises affecting the kinetic constant, and thus not modifying the topology of the futile cycle network, may as well induce transition to bimodality in the system behavior.To this aim, here we analyze three model configurations: (i) the noise-free futile cycle, namely only the first six reactions, (ii) the futile cycle with the external noise as given by N and (iii) the futile cycle with a bounded noise on the binding of E and S 0 , i.e. the formation of ES 0 , and N is absent.In Table 3 the noise-free futile cycle is given as a stoichiometry matrix and 6 mass-action reactions.The model simulated in [22] is obtained by extending the model in the table with a stoichiometry matrix containing N and four more mass-action reactions.For the sake of shortening the presentation we omit to show them here.The model with a bounded noise in a 1 is obtained by defining We simulated the above three models according to the initial condition used in [22] (S 0 , E, ES 0 , S 1 , F, F S 1 ) = (0, 20, 0, 2000, 50, 0) which is extended to account for 10 initial molecules of N , when necessary.The kinetic parameters are dimensionless and defined as k 1 = 40, k −1 = k 2 = 10000, k 3 = 200, k −3 = 100, k 4 = 5000 for the noise-free and the bounded noise case, and k 5 = k 6 = 10, k −5 = 5 and k −6 = 0.2 when the unimodal noise is considered [22].Furthermore, when the bounded noise is considered the autocorrelation is chosen as τ ∈ [k −1 5 , 1] = [0.1,1] according to the highest rate of the reactions generating the unimodal noise.
In Figure 5 a single run and averages of 1000 simulations for the futile cycle models are shown.In this case the simulation times span from 20 ÷ 30 min to 60 ÷ 80 min, thus making the choice of good parameters more crucial than in the other cases.In Figure 5 the substrate S 0 is plotted, and S 1 behaves complementarily.In top panels the noise-free (top) and the cycle unimodal noise as N (bottom).In .Stochastic models of futile cycles.Empirical probability density function for S 0 at t = 10 after 1000 simulations for the futile cycle models with the parameter configurations considered in Figure 5.
Table 4.The bistable model of gene expression in [73]: the stoichiometry matrix (rows in order R 1 , R 2 , P 1 , P 2 ) and the propensity functions.
bottom panels the cycle with bounded noise and autocorrelation τ = 0.1 in (left) and τ = 1 in (right).In both cases in the top panel the noise intensity is β = 0.5 (top) and β = 1 (bottom).The initial configuration is always (S 0 , E, ES 0 , S 1 , F, F S 1 , N ) = (0, 20, 0, 2000, 50, 0, 10) and the kinetic parameters are k 1 = 40, k −1 = k 2 = 10000, k 3 = 200, k −3 = 100, k 4 = 5000 for the noise-free and the bounded noise case, and k 5 = k 6 = 10, k −5 = 5 and k −6 = 0.2 [22].We also show in Figure 6 the empirical probability density function for the concentration of S 0 , i.e.P(X(t) = x) given the considered initial configuration, at t ∈ {2, 5, 7, 10} after 1000 simulations for the futile cycle models with the parameter configurations considered in Figure 5.The analysis of such distributions outline that for the noise-free system the distributions are clearly unimodal, whereas for noisy futile cycle, in both cases, they are bimodal.Moreover, it is important to notice that the smallest peak of the distribution, i.e. the rightmost, has a bigger variance when N is considered, rather than when a bounded noise is considered.

Bistable kinetics of gene expression
Let us consider a model by Zhdanov [13,73] where two genes G 1 and G 2 , two RNAs R 1 and R 2 and two proteins P 1 and P 2 are considered.In such a model synthesis and degradation correspond to In bottom panels averages of 1000 simulations.In all cases α R = 100 min −1 , α P = 10 min −1 , δ R = δ P = 1 min −1 , K = 100 and τ = 100 min −1 and the initial configuration is (R 1 , P 1 , R 2 , P 2 ) = (10, 0, 0, 0).The populations and the noise are plotted for the single run.
Such a reaction scheme is a genetic toggle switch if the formation of R 1 and R 2 is suppressed by P 2 and P 1 , respectively [7,14,[70][71][72].Zhdanov further simplifies the schema by considering kinetically equivalent genes, and by assuming that the mRNA synthesis occurs only if 2 regulatory sites of either P 1 or P 2 are free.The deterministic model of the simplified switch when synthesis is perturbed is where the perturbation is Here α R , δ R , α P and δ P are the rate constants of the reactions involved, term [K/(K + P i )] 2 is the probability that 2 regulatory sites are free and K is the association constant for protein P .Notice that here perturbations are given in terms of a time-dependent kinetic function for synthesis, rather than a stochastic differential equation.Before introducing a realistic noise in spite of a perturbation we perform some analysis of this model.As in [73] we re-setted model (33) in a stochastic framework by defining the reactions described in Table 4. Notice that in there two reactions have a time-dependent propensity function, i.e. a 1 (t) and a 3 (t) modeling synthesis.
In the top panels of Figure 7 we show single runs for Zhdanov model where simulations are performed with the exact SSA with time-dependent propensity function2 .We considered an initial configuration with only 10 RNAs R 1 .As in [73] we set α R = 100 min −1 , α P = 10 min −1 , δ R = δ P = 1 min −1 , K = 100 and τ = 100 min −1 ; notice that this parameters are realistic since, for instance, protein and mRNA degradation usually occur on the minute time-scale [74].We considered two possible noise intensities, i.e. α = 0.5 in left and α = 1 in right and, as expected, when α increases the number of switches increases.To investigate more in-depth this model we performed 1000 simulations for both the configurations.In the bottom panels of Figure 7 the averages of the simulations are shown.The average of our simulations evidences a major expression of protein P 2 against P 1 , for both values of α, with dumped oscillations for α = 0.5 and almost persistent oscillations for α = 1.
In Figure 8 we plot the empirical probability density function of the species concentrations, i.e.P[X(t) = x] given the considered initial configuration, at t ∈ {900, 950, 1000} min as obtained by 1000 simulations.Interestingly, these bi-modal probability distributions immediately evidence the presence of stochastic bifurcations in the more expressed populations R 2 and P 2 .In addition, the distributions for the protein seem to oscillate with period around 100, i.e. for α = 1 they are unimodal at t ∈ {900, 1000} and bi-modal at t = 950.
For the sake of confirming this hypothesis in Figure 9 the probability density function of P 2 is plotted against time, i.e. the probability of being in state x at time t, for any reachable state x and time 900 ≤ t ≤ 1000.In there we plot a heatmap with time on the y-axis and protein concentration on the x-axis; in the figure the lighter gradient denotes higher probability values.Clearly, this figure shows the oscillatory behavior of the probability distributions for both value of α and, more important, explains the uni-modality of the distribution at t = 900 and t = 1000 with α = 1, i.e. the higher variance of the rightmost peak at α = 1 makes the two modes collapse.Finally, we omit to show but, as one should expect, the oscillations of the probability distribution, which are caused by the presence of a sinusoidal perturbation in the parameters, are present and periodic over all the time window 0 ≤ t ≤ 1000.
Bounded noises.We investigated the effect of a Sine-Wiener noise affecting protein synthesis rather than a perturbation, i.e. a new ξ(t) is considered with W a Wiener process.Here simulations are performed by using the SSAn where the reactions in Table 4 are left unchanged, and the propensity functions a 1 (t) and a 3 (t) are modified to In bottom panels the averages of 1000 simulations.In all cases α R = 100 min −1 , α P = 10 min −1 , δ R = δ P = 1 min −1 , K = 100 and τ = 100 min −1 and the initial configuration is (R 1 , P 2 , R 2 , P 2 ) = (10, 0, 0, 0).The populations and the noise are plotted for the single runs.
For the sake of comparing the simulations with those in Figures 7-9, we used the same initial condition and the same values for α R , α P , δ R , δ P and K. To compare the effect of a realistic noise against the original perturbation we simulated the system with the same values i.e. the noise intensity α = 0.5 in left and α = 1 in right of the top panels in Figure 10, and in both cases τ = 100.As expected, in this case the trajectories are more scattered than those in Figure 7, and the switches are still present.However, for maximum noise intensity α = 1 time-slots emerge where the stochastic systems predicts a more complex outcome of the interaction.In fact, for t ∈ [0, 200] ∪ [800, 900] neither protein P 1 nor P 2 seem to be as expressed as in the other portions of the simulation, thus suggesting the presence of noise-induced equilibria absent when periodic perturbations are present.
To investigate more in-depth this hypothesis we again performed 1000 simulations for both the configurations, the averages of which are shown in the bottom panels of Figure 10.In this case, the simulation times, which again depend on the noise correlation, span from 3 ÷ 5 min to 30 ÷ 40 min, thus making the choice of good parameters crucial.Differently from the case in which a sinusoidal perturbation is considered, i.e. Figure 7, in this case the averages are not oscillatory, but instead show a stable convergency.Also, the final outcome seems again to predict the expression of P 2 inhibiting P 1 .To understand better this point we plotted in Figure 11 the probability density of reachable states at t = 1000 min, i.e.P[X(t) = x] given the considered initial configuration, and in Figure 12 we plotted that distribution against time for R 2 .It is worth noting that we also ranged t over [900, 1000] but since P[X(t) = x] did not change we omitted to plot it here.Again, Figure 12 is a heatmap where on the y-axis time in minutes is given, on the x-axis the possible concentration for R 2 and the lighter gradient denotes higher probability values.Notice that in this case Figure 12 represents an empirical evaluation of the solution the DCKE for this system, i.e. equation (14).Both graphics are obtained by 1000 simulations with α = 0.5 (left panels) and α = 1 (right panels).These figures show that a low-intensity noise makes the probability distribution become three-modal, i.e. notice the two rightmost peaks in Figure 11 and the white/light-blue gradients in Figure 12.Differently, when the noise intensity is higher, the two rightmost peaks almost merge, thus forming a bi-modal distribution where the smaller peak almost spreads uniformly on the state space for the variables.Notice that, in this case, the amplitude of such a peak is higher than for α = 0.5, i.e. notice the intensity of the blue gradient in Figure 12.For α = 0.5 it is possible to notice two red gradients: one approximatively for x → 200 and one for x ∈ (10, 30).The major peaks in the distribution for R 2 are for x < 10, for x ∈ (50, 100) and for x ∈ [130, 180].The probability of each of these peaks is decreasing as x increases, thus confirming the intuition of Figure 11.Similar considerations can be done when α = 1 where, as shown by Figure 11, the first dark-red area separating the first two peaks in α = 0.5 is vanished, thus forming a bi-modal instead that a three-modal probability distribution.Finally, for the sake of considering a wide range of biologically meaningful values for τ , which we recall it represents a measure of the speed of noise variation, we evaluated the solution of the DCKE for R 2 for the same configuration used in Figure 12 and τ ∈ {1, 10, 25, 100} min.We performed 1000 simulations of the model for each value of τ with α = 0.5, the value showing a more interesting behavior.In Figure 13 the probability of the reachable states at t = 1000 min is plotted.If is immediate to notice that the height of the first peak increases as τ decreases, and more precisely the distribution seems to switch from a three-modal one to a bi-modal when τ ≤ 25.In each panel of Figure 14 we plot the variation of such probability distribution for 900 ≤ t ≤ 1000.By that figure it is possible to observe that by ranging τ the dark-red gradient increases in size as far as τ decreases.This means that the amplitude between the peaks of the density strictly depends on the value of τ , thus suggesting a strong role for extrinsic noise in determining the network functionalities.In the x-axis the species concentration is represented, in the y-axis minutes are given.

Discussions
In this paper we investigated the effects of joint extrinsic and intrinsic randomness in nonlinear genetic networks, under the assumption of non-gaussian bounded external perturbations.Our applications have shown that the combination of both intrinsic and extrinsic noise-related phenomena may have a constructive functional role also when the extrinsic noise is bounded.This is in line with other researches -only focusing on either intrinsic or extrinsic noise -recasting the classical interpretation of noise as a disturbance more or less obfuscating the real behavior of a network.This work required the combination of two well-known frameworks, often used to separately describe biological systems.We combined the theory of stochastic chemically reacting systems developed by Gillespie with Langevin systems describing the bounded variations of kinetic parameters.The former shall allow considering the inherent stochastic fluctuations of small numbers of interacting entities, often called intrinsic noise, and clearly opposed to classical deterministic models based on differential equations.The latter permits to consider the influence of bounded extrinsic noises.These noises are modeled as stochastic differential equations.For these kind of systems, although an analytical characterization is unlikely to be feasible, we were able to derive a differential Chapman-Kolgomorov equation (DCKE) describing the probability of the system to occupy each one of a set of states.Then, in order to analyze these models by sampling from this equation we defined an extension of the Gillespie's Stochastic Simulation Algorithm (SSA) with a state-dependent Langevin system affecting the model jump rates.This algorithm, despite being more costly than the classical Gillespie's SSA, allows for the exact simulation of these doubly stochastic systems.
We outlined the role of extrinsic noise for some biological networks of interest.In particular, we were able to extend classical results on the validity of the Michaelis-Menten approximation to the prototypical Enzyme-Substrate-Product enzymatic reaction by drawing a Stochastic Quasi Steady State Assumption (SQSSA) for noisy reactions.Along the line of the classical deterministic or stochastic uses of the Michaelis-Menten approximation, this should permit to reduce the size of more general enzymatic networks even in presence of extrinsic bounded noises.
Moreover, we showed that in a recurrent pattern of genetic and enzymatic networks, i.e. the futile cycle, the presence of extrinsic noises induces the switching from a monomodal probability density (in absence of external perturbations) to a multimodal density.
Similarly, in the case of the toggle switch, which is inherently multistable, the presence of extrinsic noise significantly modulates the probability density of the genes concentration.In this important network motif we also investigated the role of periodic perturbations against a realistic noise.
Thus in general the co-presence of both intrinsic stochasticity and bounded extrinsic random perturbations might suggest the presence of possibly unknown functional roles for noise for these and other networks.The described noise-induced phenomena are shown to be strongly related to physical characteristics of the extrinsic noise such as the noise amplitude and its autocorrelation time.
A relevant issue that we are going to investigate in the next future is the role of the specific extrinsic bounded perturbations.Indeed, in non-genetic systems affected by bounded noises it has been shown that the effects of the perturbations depend not only on the above general characteristics of the noise, but also on its whole model [40][41][42]76].In other words the transitions of a system perturbed by a sine-Wiener noise might be quite different from those induced by another bounded perturbation, for example the Cai-Lin noise [77] or the Tsallis noise [43], also when their amplitude and autocorrelation times are equal.In other words, a single network in two different environments might show two different behaviors depending of fine details of the kind of perturbations that are present.This might also suggest that a same network might exhibit many different functions depending on its "locations".
Finally, concerning these points, we stress that these peculiar properties of bounded extrinsic perturbations make it even more important the investigations, such as those of [37], aimed at inferring by deconvolution the external noise from the experimental data, in order to infer which kind of noise affect a given network in a well determined environment.

(Figure 3 .
Figure 3. Stochastically perturbed Enzyme-Substrate-Product system.Averages of 1000 simulations for P with dotted standard deviation, for both exact and approximated Michaelis-Menten models with parameters as in Figure1(i).Here single noises with two intensities and different autocorrelations are used.The non-zero parameters are reported in top captions.

Figure 4
Figure 4. Stochastically perturbed Enzyme-Substrate-Product system.Averages of 1000 simulations for P , plotted with its standard deviation (dotted) in the fast time-scale 0 ≤ t ≤ 5 for both exact and approximated Michaelis-Menten models with parameters as in Figure1(i).Here we use independent Sine-Wiener noises with τ 1 = 1 in left and τ 1 = 5 in right panels.For every parameters configuration both low, i.e. 0.5, and high, i.e. 1, noise level intensities are used.

1 Figure 6
Figure 6.Stochastic models of futile cycles.Empirical probability density function for S 0 at t = 10 after 1000 simulations for the futile cycle models with the parameter configurations considered in Figure 5.

α = 0.5 α = 1 Figure 9 .
Figure 9. Periodically perturbed toggle switch.Empirical probability density function for R 2 plotted against time, i.e. the probability of being in any reachable state x for 900 ≤ t ≤ 1000.Lighter gradient denotes higher probability values.We used data collected with 1000 simulations of model(33) where τ = 100 and two perturbation intensities are used, i.e. α ∈ {0.5, 1} as reported in the top captions.In the x-axis the species concentration is represented, in the y-axis minutes are given.

Figure 11 .
Figure 11.Stochastically perturbed toggle switch.Empirical probability density function at t = 1000, after 1000 simulations for Zhdanov model with Sine-Wiener noise.Parameters are as in Figure 10 and two perturbation intensities are used, i.e. α ∈ {0.5, 1} as reported in the top captions.

Figure 14 .
Figure 14.Stochastically perturbed toggle switch.Empirical probability density function for R 2 plotted against time, i.e. the DCKE solution for R 2 in 900 ≤ t ≤ 1000.Lighter gradient denotes higher probability values.We used data collected with 1000 simulations of Zhdanov model with Sine-Wiener noise.From left to right τ = 1, τ = 10, τ = 25 and τ = 100.In all cases the noise intensity is α = 0.5.In the x-axis the species concentration is represented, in the y-axis minutes are given.