Effects of Collective Histone State Dynamics on Epigenetic Landscape and Kinetics of Cell Reprogramming

Cell reprogramming is a process of transitions from differentiated to pluripotent cell states via transient intermediate states. Within the epigenetic landscape framework, such a process is regarded as a sequence of transitions among basins on the landscape; therefore, theoretical construction of a model landscape which exhibits experimentally consistent dynamics can provide clues to understanding epigenetic mechanism of reprogramming. We propose a minimal gene-network model of the landscape, in which each gene is regulated by an integrated mechanism of transcription-factor binding/unbinding and the collective chemical modification of histones. We show that the slow collective variation of many histones around each gene locus alters topology of the landscape and significantly affects transition dynamics between basins. Differentiation and reprogramming follow different transition pathways on the calculated landscape, which should be verified experimentally via single-cell pursuit of the reprogramming process. Effects of modulation in collective histone state kinetics on transition dynamics and pathway are examined in search for an efficient protocol of reprogramming.

Differentiated mouse cells can be reprogrammed to induced pluripotent stem cells (iPSC) by inducing certain proteins known as Yamanaka factors (Oct4, Sox2, Klf4, and c-Myc) in the cell 1 . Though the precise mechanism of how Yamanaka factors (YF) reprogram remains elusive, clues to determining the mechanism should be obtainable from reprogramming pathways [2][3][4] . On inducing YF, marker genes of the differentiated cells are silenced in the early phase, and pluripotency genes such as Nanog become active only in the late phase, showing that the observed pathway of reprogramming is different from the pathway of differentiation. Therefore, theoretical analysis of how pathways are determined by gene regulation has been a focus of recent interest [5][6][7][8][9] .
Our understanding of gene regulation in differentiation and reprogramming has been advanced particularly by using the concept of epigenetic landscape [10][11][12][13][14][15] . In the landscape picture, stable cell states are represented by basins on the landscape while transition pathways between cell states are determined by topological connectivity among basins. Epigenetic landscape has been calculated in a variety of scenarios 6,9,[16][17][18][19][20][21] , which has shown that landscapes have multiple basins corresponding to differentiated and embryonic stem cell (ESC) or iPSC-like pluripotent states; in order to understand transition pathways, it is necessary to elucidate the distribution of basins and connectivity among them on the epigenetic landscape.
Analysis of the structure of epigenetic landscape so far has been based on the assumption that gene activity is determined by binding/unbinding of transcription factors (TF) as in the case of bacterial gene regulation. However, in differentiation and reprogramming, genes are regulated not only by TF binding/ Scientific RepoRts | 5:16746 | DOI: 10.1038/srep16746 unbinding but also by epigenetic state change including DNA methylation/demethylation, chemical modifications of histones, and the associated change in chromatin structure 22,23 . Therefore, in order to understand epigenetic landscape quantitatively, we need to develop a theoretical framework that explicitly takes epigenetic dynamics into account 16,24 . It should be noted that because of single-molecule nature of DNA, access of regulatory proteins to the gene loci is noisy, leading bursty transcription 25 . Therefore, in order to incorporate noisy genetic and epigenetic influences on the epigenetic landscape, we develop a theoretical framework based on the master equation 16 .
The epigenetic modification of a nucleosome is known to cause recruitment of modifier enzymes effecting the neighboring nucleosomes and causing them to behave similarly 26 . Theoretical [27][28][29][30][31][32] and experimental 29,30 studies have shown that this non-local interaction should bring about the collective change of many nucleosomes in a gene locus to show the discrete switching behavior. In particular, two insightful models 27,32 have been proposed on how memory arises at the epigenetic level by taking long-range interactions of nucleosomes into account 27 or by modeling short-range interactions with a Potts-like model 32 . Based on these observations, we refer to the collective histone state around a gene locus, which constitutes tens or more of similarly modified histones, as the collective histone state (CHS). In Fig. 1 we describe the coarse graining of many-histone states to CHS as a three state switch (s = 1,0,− 1): The state (i) s = 1 corresponds to the loosely packed chromatin state with histones being actively marked and the gene can express itself, (ii) s = − 1 is the tightly packed state with histones being repressively marked and the gene being silenced, and (iii) in the s = 0 state chromatin fluctuates between tightly and loosely packed structures with histones bearing neither repressive nor activating marks as in some loci of ESC 33 . We assume that binding of activator to gene regulatory region of DNA enhances the probability of the s = 1 state and binding of repressor enhances the probability of the s = − 1 state, so that the gene activity, i.e. rate of protein synthesis in the model, is regulated by TF binding/unbinding through the CHS change.
In the present paper, we focus on the role of timescale difference among different processes. Effects of the timescale or the frequency of DNA state change on gene expression have been intensively studied by using the ratio of the rate of DNA state change to the rate of protein-copy-number change as a measure 11,19,[34][35][36][37][38][39][40][41] . This measure is called adiabaticity, and the dynamics is adiabatic when this ratio is large and non-adiabatic when it is small. In the present case, timescale of DNA state change should be determined by epigenetic CHS dynamics; therefore, adiabaticity is measured by ω = / q k hs c , where q is the typical rate of change in CHS and k is the rate of protein degradation with the system being adiabatic when ω > 1 hs c and non-adiabatic when ω < 1 hs c . We show that adiabaticity of epigenetic state dynamics at each gene locus considerably affects the topological structure of the landscape and the slow non-adiabatic epigenetic CHS dynamics brings about the difference in pathway between differentiation and reprogramming in the model. Using this network model, we calculate the kinetic process of cell state transition and compare it with the observed data of reprogramming [42][43][44] . The simulated kinetic behavior predicts a Figure 1. In the present model, histones are considered to be either actively marked (A: acetylated as H3K29ac and methylated as H3K4me3), repressively marked (R: deacetylated and methylated as H3K9me3) or unmarked (U; no modification or bivalently modified). The grey circles on the nucleosomes represent active histones and black diamonds represent repressive histones. Nucleosomes tend to effect neighboring nucleosomes to modify them similarly, which allows us to define collective coarse-grained histone state (CHS). When histones are collectively active, we denote CHS as s = 1, collectively repressive as s = − 1, and collectively undefinable to be A or R as s = 0. The transition between these states depend on the repressor-binding (j) and activator-binding (m) states (subscripts of the rates q and r).
Scientific RepoRts | 5:16746 | DOI: 10.1038/srep16746 characteristic histone-state change along the pathway, which should be examined by single-cell pursuit of the reprogramming process. states. This A -B motif is ubiquitous in regulating differentiation as Cdx2-Oct4 and Gata6-Nanog, for example [45][46][47] . As in these example motifs, we regard A as a marker gene specific to a differentiated cell and B as a pluripotency gene such as Nanog. We then study reprogramming as a transition from the differentiated cell with

A network model
. Each gene in the MRSA network is regulated through CHS by assuming that the rate of protein synthesis from the gene is large only when s = 1 and rates with which s changes depend on the activator or repressor binding status. It should be noted that the MRSA model without the inclusion of CHS dynamics has earlier been used as a prototypical model of differentiation 12,13,18-20 . Along with N A,B , each gene status is described as with CHS (s = − 1,0,1), the repressor-binding state (j = 0: binding, j = 1: unbinding), and the activator-binding state (m = 1: binding, m = 0: unbinding). Protein-production rate, g sjm , is maximal for g 111 = g, and we choose g/k = 1000 to fit a typical protein-copy number of eukaryotic TF 48 . Other g sjm are chosen as g 110 = 0.8 g, g 101 = 0.1 g, and g 100 = 0. For the CHS s = 0 and − 1, we use g 011 = 0.2 g with all other g 0jm = 0 and g −1,jm = 0. Protein degradation rate is k ≈ 0.1 h −1 49 , and length of a cell cycle is about 2k −1 42 though we do not consider cell cycle explicitly. For simplicity, we adopt same values of k and g sjm for A and B. We consider that the rate constant of protein binding h is affected by other factors not involved in the present network model, which competitively or cooperatively interact with A and B when they bind to DNA; therefore, we assume the binding rate constants h a (A), h a (B) for activators and h r (A), h r (B) for repressors take different values. We assume proteins bind to DNA as a dimer 9 , so that the rates of binding are h a (α)N α (N α − 1) and h r (α)N β (N β − 1) with α,β = A, B, or B, A.
iPSC are unstable when they are cultured in medium without Lif, and stabilized when they are cultured in ES medium which contains Lif. Therefore, the iPSC state with . We set f/h = 50000 to make the ratio /( ) < α f hN 1 2 . Following the observed data of single-molecule measurement 50 , we assume that TF binding/unbinding is fast enough as f/k = 10; With such fast TF binding/unbinding, the other slower process, the CHS transition, should be a key determinant of adiabaticity of the DNA state change in the present model.  56,57 to change the CHS. The 3-state CHS switch controls the protein production rate g as shown.
Rates of stochastic CHS transitions are q jm for s = 0→ 1, r jm for s = 0→ 1, ′ q jm for s = − 1→ 0, and ′ r jm jm are chosen as real multiples of a tuning rate q. Therefore, the adiabaticity measure is ω = / q k hs c . We assume that the CHS tends to be turned active when the activator binds, and turned repressive when the repressor binds to DNA. We therefore, have  q q 11 00 and ′ ′  r r 11 00 . See Methods for further details.

Results
We first study how topology of epigenetic landscape is influenced by the CHS dynamics, and then discuss pathways and kinetics on the landscape.
Topology of epigenetic landscape. We calculated steady state distribution P s (N A , N B ) by simulating the stochastic equations described in Methods using the Gillespie algorithm 51 to derive the epigenetic landscape: ( , . 100 trajectories were used for sampling, each over 10 8 Gillespie-steps long with random initial conditions. The role of CHS switching is studied in the spectrum of adiabatic to non-adiabatic timescales. Figure 3 shows topological changes manifesting on the epigenetic landscape in both symmetric and asymmetric models. In symmetric landscapes, two distinct states appear as two basins. In the strongly adiabatic case with ω = 10 hs c ( Fig. 3a), features of CHS dynamics are averaged out, and two basins are separated by a large epigenetic barrier as in the MRSA network without CHS dynamics. Here, a path connecting two basins through the barrier is referred to as diagonal pathway. With the intermediate adiabaticity of ω = 1 hs c (Fig. 3b), the barrier is washed away due to the resonance between CHS dynamics and protein copy-number dynamics. This flat landscape should result in the widely fluctuating cell states, which disagree with the observed narrow distribution of cell states along the reprogramming pathway. In the non-adiabatic case with ω = .
0 05 hs c (Fig. 3c), both differentiated and iPSC states are stably formed, and we find the emergence of a stable intermediate state with . This novel intermediate state is connected to the differentiated as well as the iPSC states through low barrier pathways (valleys). Here, we refer to these valleys arising out of non-adiabatic CHS dynamics as epigenetic valleys or epigenetic pathway.
In asymmetric landscapes, the differentiated state has a deeper basin due to the enhanced stability. In adiabatic case with ω = 10 hs c (Fig. 3d), the iPSC basin vanishes. In the intermediate case with ω = 1 hs c (Fig. 3e), the population widely spreads to give a flat landscape as in the symmetric model, which does not support a stable cell state. In the asymmetric non-adiabatic case with ω = .
0 05 hs c (Fig. 3f), we find differentiated, iPSC, and intermediate basins, which are connected by diagonal and epigenetic pathways, and further, we find another low-lying pathway in between diagonal and epigenetic pathways, which is characteristic to asymmetric non-adiabatic landscape.
We have shown that the topology of the epigenetic landscape is decisively dependent on the adibaticity of the CHS dynamics. We should note that without the CHS dynamics, the landscape does not have epigenetic valleys or the intermediate state with , but only shows a diagonal pathway as in Fig. 3a. In the next subsection we examine the epigenetic pathway induced by the non-adiabatic CHS dynamics and its role in the reprogramming mechanism.
Pathways of transitions between cell states. We investigate the role of pathways induced by the CHS dynamics by numerically solving the master equation of the model. See Methods for the explicit form of the equation and the calculation details.
Since ES medium is used in reprogramming, we use the symmetric network model to study it. The precise mechanism of YF action is not known; therefore, we compare two possible mechanisms; (I) YF work as histone-mark erasers by changing the CHS as . Reprogrammed trajectories start from the equilibrium differentiated state. In the adiabatic case with ω = 10 hs c , we show two trajectories X 1 (t) and X 2 (t) with efficiency C 0 /k = 100 and 10, respectively (Fig. 4a). The adiabatic nature of the CHS dynamics makes the trajectories independent of the protocol γ. YF is inefficient in the case of X 2 (t) causing a reversal, in contrast X 1 (t) is able to reach the iPSC state, this signifies that ω /  C k hs 0 c for YF to be sufficiently effective. In the non-adiabatic case with ω = .
0 05 hs c (Fig. 4b), we study five trajectories X 3 (t) with γ = 0 and C 0 /k = 0.1 and X 4 (t), X 5 (t), X 6 (t), X 7 (t) with γ = 1 having efficiency C 0 /k = 5.0, 0.05, 0.1, 0.5, respectively. Starting from the differentiated state, owing to a largest efficiency, X 4 (t) gets closest to the intermediate state and surpasses the epigenetic barrier to reach the iPSC. X 5 (t) and X 6 (t) have small efficiencies, so that they reverse back to the differentiated state. X 7 (t) with a medium efficiency crosses the epigenetic barrier but does not get as close to the intermediate state as X 4 (t). On keeping γ = 1 fixed and increasing C 0 , the system goes closer to the intermediate state. By decreasing γ and keeping C 0 fixed, on the other hand, the trajectories depart from epigenetic valleys. as is apparent with the case X 3 (t). Trajectory X 3 (t) bypasses the epigenetic barrier suggesting that rapid reprogramming is realized along this pathway. Thus, the pathway of reprogramming sensitively depends on the way how YF work when the CHS dynamics is non-adiabatic. We can expect this rich behavior when kinetics of reprogramming is experimentally analyzed.
In Fig. 5, we show the time evolution of the probability distribution, ( , , ) , for the X 4 (t) (1(i)-(vi)) and X 6 (t) (2(i)-(vi)) cases with ω = .   , X(t), are drawn on the epigenetic landscape, U(N A , N B ). which exhibits reversal, the major part of ( , , ) P N N t A B remains near the differentiated state, but the minor part proceeds along the epigenetic valley and reaches the intermediate state. Importantly, this approach to the intermediate state in the non-adiabatic reprogramming process is consistent with the experimentally observed late activation of the pluripotency genes after the lineage specific genes being repressed [2][3][4] . Since only the non-adiabatic case explains the pathway through such intermediate state, and epigenetic valleys are absent with the fast adiabatic CHS dynamics, the model strongly suggests the importance of the slow CHS dynamics. This result is also consistent with the slow CHS dynamics observed by Hathaway et al. 30 : Histone modifications around the Oct4 locus in mouse ESCs are spatially week, which suggests ω ≈ . 0 05 hs c , though the turnover rate of single nucleosome is as fast as 10k 52 . In X 6 (t) at later time instances, the tail of the population passes through the intermediate and further approaches the iPSC state. In this way, even when the average of the population represented by X(t) reverses, the tail of the population reaches the iPSC state. This behavior is consistent with the experimentally observed low efficiency of reprogramming, in which only the small portion of the cell population reaches the iPSC state.
For differentiation, we use the asymmetric model because Lif is withdrawn form cultivating medium. The simulated temporal evolution of ( , , ) is shown in Fig. 6 with the non-adiabatic CHS dynamics (ω = . 0 05 hs c ). Starting from the pluripotent state by keeping C 0 = 0, the population shifts along the diagonal pathway due to the bias in the asymmetric landscape, which passes through the saddle to reach the differentiated state. Though a minor population of simulated differentiating cells proceed along the epigenetic valley, major part of cells go through the diagonal pathway. We can expect that cells behave in a collective way through cell-cell communication, which should further enhance the major pathway. Thus, the difference between differentiation and reprogramming pathways is evident with the non-adiabatic CHS dynamics. In reprogramming, two genes were set to be symmetric in the binding affinity of transcription factors to gene loci in the model, but in differentiation two genes are asymmetric, which brings about the difference in pathway between reprogramming and differentiation, but this difference is realized only in the slow non-adiabatic histone switching regime.
In Fig. 7, temporal evolution of the probability of each CHS, P(s α , t) for α = A and B, is plotted for the reprogramming process with the non-adiabatic CHS dynamics for X 4 (t). We can see that starting from the state in which Kinetics of reprogramming. We investigate the role of non-adiabatic dynamics in reprogramming kinetics by simulating latencies observed in the ensemble-level experiments of Hanna et. al. 42 . In the report of Hanna et al. 42 , N col founder cells infected with YF were placed in N col wells on a plate at t = 0 to multiply and form clones. Population of these cells in each well exponentially increased from 1 to 10 6 and saturated in t > 10 days. iPSC were detected through the Nanog expression measurement. The

P t PN N t N N A B A B
is obtained by solving the master equation and imposing an absorbing boundary condition at the iPSC state (see Methods). Assuming the cells in the effective population can be regarded as independent, we have: Figure 8a shows Q 5 (t), Q 6 (t) and Q 7 (t), which are the Q(t)s corresponding to the non-adiabatic trajectories X 5 (t), X 6 (t) and X 7 (t) of Fig. 4b. In Fig. 4b, X 5 (t) and X 6 (t) did not reach the intermediate state and revert in the epigenetic valley, which implies that the mean of ( , , ) P N N t A B does not surpass the intermediate state with small C 0 . With this parameter set, however, the tail of ( , , ) passes through the intermediate basin (Fig 5 2(i-vi)) and reaches the iPSC state, and thus brings about the slow rise of Q 5 (t) to Q 5 (t) ≈ 1 as shown in Fig. 8a. For both Q 5 (t) and Q 6 (t), Q(t) ≈ 0 initially and starts to rise at t 0 and reaches Q(t) ≈ 1 at t 1 showing that colonies have heterogeneous latencies. For Q 5 (t), we have t 0 ≈ 28k −1 and t 1 ≈ 95k −1 , which agrees with the experimental values 42 : t 0 ≈ 30k −1 and t 1 ≈ 100-200k −1 .
We find that there are two ways to accelerate reprogramming. One is to increase the efficiency C 0 : With the larger C 0 as in Q 7 (t) (Fig. 8a), we have the sharper increase of Q(t). The other is to decrease γ. With γ = 0, Q(t) increases rapidly with t 0 ≈ 8k −1 and t 1 ≈ 10k −1 for N eff = 10 5 , which is similar to the observed data with t 0 ≈ 8k −1 and t 1 ≈ 12k −1 obtained for cells in which Mbd3, a factor which binds to the methylated DNA, is silenced 43 . Thus, when YF work as histone-mark erasers with mild efficiency, reprogramming has heterogeneous latency distribution, but when YF work with high efficiency or they work as activators of pluripotency genes, reprogramming is accelerated with lesser degree of heterogeneity or is more "deterministic" 42,43 in latencies.
The experimental data of Hanna et al. on kinetics of reprogramming 42 was simulated also by Morris et al. 53 by regarding reprogramming as a diffusion process on a phenomenological model landscape having two basins at the differentiated and iPSC states. However, role of N eff was neglected in their argument 53 . In Fig. 8a, we find the slope of Q(t), which is the degree of "determinism" of reprogramming experiments 43,54 to be very sensitive to N eff . Therefore, attention is needed to compare the data with different N eff .
The γ dependence of variation in heterogeneity of latencies is apparent in R(t) shown in Fig. 8c. Increase in R(t) is much sharper for low γ. The γ dependence of the localization properties in the N A -N B space can be studied via the participation ratio: , which is large when the distribution is localized and small when delocalized. Fig. 8d shows that the distribution gets more localized with increasing γ during the reprogramming. For γ = 1 population gets accumulated around the intermediate state. Localization pattern is found to be more complex in the γ = 0 case. These features should be verifiable by single-cell level tracking.

Discussion
We have introduced a minimal model of reprogramming by integrating epigenetic modification dynamics with the gene expression mechanism, which provides a consistent view of pathways and kinetics. We showed how pathways are created on the epigenetic landscape aided by the slow epigenetic dynamics of collective histone-state modification. With the slow non-adiabatic epigenetic dynamics, landscape has epigenetic valleys which connect the differentiated and pluripotent cell states through the intermediate state. This pathway is consistent with the observed late activation of the pluripotency genes in the reprogramming process. The time course of the histone-state change and the extent of localization of distribution of cell states were simulated, which provide clues to examine the mechanism of the observed reprogramming process.
Based on the simulated temporal evolution of cell distribution, kinetics of reprogramming was calculated. The calculated kinetics is consistent with the experimental observation when YF is assumed to work as histone-mark erasers. With the non-adiabatic epigenetic dynamics, kinetics is sensitive to the precise mechanism of how YF work: When YF work more actively to promote expression of pluripotency genes, then the model suggests the accelerated kinetics with which reprogramming becomes "deterministic".
We show that when we explicitly take epigenetic dynamics into account, we can explain features of the trajectory and barriers helping experimentalists with microscopic information which is otherwise difficult to obtain. The MRSA network model coupled to CHS used here is a minimal model developed for elucidating the roles of epigenetic dynamics and it is important to apply concepts and methods developed here to more realistic networks which represent antagonistic interactions between pluripotency genes and differentiation marker genes [6][7][8][9]16 .

Methods
We first summarize in the subsection Stochastic equations how the gene-state transitions are modelled in the CHS coupled MRSA network we discussed. The details of the model parameterization are also discussed in this subsection. The master equation governing the stochastic equations is explained in the subsection Master equation. We use the proteomic field approximation 11,34,55 to reduce the dimensionality of the master equation so that the solution is tractable computationally. Details used for kinetic calculation are explained in the subsection Effective number of cells.

( ) P N i j m N i j m P N i j m P N i j m ; 9 A A A A B B B B A A A A B B B B
For the binding terms in H used in the master equation, we use an approximate form: 1. This is a reasonable approximation for copy number ≳ 10 2 . We solved the resulting proteomic field equations using the fourth order Runge-Kutta method without using the adiabatic approximation.
The trajectories ( ( ) , ( ) ) N t N t A B are calculated without any absorption condition. For calculating R(t) and related quantities, an absorption condition was imposed as follows: P(N A , N B , t) = 0 for N B > 910, and P (N A , N B , t)  Effective number of cells (N eff ). Histone states are inherited from mother to daughter cells, there is bound to be correlation among multiple cells in a colony. The effective number of cells, N eff , therefore, should be smaller than the actual number of cells in a colony. We here used N eff = 10 4 .