Gene switching rate determines response to extrinsic perturbations in the self-activation transcriptional network motif

Gene switching dynamics is a major source of randomness in genetic networks, also in the case of large concentrations of the transcription factors. In this work, we consider a common network motif - the positive feedback of a transcription factor on its own synthesis - and assess its response to extrinsic noises perturbing gene deactivation in a variety of settings where the network might operate. These settings are representative of distinct cellular types, abundance of transcription factors and ratio between gene switching and protein synthesis rates. By investigating noise-induced transitions among the different network operative states, our results suggest that gene switching rates are key parameters to shape network response to external perturbations, and that such response depends on the particular biological setting, i.e. the characteristic time scales and protein abundance. These results might have implications on our understanding of irreversible transitions for noise-related phenomena such as cellular differentiation. In addition these evidences suggest to adopt the appropriate mathematical model of the network in order to analyze the system consistently to the reference biological setting.

1 Mathematical and computational background 1.1 Bounded noises Temporal bounded noises can be synthesized either via stochastic differential equations, e.g. [WT04,CW04], or via application of a bounded function to a random process, e.g. [BC05,Dim88].
We recall the definition and some basic properties of the two noises that are used in our applications, further details are available in [d'O13].
• Cai-Lin bounded noise. Consider the following Langevin equatioṅ where η(t) is a Gaussian white noise. As it is easy to verify, if ξ(0) ∈ [−B, +B] then ξ(t) ∈ [−B, +B]. Moreover, it has zero mean and the same stationary autocorrelation of the Ornstein-Uhlenbeck process. The steady-state probability density of this bounded stochastic process is The above stationary probability density exhibits different shapes according to the values of the parameter z. For z > 0 it is unimodal and centered in 0, while for −1 < z < 0 it is bimodal with two vertical asymptotes at ξ → ±B (i.e., it is "horned").
• Sine-Wiener bounded noise. This is obtained by applying the bounded function B sin(.) to a Wiener process W (t) yielding The noise steady-state distribution is horned and equal to the one of the Cai-Lin case for z = −0.5. In Figure S1 some typical steady-state distributions from Cai-Lin and sine-Wiener approaches are depicted.
Noise autocorrelation. Regardless the type of noise considered, a fundamental concept its autocorrelation, also known as serial correlation, namely the cross-correlation of a signal (a temporal series) with itself. Informally, it is the similarity between observations as a function of the time-lag between them. Let suppose the noise has at timet a value ξ(t). Autocorrelation characteristic time τ measures the time window in which the noise has a tendency to "remember" its past history, so up tot + τ the value of ξ(t) is somewhat similar to ξ(t). When the noise is totally uncorrelated, like is the case of Gaussian white noise, the value ξ(t+dt) is totally independent from ξ(t) for any dt (in other words one can say that τ → 0).

Gillespie's algorithm for systems perturbed by bounded noises
The simulation of a set of chemical reactions described by a stoichiometry vector ν i and with propensity function a j (·) is performed by using Gillespie's algorithm [Gil77], here adapted to support time-inhomogenous propensity functions arising either from the presence of noise or from the approximation discussed in the Main Text. A detailed background on such extension can be found, e.g., in [CMd13]. The translation of a birth-death process to the algebraic representation is straightforward; see [Gil77].
When the system state at time t is x t = [x gene x prot x ξ ] , where x gene is the value of G (number of active genes), x prot is the number of/density of y/Y (number of proteins), and x ξ is the noise scalar value CL z=+0.5 CL z=-0.5 SW Figure S1. Stationary distributions of the Cai-Lin (horn-shaped) and the sine-Wiener (bell-alike) bounded noises suggest different types of perturbations, according to the parameters.
(e.g., such as one of those described in §1.1), the probability density function for the time of the next jump follows for an initial state x 0 at time t 0 , with a 0 (t + w) = j a j (x t+w ). Notice that these propensity function are computed in the state that the system will be, at times afterwards t, e.g., x t+w , which will be generally different from x t as the system is not homogenous in time. Given τ , the event of firing reaction j has density When noise is present, or when the system is hybrid -Models 2 and 3 -an evolution equation for noise/proteins (Table 2, Main Text) is coupled to common strategies to samples from such distributions, as described for instance in [CMd13]. The original Gillespie approach has scalability issues when x prot is high. In particular, the total exit rate of the process is dominated by feedback and degradation, i.e. a 0 (t) ∝ a feed (t) + a degr (t) which grows as O(x 2 prot ) when x prot is high, as x gene < 2. If we assume noise to change smoothly, we can consider the propensities to be time-homogenous and observe that τ = a 0 (t) −1 .
Hence the average time-increment of each jump is small (precisely, decreases quadratically with the number of proteins), and it might get computationally hard to sample many independent trajectories of the system, unless we switch to an hybrid representation of x prot . In this case, x prot varies via a differential equation, the total contribution of the rates to a 0 (t) is smaller so allowing the system to leap over longer jumps, with the overhead of integrating such a differential equation.
2 Biological details of the adopted model The simple model that we adopted, introduced in literature by Smolen, Baxter and Byrne in their well-known paper [SBB98], represents a single transcription factor (TF, variable Y ) that self-activates. In detail: 1. the gene producing the TF has only one promoter; 2. the active form of the TF is its dimer; 3. the dimerization rate is very fast, thus the dimer is at equilibrium with the TF and the dimerization process is not directly represented in the model; 4. the activation of the gene is thus triggered by its binding with the dimer; 5. however, also in absence of binding there is a small baseline activation rate, independent from the gene binding to dimers.
Besides the assumptions made in [SBB98], we further assume that: 6. due to the interplay with unknown signals (modelled as bounded noise) the deactivation rate stochastically fluctuates.
The mathematical translation of this biological scenario leads to a stochastic model with total activation rate in bimolecular form which, according to the the relation between reaction rate equations and propensity functions for stochastic systems -see, e.g., [Gil77] -leads to the approximation in the case of large amount of the TF as We used this form of the rate in all settings, as stochastic models where the dimerization follows law of the type Y 2 instead of Y (Y − 1) are widely used in the literature [BG05, WSW16].

Detailed analyses
In this section, we drop the bar over the symbolsc 0 ,c 2 andb 0 . This is an abuse with respect to the main text, but it is done for the sake of the simplicity of notation.

Analytical results: fast gene-switching and large protein counts
Concerning model D (fast gene switching and large protein counts), we may give some analytical results of interest. Here we assume that the parameters of the FC model are such that the system is multistable. We define the utility functions which compute, respectively, the smallest, the intermediate and the largest of the three real solutions of the equation: . This results follows from the following differential inequalities: Indeed, the above inequalities imply that y a (t) < y(t) < y b (t) where y a (t) solve the following ODEs: As it is easy to verify, these ODEs are such that for large times y a ≥ F (b * (1 + B)) and and y(0) < F (b l ) it apparently follows the quite neutral result that for large times it must be Due to the random nature of the perturbations, the existence of such a t is very likely, if not sure. Indicating with and with P st L (Y ) = P L (Y, +∞) the above results suggest that in the case The second result suggests that defining the "order parameter" y , and considering its variation with B ≥ 0, it may undergo a first order transition at b * (1 − B) = b l , because there y (B) is discontinuous. A similar result can be obtained with reference to the upper branch of the bifurcation diagram. Indeed, With respect to the unconditioned probability density P(Y, t) the above results show that its asymptotic behavior (in the functional space of the probability densities) strongly depends on its initial conditions. In other words, for small to moderate amplitude B of the noise y(t) has two distinct stochastic attractors (two distinct phases). By increasing B over a threshold B C the asymptotic behavours changes and there is a unique attractor. Of course, if the nosie amlitude is slightly over the threshold and then decrase to a value under the threshold, one passes for a scenario with a unique attractors to another one with two attractors. For example, in theory, in numerical simulations with y(0) ∈ [0, F (b l )] we would expect to observe such a first order transition, i.e. that it exists a B c (with b * (1 − B c ) = b l ) such that for 0 ≤ B < B c the stationary density P st L (Y ) is unimodal and located at low values of y, whereas for B c < B the density jumps at larger values of y.
As far as this first-order transition due to the double distinct attractors of y(t) that disappear when the parameter B increases and crosses a threshold, we recall that the non-uniqueness of the asymptotic attractors in a physical system (and, as a consequence, in its associated Fokker-Planck equation) is a landmark of genuine phase transitions, as discussed in depth by Shiino [Shi87] (see also [VdBPT94]). This phenomenon is rarely observed in non-spatial systems, and it makes our observed first-order transition quite close to 'true' phase transitions.
We note that in the statistical physics the search for multiple attractors has traditionally been done in non-linear Fokker-Planck equations, obtained by applying mean field approximations to large number of interacting nonlinear systems (e.g. by diffusion [Shi87,VdBPT94]). In our case, instead, the system is unique and non-spatial, and the

Time-varying Waddington potential
Recently Verd et al. [VCJ14] introduced in the framework of Systems Biology models affected by deterministic time-varying perturbations the concept of time-varying Waddington's potential. We believe that this potential might be extremely useful also for stochastically perturbed systems like ours. Namely, to our model: it is associated the following time-varying potential: The shape of the potential W (x, t), and the number and basin of attraction of its 'potential holes', stochastically change in the time. Thus, for example, in case of small to moderate fluctuations of b 0 (t) the corresponding irreversible transition 'low to large' ('large to low') values of y can heuristically be read as the irreversible 'capture' by large (low) value attractors of a trajectory initially confined in a potential hole centered at low values of y.

Sensitivity analysis on protein transcription
When the continuous model is in it multi-stability region for non-degenerate s results shown in the Main Text are unchanged. Consider in fact the noise-free continuous equations in model D (Table 2, Main Text), replace G(t) in the differential equation for proteins.
When ∆(s) > 0 then p(y, s) has 3 roots; numerical conditions for this to hold are This means that for any s ∈ [2.2255, 143.7882] we can find conditions to have an hysteresis plot equilibria such as the one in Figure 1. Such plots will only differ for the values of the equilibrium values (y low , y high ) and (b l , b r ), as a function of s.
Notice that we picked -on purpose -a set of parameters values in which all the 3 roots exists, as this is the interesting case where multi-stability can emerge. Thus, for any of such plots, we can always retrieve our conclusions when the noise baseline value (b 0 ) is close enough to b l -which is the situation described in the Main Text. For this reason, we did not include any quantitative statement on our predictions (like specific protein counts) being much aware on the dependency on the specific value of s.

Multiple noises on gene switching at once: nullification effects
We show that, if two independent noise sources are applied to both the process of gene switching between inactive and active states and viceversa the extrinsic noise effect is averaged out.
Let us consider the continuous case, which is model D of our manuscript; the equation for the average number of active genes can be extended to account for a noise in the activation baseline rate, c 0 in the main text, yielding G(t) new = n c 0 (t) + c 2 y 2 (t) b 0 (t) + c 0 (t) + c 2 y 2 (t) Let us consider the two noises separately. First, applying the bounded noise in b 0 , we would have that b 0 (t) = b 0 and that, for any t, b 0 (t) ∈ [b 0 (1 −B), b 0 (1 +B)] which would lead, case of the strongest possible noise intensity (B = 1) to b 0 (t) ∈ [0, 2b 0 ]. In this case the state equation for gene number fluctuates between: With our parameter setting, i.e. c 0 = 10 c 2 = 16 b 0 = 150, we will obtain for a low level of protein equilibrium y eq low 1, G b0  When the noise is applied to c 0 , the state equation for gene number fluctuates between: The same consideration of the former case lead us, for y eq low to G c0 high . For further clarification, we plot G(t) in the case of the noise source affecting both activation and inactivation of G; see Figure S2. The same conclusions can be drawn when the same noise source is considered.

From continuous to discrete molecular counts
The discretization of a model do not necessarily preserve the properties of its continuous counterpart. Different volumes -and hence different normalization constants, Fig. 7 (Main Text), and discretizationimply different reference biological systems. In general, to the best of our knowledge, it is not possible to automatically derive results on discrete-level of proteins/genes from a phase-diagram. In practical terms, this means that there is no theoretical guarantee that standard analysis on the continuous version of the model (model D) can shed any light on the network functioning in different settings (models A,B and C). This is indeed what we observe for this system, see Figure S3.

Extended comparison with related works
In this section we comment on three models of the self-activation transcriptional network motif and its relation with extrinsic Gaussian noises. In two cases this is the only form of stochasticity considered [LaLGLL09, LJ04], in a unique case this is in combination with intrinsic fluctuations but without geneswitching noise [ARLSG13]. We can observe that the predicted behaviors are independent of the initial condition, but rather dependent on the cellular volume (via normalization N A V ). Notice that for low noise, B = 0.1, oscillations are only predicted for low-volume cells, N A V = 6.02. While for higher noise B = 0.1 no matter what the volume is the dynamics oscillates across high and low protein counts.
As we mentioned in the main text, this model in absence of extrinsic noise reads as follows: but, in the literature, e.g., [SBB98,LJ04,LaLGLL09,ARLSG13], it is often written in the algebraically equivalent formẏ where R b and the sum K f + R b can be, respectively, legitimately read as a baseline and an "asymptotic" production rates. However, since it holds that c 2 it follows that the parameters R b , K f and K d must not be dealt with as they were independent. In particular, a fluctuation of the baseline rate R b cannot be deconvolved by fluctuations in the other two parameters.
Unfortunately this is what happens in literature. For example, in [LaLGLL09] the Smolen-Baxter-Byrne model was studied in the framework of the above mentioned continuous approximation with fast gene-switching and large number of proteins. The authors of [LaLGLL09] analytically investigated the consequences of white stochastic oscillations affecting the baseline production rate K d and/or the parameter K f . Thus, the model in [LaLGLL09] reads as followṡ Unfortunately, according to what we above showed, the stochastic differential equation (7) does not correspond to a biologically meaningful scenario, unless both ξ 0 (t) = ξ 1 (t) and, of course, the noises are bounded. We note here that, in such a particular case, eq. [LaLGLL09] can be read as a model of fluctuations in the parameter s.
In [LJ04] it is investigated a model where both the baseline protein production rate and the degradation rate were perturbed by white noises, yielding the following stochastic differential equatioṅ Again, isolated oscillations of the parameter R b are not meaningful. As briefly mentioned in the introduction, by means of semi-analytical methods, joining the Wentzel-Kramers-Brillouin (WKB) approximation and numerical simulations, Assaf and colleagues recently investigated the interplay between extrinsic and intrinsic noise in the circuit of a self-transcription factor with a sharp positive feedback [ARLSG13]. The biological differences between their model and ours are worth to be described in some detail.
Indeed, in their main model the extrinsic noise perturbs the production rate of the transcription factor, and the perturbation is state-dependent being active only if the state of the protein is "high". Indeed, under the implicit hypothesis that the gene switching velocity is large they assume the following probability law for the production of the transcription factor in the time interval (t, t + dt) where θ(.) is the Heavyside function and α 0 1. Thus if Y (t) is under the threshold then the production rate of the transcription factor is unperturbed.
The mechanism through which this state-dependent fluctuation of the production rate can be enacted is not specified. In absence of such specification, the asymmetry of the perturbation acting on the baseline and on the large protein synthesis rates remains unclear.
Note that in [ARLSG13] it is also briefly investigated a model where a smooth feedback is enacted (as in the Smolen-Baxter-Byrne model [SBB98,SBB99]): Thus, the perturbation adopted in [ARLSG13] for such a smooth case is equivalent to the perturbation of the parameter K f alone investigated in the paper [LaLGLL09] in the case ξ 0 (t) = 0. We note here that this equivalence can be extended to the 'sharp' model defined by 9 because in absence of extrinsic noise eq 9 is the limit case of a generalization of our model B (fast gene switching and small number of proteins) with constant b 0 . Indeed, in the hypothesis that gene switching is fast and that the activation is caused by 'Q-meric' forms (and not dimeric), proceeding as in our main text one gets (in absence of noise on b 0 ): which may be rewritten as follows: where Thus, if Q is sufficiently large eq. (12) reads: Summarizing, we may say that for both the sharp and the smooth models of [ARLSG13] the extrinsic noise acts in a biologically unrealistic way, equal (for Q = 2) or remindful (for Q 1) of that investigated in [LaLGLL09] in the case ξ 0 (t) = 0.
As far as the type of extrinsic noise is concerned, Assaf and coworkers considered an unbounded Orenstein-Uhlenbeck noise defined as followṡ where η(t) is a unitary white noise. In this way the stationary probability density function of the noise ξ is Gaussian with variance σ ex . However, from the above equations it follows that it must be One might roughly consider tolerable the error if which, however, does not seem the case in [ARLSG13]. Indeed in the simulations presented in [ARLSG13] the employed "stochastic bifurcation parameter" is the ratio between the standard deviation of the Orenstein-Uhlenbeck noise and a parameter named µ. This parameter is defined as the average of the (transitory) "quasi-stationary distribution about the high state" [ARLSG13]. Since in their simulations the high state is large (it ranges from 300 to 5000), even µ has to be large. Indeed, we simulated the sharp feedback model by means of the (exact) Gillespie algorithm and we got that the average value of the probability density function is very close to the high equilibrium state. This means that the standard deviation σ ex in most of the simulations reported in [ARLSG13] largely exceeds 0.5(1 − α 0 ).

Bounded versus unbounded forms of noise
Consider the continuous version of the model and the phase-transition that we observe. The fact that more common approaches based on unbounded extrinsic perturbations will not predict such dynamics is one of the motivations for our work. It is indeed the case that, for the parameters that we have considered, the network interplaying with unbounded white/colored noises predicts oscillations among high and low levels of protein Y . We picture evidences of this from single simulations of the system under the effect of one of the most famous models of unbounded noise, see Figure S4. We there show two single simulations of the system under the effect of an Ornstein-Uhlenbeck noise ξ generated from the following stochastic differential equation where we "compare" the Ornstein-Uhlenbeck noise variance (σ) to the intensity of the bounded noise. In support of our conclusions we cal also observe that, even for σ B * the system is predicted to stabilize to either its higher protein state, or to oscillate among equilibria -as a function of τ . We note that this happens with non-zero probability for all parameter settings of such noises, by the infinite support of the Gaussian distribution. As such, we feel that this kind of behaviours -emerging only with unbounded noises -could be considered as a modeling artifact.
For this reason, caution should be adopted to distinguish predictions more likely to be "realistic" (i.e., adhere to underlying physics of the system) from modelling artifacts. In this sense, we believe that the type of behaviour that we have shown, i.e., a phase-transition, shall be considered a peculiari noise-effect, as it is not whatsoever reproducible in canonical white/colored-noise modeling contexts.