Transient hysteresis and inherent stochasticity in gene regulatory networks

Cell fate determination, the process through which cells commit to differentiated states is commonly mediated by gene regulatory motifs with mutually exclusive expression states. The classical deterministic picture for cell fate determination includes bistability and hysteresis, which enables the persistence of the acquired cellular state after withdrawal of the stimulus, ensuring a robust cellular response. However, stochasticity inherent to gene expression dynamics is not compatible with hysteresis, since the stationary solution of the governing Chemical Master Equation does not depend on the initial conditions. We provide a quantitative description of a transient hysteresis phenomenon reconciling experimental evidence of hysteretic behaviour in gene regulatory networks with inherent stochasticity: under sufficiently slow dynamics hysteresis is transient. We quantify this with an estimate of the convergence rate to the equilibrium and introduce a natural landscape capturing system’s evolution that, unlike traditional cell fate potential landscapes, is compatible with coexistence at the microscopic level.

of the system that, unlike traditional cell fate potential landscapes, is compatible with the notion of coexistence at the microscopic level.

Introduction
In a deterministic description, binary decision making is attributed to the irreversible state transition between two mutually exclusive stable steady states in response to a signal. This state transition is usually governed by regulatory motifs with the capacity for bistability and hysteresis 1 , thus ensuring that the system does not switch back immediately when the signal is removed 2 .
The stochastic dynamic behaviour of a gene regulatory network is governed by a Chemical Master Equation (CME), which describes the time evolution of the probability distribution of the system state. The stationary solution of the CME is unique and independent on the initial state of the system 3 and therefore, incompatible with memory effects or hysteresis. The incompatibility of hysteresis with intrinsic noise in gene regulatory networks has been addressed, for example, by Lestas et al 4 . However, there are numerous works providing experimental evidence of hysteretic behaviour under significant levels of stochasticity [5][6][7][8] .
In the context of phenotypic switching and cell fate determination, three different scenarios have been distinguished and experimentally observed for binary decision making: deterministic irreversible 9-11 , stochastic reversible 12 and stochastic yet irreversible state transitioning 13 . Reversibility is understood here as the capacity of individual cells to switch back in absence of external signals. According to a pseudo-potential interpretation, dynamics are directed by a pseudo-potential landscape divided by a separatrix into two basins of attraction such that each local minimum corresponds to a specific cellular state. Stochastic irreversible transitions are found to appear when cells are initialized on (or near) the separatrix 13 .
In this article we provide a quantitative description of hysteresis and apparent irreversibility in stochastic gene regulatory networks at the single cell level as transient effects, which disappear at the stationary state. Our analysis is based on an accurate approximation of the CME. This means that our results are valid for purely stochastic regimes far from the thermodynamic limit, and thus complementary to those based on the classical linear noise approximation which are suitable for systems near the thermodynamic limit 4,14 . Since the stationary solution of the CME is unique 3 , if the solution corresponds to a bimodal distribution, state transitions at the single cells level occur necessarily in a random and spontaneous manner, switching back and forth between regions of high probability.
Fang et al 15 experimentally determined an energy potential-like landscape as the negative logarithm of the probability distribution, as well as the transition rates, based on previous theoretical studies 16 . In this contribution, we provide a theoretical basis that explains coexistence of different expression states. In fact, under the assumption of protein bursting 17 , we propose an efficient form of the CME 17, 18 that allows us to construct a meaningful probability based landscape.
Furthermore, a clear link between the characteristic kinetic parameters of regulation dynamics and the resulting landscape is established.

Results
We consider the simplest gene regulatory motif exhibiting hysteresis, a single gene with positive self-regulation (see Fig. S1). In its deterministic description, the evolution of the amounts of mRNA and protein X (m and x, respectively) for the self-regulatory gene network is given by the set of ODEs: where γ m and γ x are the mRNA and protein degradation rates, respectively. k m c(x) is the transcription rate, that is essentially proportional to the input function c(x) which collects the expression from the activated and inactivated promoter states. This function incorporates the effect of protein self-regulation and takes the form 19,20 : with ρ(x) being a Hill function 21 that describes the ratio of promoter in the inactive form as a function of bound protein: The above expression, can be interpreted as the probability of the promoter being in its inactive state, where K = k off /k on is the equilibrium binding constant and H ∈ Z\{0} is an integer (Hill coefficient) which indicates whether protein X inhibits (H > 0) or activates (H < 0) expression.
Finally, expression (3) includes basal transcription or leakage with a constant rate ε = k ε /k m (see Fig.1) typically much smaller than 1. The parameters of the Hill function employed along the paper are H = −7 (the value taken from 22 ) and K = 100, whereas ε = 0.05. Unless other value is indicated, we use a = 54. Assuming that mRNA degrades faster than protein X we have that m * = k m c(x)/γ m and model (1) reduces to: where τ = tγ x , a = k m /γ x and b = k x /γ m .
The self-regulatory network described by the deterministic equations (1-2) shows bistability and hysteresis (see Fig. 1). For a range of the control parameter, b, the system evolves towards one stable state or another depending on the initial conditions. We therefore say that the system has memory, since steady state values provide information about the system's past. In systems with hysteresis (dependency of the state of the system on its past), forward and reverse induction experiments follow different paths resulting in a hysteresis loop (the system switches back and forth for different values of the control parameter) 23 .
Gene expression is inherently stochastic. Taking into account that mRNA degrades faster than protein X in most prokaryotic and eukaryote organisms 24 , protein is assumed to be produced in bursts 19, 20, 25, 26 at a frequency a = k m /γ x , (see equation (5)). From this assumption, it follows 20, 26 that the temporal evolution of the associated probability density function p : R + × R + → R + can be described by a Partial Integro-Differential Equation (PIDE) of the form: parameter b below a given threshold, there is a unique stable steady state of low protein x towards which the system evolves independently of the initial conditions. For input signals above a second threshold, the system evolves towards a unique stable steady state of high x. For signal values within both thresholds, the system is bistable, and evolves towards one stable state or another depending on the initial conditions. In the bistability region, enclosed by two saddle-node bifurcations, three different steady states coexist (stable and unstable branches are depicted using solid and dotted lines, respectively) for a given b.
where x and τ correspond with the amount of protein and dimensionless time, respectively. The latter variable is associated to the time scale of the protein degradation, as in the previous deterministic description. In addition, ω(x − y) is the conditional probability for protein level to jump from a state y to a state x after a burst, which is proportional to: with b, as in equation (5), representing the burst size. The stationary form of the one dimensional equation (6) has analytical solution 19,20 It has been shown that the equilibrium solution associated to a CME is unique and stable 3 . This is also the case for the Friedman equation In order to compute an estimate of the convergence rate to equilibrium we make use of entropy methods 27, 28 and define the entropy norm as G is a convex function in u, that in this study has been chosen to be H . According to Pájaro et al 28 and Cañizo et al 27 , G satisfies the following differential inequality: with η being a positive constant related to regulation (parameters H and K), as well as the transcription-translation kinetics (a, b). The smaller η, the slower its convergence towards the corresponding equilibrium solution. Computing η requires a full simulation of (6) until the system reaches the equilibrium distribution for each parameter on a given range, what is computationally involved. In this work, the PIDE model (6) is solved by using the semilagrangian method implemented in the toolbox SELANSI 18 .
Alternatively, and in order to avoid simulation burden, we provide a truncation method to approximate the rate of convergence. The method makes use of the discrete jump process representation (see Fig. S3), which is a precursor of Friedman PIDE model, by making the protein amount a continuous variable 17 . Our method (See SI) provides a good approximation of the convergence rate η by the negative eigenvalue with smallest absolute value of the state change matrix M associated to the discrete jump process, truncated to an N maximum number of proteins, which we refer to as λ 1 . This proof of concept has served to clarify how hysteresis, as it is known in deterministic nonlinear systems (i.e. as a long term stationary phenomenon) has not an equivalence in a microscopic world governed by a CME. For stochastic systems, hysteretic behaviour is a transitory phenomenon, i.e. it can be only obtained under transients that may resemble stationary solutions due to the extremely slow dynamics at which bimodal distributions evolve.
Nonetheless, some correspondence can be drawn between the most frequently visited states on a microscopic system and the stable states on the deterministic counterpart. As it has been discussed in Pájaro et al 20 , the extreme states of a stationary bimodal distribution, namely those that include the highest and lowest probable states reached, satisfy: where ρ(x) is defined in (4). Making zero the right hand side of equation (5) and re-ordering terms, the set of all possible equilibria satisfies: Both expressions (9) and (10)   showing hysteretic behaviour at the level of gene regulatory networks: if the system is governed by the CME, hysteresis is necessarily transient. Note that for slow dynamics (high a and low b values), the time needed to reach the stationary state might be longer than the natural timescales of relevance to the process. This is in accordance with previous studies reporting large mean passage times 14 , and also with Wu et al. 13 where they engineer a synthetic switch with stochastic yet irreversible transitions (the same mutually inhibitory gene regulatory motif is analyzed in the supplementary material using our PIDE approach).
The characterization of a cell response as hysteretic or non-hysteretic is important. For example, in a a recent study concerning epithelial to mesenchymal transition (EMT), a process through which epithelial cells transdifferentiate into a mesenchymal cell fate, the authors characterize two types of responses, hysteretic and non-hysteretic EMT, and report the notable influence of hysteresis on the metastatic ability of cancer cells 31 . Essentially, what the picture shows is that, rather than a parameter-dependent preferential state among two stable ones, there are two highly probable states that coexist for a given parameter region on a cell population.
Invoking pseudo-potential concepts to interpret dynamics in GRN under fluctuations 13 , although attractive from an intuitive point of view, may be misleading since it cannot capture the notion of coexistence. By coexistence we mean that two different protein expression levels coinciding with the peaks of the bimodal distribution coexist on a cell population (assuming no cell to cell variability on the initial conditions).
The pseudo-potential landscape is not easy to compute either, specially when increasing the number of proteins expressed. Alternatively, we can use the stationary solution of (6) to construct on the natural framework of probability distributions, a landscape informing of the possible transitions or evolution of the underlying microscopic system. As we illustrate in the example discussed in the supplementary material, its computation can be extended in a straightforward manner to larger dimensional protein spaces. This can be of used to efficiently identify most prevalent phenotypes coexisting on a given cell population.
Data Availability: All relevant data needed to reproduce the results are included in the text and supplementary information.   Figure S1: Self-regulatory transcription-translation mechanism. The promoter is assumed to switch between active (DNA on ) and inactive (DNA off ) states, with rate constants k on and k off per unit time, respectively. The transition is assumed to be controlled by a feedback mechanism induced by the binding/unbinding of a given number of X-protein molecules. Transcription of messenger RNA (mRNA) from the active DNA form, and translation into protein X are assumed to occur at rates (per unit time) k m and k x , respectively. k ε is the rate constant associated with transcriptional leakage. The mRNA and protein degradations are assumed to occur by first order processes with rate constants γ m and γ x , respectively. * Author to whom correspondence should be addressed. E-mail: antonio@iim.csic.es · · · n − 1 n n + 1 · · · g n i g n n − 1 g n + 1 n g i n γ x (n + 1) γ x (n) Figure S3: Jump process representation of one protein produced in bursts, where one state n can be reached from lower states 0 ≤ i < n with different transition probability functions g n i . Equivalently, from the state n the protein number can jump to higher states i with transition probability function g i n . The degradation follows a one step process (i. e. from state n to state n − 1).

A Rates of convergence (truncation method)
Here we describe a truncation method to compute the rates of convergence. Let P : R + × N → [0, 1], be the probability of having n proteins at time t. The time evolution of P(t, n) is given by the following CME with jumps, that reads: where the transition probability g j i is proportional to the production rate of messenger RNA, so that: In order to obtain an approximation of the convergence rate of the PIDE model 3 towards the stationary state, we use the truncated form of the discrete equation (S1). Let N be the maximum possible number of proteins. Then, equation (S1) can be written in matrix form as: where the matrix M reads: with the elements of the diagonal d i being of the form: equivalently: The steady state is given by the null space of matrix M, which is spanned by the normalized eigenvector associated to the unique zero eigenvalue, as the associated eigenspace has dimension one. Actually, since the graph associated to matrix M (Fig. S3) has one trap, all the eigenvalues are negative except one (which is zero) 2 . By λ 1 , we denote the negative eigenvalue closer to zero, i.e the one with smallest absolute value.

B Mutual inhibitory gene regulatory motif in yeast
We consider the mutual repression gene network in 1,4 , where the Lacl promoter is repressed by the protein expressed by the TetR promoter and vice versa, and ATc is used to inhibit the expression of TetR. Let us define x = (x 1 , x 2 ) with x 1 and x 2 being the amounts of Lacl and TetR respectively, and A be the amount of ATc. We use the following input functions to accommodate the network to the PIDE formulation 3 : where the parameters n t = 1.56, n l = 3.35, k t = 11, k l = 264, k AT c = 0.94, and the degradation rate of the proteins γ i x = 0.002 min −1 are taken from 4 . In 1 we find C rl = C rt = 0.005 and n t m ≈ 11.5, so we consider that m ≈ 7.37. Finally we set A = 4 (because hysteresis was observed for a range of AT c between 0 and 250). We take γ i m , k i m and k i x such that The transient distributions are depicted in Fig. S6. Fig. S7 represents the corresponding marginal distribution for the same snapshots. As it is shown, the distribution at 50 h resembles an stationary distribution, since no significant differences are observed with those obtained at t=100 and even at t=150. However, comparing those distributions with the stationary distribution (see also third row in Fig. S8), we clearly conclude that the system is not at the stationary state at all. Thus, the corresponding dose-response curve at t=50 should present a transient hysteresis phenomenon.
Note that as shown in Figures S8 and S9, even snapshots taken at much longer times (e.g. 1500 h) still differ significantly from the stationary solution.
The results for t = 50 are coherent with the observation by 4 that if a trajectory starts clearly within one of the "basins of attraction" remains there for a long time. Note that the time needed to reach the stationary state might be longer than the natural timescales of relevance to the process. This is in accordance with 4 where the transitions are characterized as stochastic yet irreversible.     .