Intrinsic phenotypic stability of a bi-stable auto regulatory gene

Even under homogenous conditions clonal cells can assume different distinct states for generations to follow, also known as epigenetic inheritance. Such long periods of different phenotypic states can be formed due to the existence of more than one stable state in the molecule concentration, where the different states are explored through molecular fluctuations. By formulating a single reaction variable representing the birth and death of molecules, including transcription, translation and decay, we calculate the escape time from the phenotypic states attained from autocatalytic synthesis through a Fokker- Planck formulation and integration of an effective pseudo-potential. We calculate the stability of the phenotypic states both for cooperative binding feedback and dimer binding feedback, resulting in non-linear decay.


Cooperative binding feedback
Below we have listed the biochemical steps leading to positive autoregulation where two proteins bind cooperatively on two sites located at the promoter to assist RNA polymerase to bind and initiate transcription.

Transcription,
M → M + 1 , a 1 = ΩF(p) , ν ν ν 1 = 1 0 Translation, P → P + 1 , a 2 = k p M , ν ν ν 2 = 0 1 mRNA decay, M → M − 1 , a 3 = γ m M , ν ν ν 3 = −1 0 Protein decay, P → P − 1 , a 4 = γ p P , ν ν ν 4 = 0 −1 The function F(p) = k m,0 + k m,a p 2 K 2 + p 2 (S1) The transcription rate of the free promoter is k m,0 and the transcription rate when two proteins are present to assist the RNA polymerase to bind is k m,0 + k m,a . The function F(p) is the average occupancy of the protein complex on the promoter as a function of the protein concentration p = P/Ω. We have assumed that the cooperative binding is strong and that both proteins are needed to assist the RNA polymerase to bind on the promoter, otherwise a linear term in p should be present both in the denominator and numerator of F(p). K is the dissociation constant of the protein complex to its binding site on the promoter and gives the inverse strength of the binding.

Dimer binding feedback
Below we have listed the biochemical steps giving positive autoregulation by protein dimers binding on the promoter. Two proteins form a dimer in solution and can then bind at the promoter to assist RNA polymerase to bind and initiate transcription.
where k m,0 is transcription rate of the free promoter and the transcription rate when a protein dimer is present to assist the RNA polymerase to bind is k m,0 + k m,a . The dissociation constant K 2 gives the (inverse) binding strength of dimers on the promoter and the function F(p 2 ) is the average occupancy of dimers on the promoter as a function of the number of dimers in solution.
Most calculations and derivations that will follow are performed for dimer mediated feedback since the results for the cooperative promoter feedback can be obtained as a special case by the right choice of parameters. Stochastic simulations of the above two type of feedback systems, cooperative binding feedback and dimer binding feedback, are implemented as a reference for the analytical derivations.

Large molecule number limit calculations
In the large molecule limit, Ω → ∞, the concentrations (m, p, p 2 ) = (M, P, P 2 )/Ω become continuous and reactions occur continuous in time. The concentrations of the molecules in equation S2 will in the large molecule limit evolve according to the following set of ordinary differential equations (ODE) When mRNA turnover is fast compared to other reactions, dm dt = 0, we have that m = F(p 2 )/γ m . Then we have from equation S3 for the total concentration p T = p + 2p 2 , If monomers and dimers are (rapidly) partitioned according to the dissociation constant K 1 = k −1 /k 1 , then p 2 = p 2 /K 1 . We can now formulate a closed form differential equation for any of the concentrations p, p 2 or p T since each can be expressed 2/7 from the other two. For the bifurcation analysis that follows it is most convenient to write the equation in terms of a scaled monomer concentration ρ = p K , with K 2 = K 1 K 2 and scaled time τ = γ p t The subscript ρ and ρ T for the function Φ are denoting whether the function generates the dynamics of the monomer concentration ρ or the total concentration ρ T for dimer binding feedback. The two descriptions are equivalent since all concentrations are uniquely determined, e.g. knowing ρ uniquely gives ρ T and ρ 2 . Similar calculations as above for the cooperative binding feedback gives that We can see that this result is also obtained by letting K 1 → ∞ (i.e. no dimers) which gives α = 0 in equation S5. Note that for the case of dimers mediating the feedback -no dimers would result in no feedback. However, the operation can be performed to obtain the result for cooperative binding feedback, even though feedback is then physically mediated by monomers. We will use the result later as a quick route of deriving results for cooperative binding feedback.

Stability and bifurcation analysis
Assume that we have a potential function U(x) such that a change in the concentration x is driven by the drift term, Φ(x) = − dU dx . For the dimer binding feedback the drift (as a function on monomer proteins) equals the RHS of equation S5, Φ ρ.DIM (ρ) or equivalently Φ ρ T .DIM (ρ). For cooperative binding feedback the drift equals the RHS of equation S6, Φ CP (ρ). The observed concentration will change until the system is resting at a minima in U(x) where its derivative is zero and its second derivative positive. This is called a stable steady state. An unstable steady state is where U(x) is having a maximum. For the rate parameters where the system displays bi-stability there are always two stable stationary solutions separated by an unstable one and the region of bi-stability is given by the interval (ρ * ,1 , ρ * ,2 ) where both endpoints are solutions of Φ(ρ * ) = 0 and Φ (ρ * ) = 0. The region of bi-stability is given by where ρ * are the monomer concentrations separating the bi-stable region. The region of bi-stability in (S a , S 0 ) space can for each choice of α be visualised by plotting the parametric curves of S 0 and S a against ρ * . Moreover, it does not matter if we express the RHS of equation S7 in terms of monomer, dimer or total concentration -the region of bi-stability in (S a , S 0 ) space is the same. The bifurcation for the cooperative binding feedback is obtained by letting α = 0.

Finite molecule number calculations and intrinsic stability
In the large molecule limit, where reactions are assumed to become continuous, the concentrations will converge to a (stable) stationary solution of the differential equations S3 (and equation S5) and stay there, unless the rate constants will change and another stationary solution is assumed. However, since the number of molecules are finite, the time between events of synthesis and decay will not be continuous but random and exponentially distributed with mean following the reaction propensities.
The deviation from the stationary number of molecules that each reaction results in are also related to the total number of molecules since the smallest state jump is of the order of one molecule. Thus, for all systems consisting of a finite number of molecules we need to have a description of the noise and deviation from the stationary concentration. Such a description can be expressed differently, e.g. via a chemical master equation describing the probability of the different chemical states in time, or (as here) translating the chemical reaction system to a stochastic differential equation (SDE). We will use the former route since we are interested in calculating the escape rate from the stable states to compute their intrinsic stability. Moreover, we can express the stochastic properties of the system using any of the concentrations ρ, ρ 2 or ρ T as long as both the drift and noise term are describing the noise of the same quantity. The part describing the non-stochastic contribution to the evolution of the total concentration we already have, which is given by the right hand side of equation S5 (with α = 0 for cooperative binding 3/7 feedback). The stochastic contribution to the SDE we have to calculate, but before we do so we start with formulating the SDE and its corresponding Fokker-Planck equation (for the total protein concentration), Where Π(ρ T , τ) is the probability distribution of the scaled total protein concentration, the drift term Φ(ρ T ) the right hand side of equation , and W (τ) is a Wiener Process. Now we need to quantify the diffusive noise term θ (ρ T ) = B 2 (ρ T ) 2 . For a multivariate Ornstein-Uhlenbeck process defined by the stochastic differential equation dx(t) = −A A Ax(t)dt + B B BdW (t), where A A A and B B B are constant matrices, the covariance matrix C C C at stationary conditions satisfies a fluctuation-dissipation relation If x is sufficiently close to a stationary point x s we can for a more general system linearise the equations, e.g. S3, such that the fluctuation-dissipation relation S9 (approximately) holds. The matrix A A A is then the Jacobian matrix evaluated at the stationary point x s , and if we express the reactions in the two coordinates mRNA, m, and total protein, p T , we obtain where f m = F(p 2 ) − γ m m and f p = k p m − (γ p p + 2γ p 2 p 2 ) which gives that where β = The covariance matrix C C C can now be calculated from the fluctuation dissipation equation S9, which gives the variance of the total concentration of proteins at a stationary point as Now, assuming that mRNA turnover and protein dimerization are fast compared to the other time scales in the system we will collapse the reaction system listed in equations S2 with only synthesis and decay of proteins and its impact on the total number of proteins. One species, P T , and three reactions, where q is number of total proteins added, i.e. the stoichiometry, per synthesis is now including both transcription and translation. We can think of q as a parameter that accounts for the fact that the intrinsic stochasticity of protein synthesis propagates from transcription to translation and, if mRNA is short lived and translation fast, the multiple proteins translated per mRNA can be viewed as a stoichiometric effect on the synthesis events appearing with the rate of transcription initiation. If we want to collapse the reactions of equation S2 with the reduced set of reactions listed in S14 and recover the stochastic part we should 4/7 match the covariance of the total concentration of proteins. For a scalar system the fluctuation dissipation relation, equation S9, becomes an algebraic equation and the variance is given by S15) and the diffusive component Now, setting C = C 22 , gives q as When γ m γ p (1 + 2β ∂ p 2 ∂ p T ) we obtain q = 1 + 2 k p γ m . For the stochastic differential equation and its Fokker-Planck equation S8 describing the scaled total concentration ρ T = p T /K we obtain the noise term B(ρ T ) by combining equations S16 and S17 Now, the term B(ρ T ) (and θ (ρ T )) is strictly speaking only valid at a stable stationary point, (p, p 2 , p T ) S i.e the two stable points of a bi-stable system, since we have done the calculations under such assumptions. However, we approximate B(ρ T ) to have the same dependence on ρ T for all non-stationary points as well.
For the dimer binding feedback the SDE and its Focker-Planck equation can be expressed in closed form in terms of the total concentration, since ρ 2 = K For the cooperative binding feedback we have the drift Φ(ρ) is given by RHS of equation S6. The diffusion term θ (ρ) in the Fokker-Planck equation S8 is obtained by letting K 1 → ∞, giving ρ T = ρ and ρ 2 = 0 in equation S18

Escape time calculation
Using Eq. S8 for probability distribution of the scaled total protein concentration, the probability of this variable that remains in region R at time τ, G (ρ T , τ) = R Π(x, τ | ρ T , 0)dx, obeys the equation integrating Eq. S21 over time interval (0, ∞) gives the differential equation governing escape time, gives the escape time (in real time after division by γ p , since t = τ/γ p ), and the quasi-potential for dimer binding feedback is whereρ T = ρ T + K 1 8K , δ 1 =