Stochastic epigenetic dynamics of gene switching
SStochastic epigenetic dynamics of gene switching
Bhaswati Bhattacharyya, ∗ Jin Wang, † and Masaki Sasai ‡ Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan Department of Chemistry, Physics and Applied Mathematics,State University of New York at Stony Brook, Stony Brook, New York, 11794-3400, USA (Dated: October 20, 2020)Epigenetic modifications of histones crucially affect the eukaryotic gene activity, while the epige-netic histone state is largely determined by the binding of specific factors such as the transcriptionfactors (TFs) to DNA. Here, the way how the TFs and the histone state are dynamically correlatedis not obvious when the TF synthesis is regulated by the histone state. This type of feedback regula-tory relations are ubiquitous in gene networks to determine cell fate in differentiation and other celltransformations. To gain insights into such dynamical feedback regulations, we theoretically ana-lyze a model of epigenetic gene switching by extending the Doi-Peliti operator formalism of reactionkinetics to the problem of coupled molecular processes. The spin-1 and spin-1/2 coherent staterepresentations are introduced to describe stochastic reactions of histones and binding/unbinding ofTF in a unified way, which provides a concise view of the effects of timescale difference among thesemolecular processes; even in the case that binding/unbinding of TF to/from DNA are adiabaticallyfast, the slow nonadiabatic histone dynamics give rise to a distinct circular flow of the probabilityflux around basins in the landscape of the gene state distribution, which leads to hysteresis in geneswitching. In contrast to the general belief that the change in the amount of TF precedes the histonestate change, the flux drives histones to be modified prior to the change in the amount of TF inthe self-regulating circuits. The flux-landscape analyses shed light on the nonlinear nonadiabaticmechanism of epigenetic cell fate decision making.
I. INTRODUCTION
Epigenetic pattern formation associated with chemicalmodifications of histones provides long-term memory ofgene regulation, which plays a critical role in cell differ-entiation, reprogramming, and oncogenesis [1, 2]. Theeffects of epigenetic modifications have been discussedtheoretically [3–12]; however, their quantitative dynam-ics still remain elusive. Statistical mechanical analyseshave suggested that histones in an array of ∼ –10 interacting nucleosomes (i.e., particles of histone-DNAcomplex) show collective changes in their modificationpattern [13–15], and such collective modifications havebeen indeed found in loops or domains of chromatin [16–18]. It has been traditionally considered that concen-tration of transcription factor (TF) largely determinesthe histone state [19]; however, the relation between TFand the collective histone state is not obvious when bothTF and the histone state are regulated with a networkof feedback loops. In particular, the mechanism of howhistone modifications are induced prior to the gene ac-tivation in developmental processes remains a mystery[20–22]. Therefore, to get physical insights into these dy-namical regulatory systems, the relation between TF andthe histone state needs to be tested by physical models.In physical modelling, particularly important is torepresent the collective histone state as a dynamicallyfluctuating variable to examine the effects of TF bind- ∗ [email protected] † [email protected] ‡ [email protected] ing/unbinding on the histone state dynamics. Here, weneed to take account of the effects of timescale differ-ence among various molecular processes; when a certainmolecular process is much faster than the other pro-cesses, the fast process can be regarded as being averagedas in quasi-equilibrium in the adiabatic approximation,while the slowest process could be regarded as station-ary in the nonadiabatic limit. The effects of adiabatic-ity/nonadiabaticity on the simple bacterial gene dynam-ics have been intensively investigated from the statisticalphysical viewpoint with theoretical [5, 23–30] and exper-imental [31, 32] methods, which revealed the large fluctu-ation of gene activity in the middle-range regime betweenadiabatic and nonadiabatic limits. However, the effectson the more complex eukaryotic genes have remainedchallenging [3, 4, 12, 33]. In the present paper, we in-vestigate the problem of adiabaticity/nonadiabaticity ineukaryotic genes by explicitly taking account of the de-gree of collective histone state transitions, which are themechanism absent from bacterial genes but play a centralrole in eukaryotic gene switching.A straightforward way to analyze the physical mod-els of gene regulation is to simulate their kinetics witha Monte-Carlo type numerical method. However, sucha calculation does not necessarily provide a global un-derstanding and physical picture directly. A clearer pic-ture would be obtainable when the stochastic kinetics isdescribed with the chemical Langevin equation, whichemerges in the high molecular number or large volumelimit as a continuous description of stochastic chemi-cal reactions. The Langevin dynamical method leads tothe combined description of probability flux and land-scape, where the landscape represents the distribution a r X i v : . [ q - b i o . M N ] O c t of states generated via nonlinear interactions and theflux reflects the nonequilibrium driving force of transi-tions among states. The combined flux-landscape anal-yses have been useful to obtain a global and physicalperspective of various complex systems [34–36]; however,in the present problem of gene switching, the continuouschanges of molecular concentrations are coupled with thediscrete processes of TF binding/unbinding and transi-tions of histone states. Therefore, to describe such cou-pled discrete and continuous dynamics, we need to con-sider the coupled multiple landscapes; the Langevin dy-namics are the motions on individual landscapes and thediscrete molecular state changes are transitions amonglandscapes [29, 30]. In this problem, even utilizing theflux-landscape method on individual landscapes, it is stilldifficult to obtain a global picture when the local stochas-tic transitions among multiple landscapes are frequent.To overcome this difficulty, we here develop a theoreticalmethod by extending the Doi-Peliti operator formalism[37–39] of chemical reaction kinetics. By introducing thespin-1 and spin-1/2 coherent-state representations, thecombined discrete and continuous description is trans-formed to a unified continuous description with expandeddimension. Then, the coupled molecular processes aredescribed as continuous dynamics on a single landscapein the higher dimensional space, which leads to a globalpicture and quantification of gene switching dynamics.Using this extended Doi-Peliti method, we show that thenonadiabatic dynamics of histone state transitions giverise to a nontrivial temporal correlation and hysteresis ineukaryotic gene switching. II. A PHYSICAL MODEL OF EUKARYOTICGENE SWITCHINGA. Self-regulating genes
Since the histones modified with the gene-repressingmarks and the histones having the gene-activating marksare not directly transformed to each other but are trans-formed through the multiple steps of the mark-deletionand the mark-addition or the replacement with the newlysynthesized unmarked histones, it is natural to describethe individual histones with three states; gene-repressing,unmarked, and gene-activating states [13, 14]. Similarly,the quantitative experimental data of transitions amongcollective histone states, which emerge through the co-operative interactions among many nucleosomes, havebeen fitted by the three-state transition models [40, 41].Therefore, in the similar way to those previous mod-els, we classify the collective histone states around apromoter site of a modeled gene into three states: thegene-activating state with histones marked as H3K9acor others ( s = 1), the gene-repressing state marked asH3K9me3 or others ( s = − s = 0). We write the transition rate from state s (cid:48) to s as r ss (cid:48) . The chromatin chain with s = 1 takes an openstructure, which facilitates access of the large-sized tran-scribing molecules to DNA to enhance the gene activity,whereas the s = − σ = 1) or unbinds fromDNA ( σ = −
1) with binding rate h and unbinding rate f . These rates are insensitive to the chromatin packingdensity when the size of TF is as small as the pioneerTFs which can diffuse into the compact chromatin space[42, 43]. Here, for simplicity, we consider that the TFis a pioneer factor as Sox2 or Oct4 in mammalian cells[43] and its binding/unbinding rates are not affected sig-nificantly by the chromatin compactness or the histonestate. However, the bound TF should recruit histonemodifier enzymes, so that binding of a single TF nu-cleates the histone state change, which is expanded andpropagates along the DNA sequence to induce the col-lective histone state change as observed in engineered[40, 41, 44] and model cells [45, 46]. Therefore, bindingof the TF modifies r ss (cid:48) as r ss (cid:48) = r ss (cid:48) + δ σ ∆ r ss (cid:48) , where δ σ is a Kronecker delta. The rate of collective changein many histones, r ss (cid:48) , should be smaller in general thanthe rate of single-molecule binding/unbinding, h or f ,which enables histones to retain the effect of TF-bindingas memory; ∆ r ss (cid:48) > s > s (cid:48) when the TF is anactivator while ∆ r ss (cid:48) > s < s (cid:48) when the TF is arepressor.As a prototypical motif of gene circuits, we consider aself-regulating gene as in Fig. 1; a dimer of the productprotein is the TF acting on the gene itself. Here, by as-suming that dimerization is much faster than the otherreactions, the TF-binding rate should be h = h p , where p = n/ Ω is the protein concentration, n is the proteincopy number in the nucleus, Ω is a typical copy numberof the protein in the nucleus when the transcription fromthe gene is active, and h is a constant. The protein pro-duction rate depends on both σ and s as g σs . Here, g is the largest and g − − is the smallest of g σs ≥ g − is largest and g − is small-est for a repressor TF. With the approximation that theprotein degradation depends on the total copy number n , the degradation rate is kn with a constant k . In thisway, the rates of the histone state transitions are affectedby the binding status of the TF, while the binding rateof the TF is regulated by the amount of the TF, and thesynthesis rate of the TF is affected by the histone state.Therefore, the histone state is a part of the feedback loopof the regulation. Because the rates of the histone statetransitions are smaller than the TF binding/unbindingrates in general, the histone state constitutes a relativelyslowly varying part in the feedback loop.Self-activating motifs discussed in the present paperare ubiquitous in cells. In mouse embryonic stem cells(mESCs), for example, the genes necessary for sustain- FIG. 1. Schematic of a self-regulating gene. Close-up viewsof chromatin at around the promoter region are illustratedat the bottom, where histones (brown) are actively marked(red) or repressively marked (green). Histones in the chro-matin are collectively marked through their mutual interac-tions. When the activating mark is dominant, chromatin hasan open structure, which enhances the gene transcription ac-tivity ( s = 1), and when the repressive mark is dominant,condensed chromatin suppresses the gene activity ( s = − s = 0). The transcription factor (TF) is a dimer of the prod-uct protein (blue oval). The bound TF recruits histone modi-fiers (black), which trigger the collective histone modificationas was observed in engineered cells [40, 41, 44] by perturb-ing the transition rates r ss (cid:48) among the histone states. Theprotein production rate g σs depends both on the TF bindingstatus, σ , and on the histone state, s . See the text for thedefinition of rates denoted on the arrows. ing pluripotency are activating each other. A Sox2-Oct4heterodimer binds on the gene loci of Sox2 and
Oct4 andactivates themselves [47, 48]. The self-activator gene inthe present paper can be regarded as a simplified modelof this
Sox2 - Oct4 system when these two genes are de-scribed as the strongly correlated loci.
B. Operator formalism of reaction kinetics
To describe the reactions in Fig. 1, we use the op-erator formalism of Doi and Peliti [37–39], which hasbeen applied to the problems of gene regulation with-out explicitly considering the histone state s [29, 33, 49–51] and to the problem of histone state without consid-ering the TF-binding status σ [15]. Now, taking intoaccount those processes having different timescales in aunified way, we consider the probability that the pro-tein copy number is n at time t , P σs ( n, t ). We de-fine a six dimensional vector ψ ( t ), whose component is ψ ( t ) σs = (cid:80) n P σs ( n, t ) | n (cid:105) . The operators a and a † areintroduced as a † | n (cid:105) = | n + 1 (cid:105) , a | n (cid:105) = n | n − (cid:105) , and [ a, a † ] = 1. Then, by assuming that all the reactionsare Markovian, the master equation of the reactions is ∂∂t ψ = − H ψ , with H being a six dimensional opera-tor, H = k ( a † a − a ) + G (1 − a † ) + J ⊗ K, (1)where is a unit matrix, G is a diagonal matrix, G σs,σ (cid:48) s (cid:48) = δ σs,σ (cid:48) s (cid:48) Ω g σs , and J and K are transition ma-trices for σ and s , respectively; J = (cid:18) − f hf − h (cid:19) ,K = − r r r − r − r − r − r − − r − , (2)with h = h ( a † a ) .We should note that H is non-Hermitian reflectingthe nonequilibrium features of the processes in gene ex-pression. From Eq. 1, we can formally write the transi-tion probability matrix P ( n f , τ | n i ,
0) between the state n = n i at t = 0 and the state n = n f at t = τ as P ( n f , τ | n i ,
0) = 1 n f ! (cid:104) n f | exp (cid:18) − (cid:90) τ dt H (cid:19) | n i (cid:105) . (3) C. Continuous dynamics in the higher dimensions
The temporal development of ψ ( t ) can be representedin a path-integral form by using the coherent-state rep-resentation, | z (cid:105) = e a † z | (cid:105) , with a complex variable z ( t ).The transition paths in the σ and s space can be repre-sented by using the spin-1/2 and spin-1 coherent stateswith spin angles θ ( t ) and α ( t ) and their conjugate vari-ables φ ( t ) and β ( t ); σ ( θ, φ ) = e iφ/ (cos θ σ + e − iφ/ (sin θ σ − s ( α, β ) = e iβ (cos α s + 2(cos α α s + e − iβ (sin α s − , (4)with σ i = ( δ i , δ i − ) T for TF binding/unbinding and s i = ( δ i , δ i , δ i − ) T for the histone degree of freedom.Considering the non-Hermiticity of H , we introduce theconjugate vectors, (cid:104) z | = (cid:104) | e az ∗ , ˜ σ = e − iφ/ σ T + e iφ/ σ T − , ˜ s = e − iβ s T + s T + e iβ s T − , (5)where z = ψe − iχ , ψ = | z | , and z ∗ = e iχ . With pairsof vectors defined in Eqs. 4 and 5, we have the identitymatrix = I z ⊗ I σ ⊗ I s as I z = 12 π (cid:90) π dψ (cid:90) π − π dχ | z (cid:105)(cid:104) z | e − ψ ,I σ = 12 π (cid:90) π sin θdθ (cid:90) π dφ ( ˜ σσ ) ,I s = 34 π (cid:90) π sin αdα (cid:90) π dβ (˜ ss ) . (6)Using Eq. 6, Eq. 3 is represented in a path-integral formas P ( n f , τ | n i , const. (cid:90) DαDβDθDφDψDχ exp (cid:18) − (cid:90) L dt (cid:19) , (7)where L is an effective “Lagrangian”, L = L χψ + L φθ + L βα , and L χψ = iχ dψdt + Ω(1 − e iχ ) (cid:104) g − ( θ ) sin α g ( θ ) sin α α g ( θ ) cos α (cid:105) +(1 − e − iχ ) kψ, L φθ = i φ ddt cos θ + h Ω − ψ (1 − e iφ ) sin θ f (1 − e − iφ ) cos θ , L βα = iβ ddt cos α +(1 − e iβ ) (cid:104) r − ( θ ) sin α r ( θ ) sin α α (cid:105) +(1 − e − iβ ) (cid:104) r − ( θ ) sin α α r ( θ ) cos α (cid:105) , (8)with g s ( θ ) = g s cos θ + g − s sin θ and r ss (cid:48) ( θ ) = r ss (cid:48) +∆ r ss (cid:48) cos θ .Then, by retaining up to the 2nd order terms of χ , φ ,and β in Eq. 7 (i.e., using the saddle-point approxima-tion) in a similar way to the method in [29], we obtain theLangevin equations describing fluctuations in the proteinconcentration p = ψ/ Ω, the TF binding status ξ = cos θ ,and the histone state ζ = cos α ; dpdt = F + p ( θ, α ) − kp + η p ,dξdt = 2 h p ξ s − f ξ c + η ξ ,dζdt = F + ζ ( θ, α ) − F − ζ ( θ, α ) + η ζ , (9)where F + p = ξ c (cid:2) g ζ + 2 g ζ c ζ s + g − ζ (cid:3) + ξ s (cid:2) g − ζ + 2 g − ζ c ζ s + g − − ζ (cid:3) ,F + ζ = 2 ( r + ξ c ∆ r ) ζ c ζ s + ( r − + ξ c ∆ r − ) ζ ,F − ζ = ( r + ξ c ∆ r ) ζ + 2 ( r − + ξ c ∆ r − ) ζ c ζ s , (10) with ξ c = cos ( θ/ ξ s = sin ( θ/ ζ c = cos ( α/ ζ s = sin ( α/ η p , η ξ , and η ζ aremutually independent Gaussian noises with (cid:104) η x (cid:105) = 0, (cid:104) η x ( t ) η x ( t (cid:48) ) (cid:105) = B x δ ( t − t (cid:48) ) for x = p , ξ , or ζ as B p = (cid:0) F + p ( θ, α ) + kp (cid:1) / Ω ,B ξ = 4 h p ξ s + 4 f ξ c ,B ζ = F + ζ ( θ, α ) + F − ζ ( θ, α ) . (11)We should note that p = n/ Ω was an almost contin-uous variable for a large Ω. Hence, the original cou-pled dynamics of a nearly continuous variable p and thediscrete variables, σ and s , were transformed into theBrownian dynamics in the 3-dimensional (3D) space ofthe continuous variables, p , ξ , and ζ , with TF bound( ξ = 1)/unbound ( ξ = −
1) and histone state activating( ζ = 1)/neutral ( ζ = 0)/repressing ( ζ = − p , ξ , and ζ outside the origi-nally defined range of values, p ≥ − ≤ ξ ≤
1, and − ≤ ζ ≤
1. In order to reduce this out-of-range fluctua-tions, we add soft-wall forces, w p , w ξ , and w ζ , to Eq. 9 innumerical simulations as dpdt = F + p ( θ, α ) − kp + w p + η p , dξdt = 2 h p ξ s − f ξ c + w ξ + η ξ , and dζdt = F + ζ ( θ, α ) − F − ζ ( θ, α ) + w ζ + η ζ with w p = (cid:26) − c ( p − . for p ≤ .
10 otherwise ,w ξ = − c ( ξ + 0 . for ξ ≤ − . − c ( ξ − . for ξ ≥ .
90 otherwise ,w ζ = − c ( ζ + 0 . for ζ ≤ − . − c ( ζ − . for ζ ≥ .
90 otherwise , (12)with a constant c > D. Adiabatic approximation of TFbinding/unbinding
In a single-molecule tracking experiment of the TFin mESCs, the observed timescale of binding/unbindingwas ∼ / min [43], whereas the observed degradationrate k of TF was ∼ . / hour [52, 53]. Though a sin-gle histone can be replaced in ∼ hour, many histonescollectively change with the rate ∼ / day [40]. There-fore, the estimated ratios are h/k ∼ f /k = O (10 )and r ss (cid:48) /k = O (1) ∼ O (10 − ). With such large h and f , TF-binding/unbinding reactions can be treatedas adiabatic : TF-binding/unbinding are regarded asin quasi-equilibrium as (cid:104) dξ/dt (cid:105) = 0, leading to ξ s =1 / (cid:2) ( h /f ) p + 1 (cid:3) . With this adiabatic approximation,the 3D calculation in Eq. 9 for ( p, ξ, ζ ) is reduced to the2D calculation for ( p, ζ ). On the other hand, the rateof collective histone change is small; hence, the histonedynamics remains nonadiabatic .Validity of the adiabatic treatment of TF bind-ing/unbinding can be checked by comparing the simu-lated results of the present model with the experimentaldata of Hathaway et al. [40]. Hathaway et al. developeda technique to forcibly bind a chromo-shadow domain ofHP1 (csHP1) to a DNA site near the promoter of Oct4 in mESCs. The bound csHP1 nucleated the repressivelymarked histones and the histones around the promoterregion were collectively transformed to the repressivestate in several days. In Fig. 2, this observation is com-pared with the calculated results obtained by applyingthe adiabatic approximation of TF binding/unbindingto Eq. 9. We assumed that csHP1 binding intensivelysuppresses the acetylation and other activating modifi-cations of histones, so as to reduce the corresponding r ss (cid:48) /k from O (10 − ) to O (10 − ). The simulated re-sults reproduce the relaxation of the collective histonestate toward the repressive state after csHP1 binding.Here, the experimental histone state was quantified as ζ ( t ) = H active ( t ) /H active (0) − H repress ( t ) /H repress (day5),where H active ( t ) and H repress ( t ) are the observed frac-tions of histones marked as H3K27ac and H3K9me3, re-spectively, in the ∼
10 kb region around the
Oct4 gene.
Oct4 was transcriptionally active for t <
0, but its his-tone state was modified with the bound csHP1 for t ≥ III. LANDSCAPES, CIRCULAR FLUXES, ANDTEMPORAL CORRELATION
Either with adiabatic or nonadiabatic treatment of TF-binding/unbinding kinetics, the flux-landscape approach[29, 30, 34, 36] provides a concise view of Eq. 9. Withthe adiabatic approximation for example, the landscape, U ( p, ζ ), is obtained from the calculated stationary proba-bility density distribution as U ( p, ζ ) = − ln P ( p, ζ ). Theprobability flux J is obtained as a 2D vector field in theadiabatic TF-binding/unbinding case from the Fokker-Planck equation corresponding to Eq. 9, ∂∂t P ( p, ζ, t ) = −∇ · J ( p, ζ, t ), where ∇ = ( ∂ p , ∂ ζ ) and J = ( J p , J ζ )with J p = [ F + p ( θ, α ) − kp + w p ] P − ∂∂p [ B p P ] and J ζ = [ F + ζ ( θ, α ) − F − ζ ( θ, α ) + w ζ ] P − ∂∂ζ [ B ζ P ]. Of note,even in the stationary state, J can be nonzero when ∇ · J = 0. This divergence-less circulating flux is a hall-mark of the breaking of the detailed balance reflectingthe biased thermal/chemical energy flow such as the nu-cleotide consumption in protein synthesis and the heatdissipation [29, 35, 36].Fig. 3 shows the calculated U and J at the stationarystate in the case of adiabatic TF-binding/unbinding. Foran activator TF, when the TF-binding affinity is small(small h /f ) (Fig. 3a), U has a single basin at p ≈ ζ ≈ − ∼ h /f ) (Fig. 3b), U has two FIG. 2. Relaxation of a self-activating gene to the repres-sive state. Eq. 9 was numerically calculated by applying theadiabatic approximation of TF binding/unbinding. The pro-tein concentration p ( t ) (black line) and the histone state ζ ( t )(red line) were derived by averaging 1000 simulated trajecto-ries. See text for the experimental estimation of ζ ( t ) (greendots) from the data of [40]. k was set to k = 2 day − ,and the other rates were defined in units of k ; for t < r = r − = 0 . r = r − = 0 .
1, and ∆ r = ∆ r − = 1,and for t ≥ r = r − = r = 0 . r − = 0 . r = ∆ r − = 0 .
01. For both t < t ≥ r = ∆ r − = 0, h /f = 200, g = 1, g − = g = 0 . g − = g − = g − − = 0, Ω = 100, and c = 20. coexisting basins in the off state and at p ≈ . ζ ≈ h /f ) (Fig. 3c), a dominant basin is found at theon state. Thus, the binding affinity of the activator TFdetermines the distribution of stable states and works asa switch of gene states. When the TF is a repressor, wefind a single basin at intermediate p and ζ for a widerange of binding affinity (Fig. 3d).In all the cases shown in Fig. 3, we find a flux J glob-ally circulating around a basin or traversing betweenbasins though the deterministic part of Eq. 9 is non-oscillatory; the oscillatory feature of the flow becomesevident through the stochastic on-off switching fluctua-tions as a stochastic resonance effect. When the epige-netic effect is absent, the flux is diminished in the adi-abatic limit [29, 30]. However, here with nonadiabaticepigenetic histone dynamics, the circular flux is signifi-cant even in the limit of adiabatic TF-binding/unbindingbecause the timescales in p and ζ are near to each other,so that the two processes are non-separable. The evidentcirculating flux suggests a temporal correlation betweenthe histone modification and the gene activity change.By following the flux direction along the off-to-on (theon-to-off) path, the histone state first tends to becomeactivating (repressing) followed by increase (decrease) inprotein concentration.This temporal correlation is confirmed by calculating FIG. 3. Landscape U ( p, ζ ) and probability flux J ( p, ζ ) cal-culated on the 2D plane of the protein concentration p andthe histone state ζ in the adiabatic approximation of TF-binding/unbinding. U is shown with contours and J withyellow arrows. (a–c) Self-activating and (d) self-repressinggenes. (a) h /f = 2, (b) h /f = 20, (c) h /f = 200, and(d) h /f = 2. The rate parameters are scaled in units of k by setting k = 1; g = 1, g − = g = 0 . g − = g − = g − − = 0, ∆ r = ∆ r − = 1, ∆ r = ∆ r − = 0 for (a)–(c)and g − = 1, g = g − = 0 . g = g − = g − − = 0,∆ r = ∆ r − = 0, and ∆ r = ∆ r − = 1 for (d). For(a)–(d), Ω = 100, c = 20, and r = r − = r = r − = 1. the normalized difference between the two-time cross cor-relations, A ( t ) =[ (cid:104) δζ ( τ ) δp ( τ + t ) (cid:105)− (cid:104) δp ( τ ) δζ ( τ + t ) (cid:105) ] / | (cid:104) δζ ( τ ) δp ( τ ) (cid:105) | , (13)with δζ = ζ − (cid:104) ζ (cid:105) and δp = p − (cid:104) p (cid:105) , where (cid:104)· · · (cid:105) isthe average over τ and the calculated trajectories. Wecan write A ( t ) ∝ (cid:104) det [ q ( τ ) , q ( τ + t )] (cid:105) with a 2D vector q ( τ ) = ( ζ ( τ ) , p ( τ )). A positive value of A ( t ) at t > ζ tends to increase(decrease) p at a later time t . For activator (Fig. 4a)and repressor (Fig. 4b) TFs, A ( t ) is plotted for variousvalues of r ss (cid:48) /k = u , showing that A ( t ) has a positive-valued peak at t u = 1 /u ∼ . /u . We find that the peakis evident even when the histone dynamics are as slow as u <
1, which indicates that the prior change in the his-tone state to the gene activity is not owing to the fasterrate of reactions in histones but is due to the circular fluxgenerated by the nonadiabatic dynamics of histones. Thedelayed influence of ζ on p should lead to the differenton-to-off and off-to-on paths, inducing hysteresis in theswitching dynamics.The hysteresis is shown by calculating the optimalpaths of transitions. An optimal path is obtained by firstsetting its start and end points and then minimizing theeffective action in the path-integral formalism of kinet-ics connecting those points [15, 33, 36, 54–57]. Figs. 4c FIG. 4. Temporal correlation and optimal paths in self-regulating genes. (a, c) Self-activating and (b, d) self-repressing genes. (a, b) The normalized difference betweentwo-time cross correlations, A ( t ), is plotted as a function of t in units of 1 /k . (c, d) The optimal paths for the on-to-off (green) and off-to-on (purple) directions are superposedon the flux J on the 2D plane of the protein concentra-tion p and the histone state ζ . Calculated using the adia-batic approximation of TF binding/unbinding. The rates ofthe histone state change are scaled by the parameter u as r = r − = r = r − = u with ∆ r = ∆ r − = u in(a, c) or ∆ r = ∆ r − = u in (b, d). In (a) and (b), A ( t )is plotted with u = 1 . u = 1 (red), and u = 0 . u = 1. The other parameters of (a,c) are same as in Fig. 3b and those of (b, d) are same as inFig. 3d. and 4d show the paths calculated by the simulated an-nealing of the action with the algorithm of [36]. Thusobtained off-to-on and on-to-off paths are indeed differ-ent from each other, both of which are consistent in theirroute orientations with the circulating probability flux.Because equilibrium kinetic paths should pass throughthe same saddle point of the landscape in both directionswithout showing the hysteresis, the calculated hysteresisis a manifestation of the nonequilibrium feature, and thearea formed by the loop of paths gives a measure of theheat dissipation [58].When the TF binding/unbinding is slow, we need tosolve Eq. 9 nonadiabatically. Slow nonadiabatic bind-ing/unbinding were examined recently by tuning thebinding rate in bacterial cells experimentally [31, 32].Fig. 5a shows the landscape and flux in such a nona-diabatic binding/unbinding case with the intermediatebinding affinity of the activator TF. We find two basinsfor the on and off states and a distinct circular flux be-tween them. The qualitative feature is same as in theadiabatic TF-binding/unbinding case; however, here, wefind the correlated binding/unbinding behavior with thehistone-state change, which should generate hysteresis inthe 3D space. This hysteresis can be found in the calcu- FIG. 5. Landscape, probability flux, and optimal paths of aself-activating gene calculated in the 3D space of the proteinconcentration p , the TF binding status ξ , and the histonestate ζ , in the nonadiabatic kinetics. (a) Flux J is superposedon the landscape U with 1 < U ≤ − . < U ≤ U ≤ − . p ≈ ξ ≈ −
1, and ζ ≈ ∼ − p ≈ ξ ≈ ζ ≈
1. (b) Flux and optimal paths of on-to-off(green) and off-to-on (purple) directions. The rate parametersare scaled in units of k with h = 2 and f = 0 .
1. The otherparameters are same as in Fig. 3b. lated optimal paths in the 3D space (Fig. 5b).
IV. DISCUSSION
The present flux-landscape analyses provided a newview that the nonadiabatic circular flux generates non-trivial temporal correlation, hysteresis, and dissipationin eukaryotic gene switching. These analyses provide aclue to resolving the ‘chicken-and-egg’ argument on thecausality between the histone state and TF binding. Thelandscape analyses showed that the stability of each hi-stone state is determined by the binding affinity of TF,which is in accord with the general belief that the specificTF binding is the cause and the histone-state change isthe result [19]. However, unexpectedly, the histone statetends to change in self-regulating gene circuits prior tothe change in the amount of TF, which induces hysteresisin the switching dynamics; the histone-state fluctuationcan be a trigger for switching the feedback loop. Thistemporal correlation should be tested experimentally byanalyzing A ( t ) in single-cell observations. This type ofexperiments should be possible when the amount of theproduct TF is measured by a co-expressing fluorescentprotein and the live-cell histone state is monitored si-multaneously by the technique of fluorescently labeledspecific antigen binding [21, 59]. The flux structure andtimescales should be also tested experimentally by exam-ining the irreversibility in temporal correlations [60].A possible test of the present model is to quantita-tively monitor the response of somatic cells to the incor-poration of exogenous genes such as Yamanaka factors,which include Sox2 and
Oct4 [61]. By introducing Ya-manaka factors, the differentiated somatic cells can turn into the pluripotent cells. This reprogramming of cells
FIG. 6. Landscape U ( p, ζ ) and probability flux J ( p, ζ ) withthe exogenous TF contribution p ex . (a) p ex = 0 . p ex = 0 .
6. Calculated on the 2D plane of the protein concen-tration p and the histone state ζ in the adiabatic approxima-tion of TF-binding/unbinding. The parameters are the sameas in Fig. 3a. may start from the binding of exogenous pioneer factors,Sox2 and Oct4, to the loci of endogenous genes, Sox2 and
Oct4 . By writing the concentration of exogenousSox2 and Oct4 proteins as p ex , the total concentrationof a Sox2-Oct4 heterodimer TF should be approximatelyproportional to ( p + p ex ) in the present model; hence,the TF binding rate should become h = h ( p + p ex ) .Shown in Fig. 6 are the landscapes and fluxes for thesmall (Fig. 6a) and large (Fig. 6b) values of p ex with theparameters of Fig. 3a, with which the landscape is dom-inated by the off-state when p ex = 0. We find that thelandscape is shifted upon introduction of exogenous fac-tors with p ex > Nanog induces a large fluctuation in mESCs [3]. Itis intriguing to develop a method to explain the degreeof freedom of super-enhancer formation/dissolution andthe associated large-scale chromatin structural change byextending the present scheme. Thus, the flux-landscapeapproach should provide physical insights into variousproblems of nonadiabatic stochastic switching.
ACKNOWLEDGEMENTS
B. B. and M. S. thank Dr. S. S. Ashwin for fruitfuldiscussions. This work was supported by JST-CRESTGrant JPMJCR15G2, the Riken Pioneering Project, andJSPS-KAKENHI Grants JP19H01860, 19H05258 and20H05530. J. W. thanks NSF-PHY-76066 for supports. [1] R. Bonasio, S. Tu, and D. Reinberg, Science , 612(2010).[2] W. A. Flavahan, E. Gaskell, and B. E. Bernstein, Science , eaal2380 (2017).[3] M. Sasai, Y. Kawabata, K. Makishi, K. Itoh, and T. P.Terada, PLoS Comput Biol , e1003380 (2013).[4] S. S. Ashwin and M. Sasai, Sci Rep , 16746 (2015).[5] H. Feng, B. Han, and J. Wang, J Phys Chem B ,1254 (2011).[6] H. Feng and J. Wang, Sci Rep , 550 (2012).[7] C. Li and J. Wang, J R Soc Interface , 20130787 (2013).[8] C. C. Chen, S. Xiao, D. Xie, X. Cao, C. X. Song,T. Wang, C. He, and S. Zhong, PLoS Comput Biol ,e1003367 (2013).[9] C. Chen and J. Wang, Sci Rep , 20679 (2016).[10] N. Folguera-Blasco, R. P´erez-Carrasco, E. Cuy´as, J. A.Menendez, and T. Alarc`on, PLoS Comput Biol ,e1006592 (2019).[11] X. J. Tian, H. Zhang, J. Sannerud, and J. Xing, ProcNatl Acad Sci U S A , E2889 (2016).[12] C. Yu, Q. Liu, C. Chen, J. Yu, and J. Wang, PhysicalBiology , 051003 (2019).[13] I. B. Dodd, M. A. Micheelsen, K. Sneppen, and G. Thon,Cell , 813 (2007).[14] H. Zhang, X. J. Tian, A. Mukhopadhyay, K. S. Kim, andJ. Xing, Phys Rev Lett , 068101 (2014).[15] A. Sood and B. Zhang, Phys Rev E , 062409 (2020).[16] J. R. Dixon, S. Selvaraj, F. Yue, A. Kim, Y. Li, Y. Shen,M. Hu, J. S. Liu, and B. Ren, Nature , 376 (2012).[17] S. S. Rao, M. H. Huntley, N. C. Durand, E. K. Sta-menova, I. D. Bochkov, J. T. Robinson, A. L. Sanborn,I. Machol, A. D. Omer, E. S. Lander, and E. Lieber-man Aiden, Cell , 1665 (2014).[18] A. N. Boettiger, B. Bintu, J. R. Moffitt, S. Wang,B. J. Beliveau, G. Fudenberg, M. Imakaev, L. A. Mirny,C. Wu, and X. Zhuang, Nature , 418 (2016).[19] M. Ptashne, Proc Natl Acad Sci U S A , 7101 (2013).[20] M. Samata, A. Alexiadis, G. Richard, P. Georgiev,J. Nuebler, T. Tanvi Kulkarni, G. Renschler, M. F. Basili-cata, F. L. Zenk, M. Shvedunova, G. Semplicio, L. Mirny,N. Iovino, and A. Akhtar, Cell , 127 (2020).[21] Y. Sato, L. Hilbert, H. Oda, Y. Wan, J. M. Heddleston,T.-L. Chew, V. Zaburdaev, P. Keller, T. Lionnet, N. Vas-tenhouw, and H. Kimura, Development , dev179127(2019).[22] B. Zhang, X. Wu, W. Zhang, W. Shen, Q. Sun, K. Liu,Y. Zhang, Q. Wang, Y. Li, A. Meng, and W. Xie, MolCell , 673 (2018).[23] J. E. M. Hornos, D. Schultz, G. C. P. Innocentini,J. Wang, A. M. Walczak, J. N. Onuchic, and P. G.Wolynes, Phys Rev E , 051907 (2005).[24] A. M. Walczak, J. N. Onuchic, and P. G. Wolynes, ProcNatl Acad Sci U S A , 18926 (2005).[25] D. Schultz, E. B. Jacob, J. N. Onuchic, and P. G.Wolynes, Proc Natl Acad Sci U S A , 17582 (2007).[26] M. Yoda, T. Ushikubo, W. Inoue, and M. Sasai, J ChemPhys , 115101 (2007).[27] Y. Okabe, Y. Yagi, and M. Sasai, J Chem Phys ,105107 (2007).[28] P.-Z. Shi and H. Qian, J Chem Phys , 065104 (2011). [29] K. Zhang, M. Sasai, and J. Wang, Proc Natl Acad SciU S A , 14930 (2013).[30] C. Chen, K. Zhang, H. Feng, M. Sasai, and J. Wang,Phys Chem Chem Phys , 29036 (2015).[31] Z. Jiang, L. Tian, X. Fang, K. Zhang, Q. Liu, Q. Dong,E. Wang, and J. Wang, BMC Biol , 49 (2019).[32] X. Fang, Q. Liu, C. Bohrer, Z. Hensel, W. Han, J. Wang,and J. Xiao, Nat Commun , 2787 (2018).[33] B. Zhang and P. G. Wolynes, Proc Natl Acad Sci U S A , 10185 (2014).[34] X. Fang, K. Kruse, T. Lu, and J. Wang, Rev Mod Phys , 045004 (2019).[35] J. Wang, L. Xu, and E. Wang, Proc Natl Acad Sci U SA , 12271 (2008).[36] J. Wang, K. Zhang, L. Xu, and E. Wang, Proc NatlAcad Sci U S A , 8257 (2011).[37] M. Doi, J Phys A , 1465 (1976).[38] L. Peliti, J Physique , 1469 (1985).[39] D. C. Mattis and M. L. Glasser, Rev Mod Phys , 979(1998).[40] N. A. Hathaway, O. Bell, C. Hodges, E. L. Miller, D. S.Neel, and G. R. Crabtree, Cell , 1447 (2012).[41] L. Bintu, J. Yong, Y. E. Antebi, K. McCue, Y. Kazuki,N. Uno, M. Oshimura, and M. B. Elowitz, Science ,720 (2016).[42] A. Soufi, G. Donahue, and K. S. Zaret, Cell , 994(2012).[43] J. Chen, Z. Zhang, L. Li, B. C. Chen, A. Revyakin,B. Hajj, W. Legant, M. Dahan, T. Lionnet, E. Betzig,R. Tjian, and Z. Liu, Cell , 1274 (2014).[44] M. Park, N. Patel, A. J. Keung, and A. S. Khalil, Cell , 227 (2019).[45] I. M. Hall, G. D. Shankaranarayana, K. Noma, N. Ayoub,A. Cohen, and S. I. S. Grewal, Science , 2232 (2002).[46] L. N. Rusch´e, A. L. Kirchmaier, and J. Rine, Mol BiolCell , 2207 (2002).[47] S. Okumura-Nakanishi, M. Saito, H. Niwa, andF. Ishikawa, J Biol Chem , 5307 (2005).[48] S. Masui, Y. Nakatake, Y. Toyooka, D. Shimosato,R. Yagi, K. Takahashi, H. Okochi, A. Okuda, R. Ma-toba, A. A. Sharov, et al. , Nat Cell Biol , 625 (2007).[49] M. Sasai and P. G. Wolynes, Proc Natl Acad Sci U S A , 2374 (2003).[50] A. M. Walczak, M. Sasai, and P. G. Wolynes, BiophysJ , 828 (2005).[51] J. Ohkubo, J Chem Phys , 044108 (2008).[52] J. L. Chew, Y. H. Loh, W. Zhang, X. Chen, W. L. Tam,L.-S. Yeap, P. Li, Y.-S. Ang, B. Lim, P. Robson, andH.-H. Ng, Mol Cell Biol , 6031 (2005).[53] M. Thomson, S. J. Liu, L. N. Zou, Z. Smith, A. Meissner,and S. Ramanathan, Cell , 875 (2011).[54] E. Aurell and K. Sneppen, Phys Rev Lett , 048101(2002).[55] D. M. Roma, R. A. Flanagan, A. E. Ruckenstein, A. M.Sengupta, and R. Mukhopadhyay, Phys Rev E ,0119021 (2005).[56] J. Wang, K. Zhang, and E. Wang, J Chem Phys ,125103 (2010).[57] H. Feng, K. Zhang, and J. Wang, Chem Sci , 3761(2014).[58] H. Feng and J. Wang, J Chem Phys , 234511 (2011). [59] Y. Hayashi-Takanaka, K. Yamagata, T. Wakayama, T. J.Stasevich, T. Kainuma, T. Tsurimoto, M. Tachibana,Y. Shinkai, H. Kurumizaka, N. Nozaki, et al. , Nucl AcidsRes , 6475 (2011). [60] Q. Liu and J. Wang, Proc Natl Acad Sci U S A , 923(2020).[61] K. Takahashi and S. Yamanaka, Cell126