Modeling the Dynamics of Bivalent Histone Modifications

Epigenetic modifications to histones may promote either activation or repression of the transcription of nearby genes. Recent experimental studies show that the promoters of many lineage-control genes in stem cells have “bivalent domains” in which the nucleosomes contain both active (H3K4me3) and repressive (H3K27me3) marks. It is generally agreed that bivalent domains play an important role in stem cell differentiation, but the underlying mechanisms remain unclear. Here we formulate a mathematical model to investigate the dynamic properties of histone modification patterns. We then illustrate that our modeling framework can be used to capture key features of experimentally observed combinatorial chromatin states.


Introduction
Histones can undergo various types of covalent modifications, such as methylation and acetylation, which serve as an additional layer of transcriptional control by mediating the chromatin accessibility and by recruiting regulatory proteins [1,2]. Experimental studies using chromatin immunoprecipitation followed by massively parallel sequencing (ChIP-seq) have suggested that different cell types can be characterized by different histone modification patterns [3].
The molecular mechanisms underlying chromatin state establishment, maintenance, and heritability remain incompletely understood. A number of mechanisms are implicated [4], including (1) sequence-specific recruitment through interactions between chromatin regulators and DNA binding factors; (2) recruitment of chromatin regulators to existing histone marks; (3) histone marks deposited by transcriptional machineries; (4) RNA mediated recruitment; and (5) stochasticity associated with DNA replication. However, any single mechanism alone is insufficient for chromatin state establishment [4,5].
One of the best characterized chromatin states is a bivalent domain, a segment of the nucleosome array, in which H3K4me3 (an active mark) and H3K27me3 (a repressive mark) coexist on most individual nucleosomes within the domain [6]. Bivalent domains are thought to be an important feature of stem cells. For example, bivalent domains have been discovered in the promoters of most lineage-control genes in embryonic stem cells, and most of these domains become monovalent upon cell differentiation [3,[6][7][8][9]. Also, a recent study observed that gene activation in the differentiation process occurs in conjunction with the decay of repressive marks in bivalent domains [10]. In particular, one prominent proposal [6] for the function of bivalent domains is that the H3K27me3 marks act to repress the lineage-control gene in stem cells, while the H3K4me3 marks poise these genes for activation upon cell differentiation. Thus this proposal suggests that activation of these genes in differentiated cells is determined by the existence of bivalent domains in stem cells. These findings indicate the importance of bivalent domains and motivate further study in order to illuminate the underlying principles and mechanisms involved in their formation and evolution.
It has been proposed that the formation of chromatin domains is consistent with a model that includes not only the chemical interactions between histone marks, but also nucleation sites where domains are more likely to form [11]. The dynamics of histone modifications have been studied both theoretically and experimentally for some time [11][12][13][14][15]. In general, histone methylation marks are catalyzed by a variety of methyltransferase enzymes which may act singly or cooperatively. For example, H3K27me3 marks are catalyzed by Ezh2, a core member of the Polycomb group proteins. In addition to the normal stochastic conversion which would be expected from each of these individual enzymes, there is also a feedback process between the histone marks and the enzymes [16]. Existing H3K27me3 marks may attract Polycomb group complexes, which enhance nearby methylation [17,18]. A similar recruitment mechanism has also been suggested for H3K4me3 via Trithorax protein complexes (TrxG) [19]. In addition, there exists experimental evidence supporting a negative feedback mechanism between H3K4me3 and H3K27me3 marks via the action of histone demethylases [20][21][22][23][24].
Certain specific DNA sequences may serve as the docking sites of modification enzymes and may therefore be associated with enhanced local attraction of histone marks [4,25]. We refer to these as nucleation sites. For example, CpG islands are strongly enriched in bivalent domains in human and mouse embryonic stem cells [19], and appear to be required for Polycomb binding in certain cases [26].
Recently, in silico methods have provided important additional insights for chromatin state inheritance. Major contributions have been made by Dodd el al. [12] and Sedighi and Sengupta [27]. These paper considered 1-dimensional lattice models in which nucleosomes are allowed to have active or repressive modifications that evolve stochastically and by recruitment. They found that a bistable state with either mostly active nucleosomes or mostly repressive nucleosomes can appear and be heritable, consistent with experimental observations. Subsequently, Hodges and Crabtree [11] found that adding a nucleation site into a model of the above type produces a bounded chromatin domain. Also, in a more recent paper, Binder et al. [28] proposed a model describing binding of catalytic enzymes to DNA and their interaction with histone marks with one aim being explaining length distributions of modified chromatin regions. These past studies are limited to a single type of histone mark on a nucleosome, whereas it is well-known that gene regulation is governed by combinatorial patterns of multiple histone marks [2,29]. In this paper, we extend previous studies by presenting an approach to model the dynamics of combinatorial chromatin states. This is achieved by allowing each individual nucleosome to carry both active and repressive marks simultaneously.
In the next section we describe our model. Then, in the Results section, we apply this model to investigate the dynamics of histone modification patterns with the focus on bivalent domains. Discussion and Conclusions are given at the end of the paper.

General Framework of our Model
We consider a 1D lattice of N nucleosomes, where there is a nucleosome at each lattice site i~1,2,:::,N. An actual nucleosome consists of 8 histone protein molecules, that can be regarded as two identical groups of four each. In what follows we only consider the state of one of these four histone group members, namely the H3 histone, which is specifically related to bivalency. Thus, in our model, we represent the state of a nucleosome as being determined by the states of its two H3 histone copies. There are two modification sites in each H3 histone, one which may have an active mark (such as H3K4me3) and the other which may have a repressive mark (such as H3K27me3). Thus, there are 16 possible states of a nucleosome, and each of which is determined by 4 histone modification sites (see Figure S1). As shown in Text S1, this, together with the restriction obtained from experiment [30] that active and repressive marks do not occur simultaneously on the same H3 histone, leads to the six physically distinct nucleosome states depicted in Fig. 1. In Fig. 1 the circle represents a nucleosome and the vertical ellipses represent H3 histones. The lower case letters within each ellipse represent the states of the two modification sites of the H3 histone (u = unmodified, a = modified by an active mark, r = modified by a repressive mark). For convenience, we assign the symbols UU, AA, RR, AU, UR, and AR to the six possible states. From now on, when we say 'histone' it is to be understood that we mean an H3 histone. We note that the state AR will play a prominent role in subsequent considerations in the Results section, and we will call a nucleosome in this state a 'bivalent nucleosome'.
We then allow each nucleosome state to evolve according to a discrete time (t) model, in which from time t to time tz1, a nucleosome state changes from state s to state s 0 with probability p ss 0 . Since the time step t?tz1 is regarded as small, we assume that, at most, only one modification site may change on each time

Reduced Model
The general framework above can lead to a relatively complex class of models and has many parameters. Thus, for the simulations that we report in this paper, we have adopted the somewhat modest goal of illustrating different types of dynamics that can arise when different nucleosome states interact and compete. With this goal in mind, we now seek an illustrative, but still somewhat plausible, reduction of our general 6-state model. Our reduction is based on the assumption, motivated in Text S1, that the occurrence of nucleosome states having either active marks on both histones (AA in Fig. 1) or repressive marks on both histones (RR in Fig. 1) are unlikely. Thus we consider the idealized case where AA and RR states do not occur. Hence each nucleosome of the reduced model is in only one of 4 nucleosome states, namely AU, UR, AR and UU (see Fig. 3A). Referring to Fig. 1, we see that the four states have the following meanings.
AU: One histone has an active mark and the nucleosome's other three sites are unmodified.
UR:One histone has a repressive mark and the nucleosome's other three sites are unmodified.
AR:One histone has an active modification, while its other site is unmodified. The other histone has a repressive modification, while its other site is unmodified.
UU:All four sites of the nucleosome are unmodified.

Model Dynamics
During a cell cycle, we consider the time t states of modeled nucleosomes on our one dimensional lattice and update these states to new states at time tz1 through two probabilistic processes that we call ''recruitment conversion'' and ''exchange conversion''. At the conclusion of a cell cycle, ''replication'' occurs, following which a new cycle begins. N Recruitment. This refers to the recruitment of histone marks to a nucleosome through interaction with neighboring nucleosomes. Recruitment at a site i depends on the states of the nucleosomes in an interval of length 2l centered at i, and we refer to l as the range of recruitment. We define f i X as the fraction of nucleosomes in this interval which carry a type-X histone mark, where the subscript f i X is X~A,R. If lƒiƒN{l, then the recruitment range will span 2lz1 nucleosomes on our lattice. However, if i is too close to the beginning or the end of the lattice (i.e., 1ƒivl or N{lviƒN, respectively), then the recruitment range will include 'phantom' sites j not on the lattice (jv1 and jwN, respectively), and for the purpose of determining f i X , we consider such phantom sites j to be in the UU state. The probability of recruitment conversion from U to X at site i is taken to be given by f i X r UX , where r UX is a constant describing the strength of the recruitment interaction. On the other hand, the probability of recruitment conversion from X to U (i.e., mark removal) depends on the concentration of histone marks which are opposite (rather than similar) to X (where we regard A and R as opposites). In this case, the conversion probability is taken to be given by Fig. 3B). Note that in our model we allow r XU to differ from r UX because different enzymes are recruited for the addition and removal of histone marks. N Exchange. Unlike the recruitment process, the exchange process refers to histone modifications which occur spontaneously, independent of the states of nearby nucleosomes. The probabilities for exchange conversion are denoted by p UA , p UR , p AU , and p RU (see Fig. 3C). In particular, we think of p AU and p RU as corresponding to the histone turnover process, and p UA and p UR as corresponding to processes involving nucleation sites (See Table 1).
N DNA replication. When DNA replication occurs, we imagine that in the real situation the parental nucleosomes are randomly assigned to one of the two daughter strands at the same site as that which they occupied on the parental strand, while the corresponding site on the other strand is assigned an unmodified nucleosome (i.e., a nucleosome in the UU state). This scenario is supported by an experimental observation [31]. In our model, we do not follow both daughter strands. Rather we follow just one. Thus, with probability 1/2, our model replication process randomly replaces each nucleosome with an unmodified (UU) nucleosome. This model DNA replication occurs periodically with a period equal to the 'cell cycle time' t. This is similar to how replication is modeled in [12].
In accord with the above recruitment and exchange processes, during a cell cycle, our model gives appropriate equations for the probabilities P i XY (tz1) that nucleosome i is in state XY~UU,AU,UR,AR at time tz1, given the state of the lattice at time t. After the probabilities P i XY (tz1) are determined the state (UU, AU, UR or AR) of each nucleosome i is randomly chosen according to the probabilities P i XY (tz1), thus determining the state at time tz1.
if nucleosome i is not in state XY , our model equations for the probabilities are Consistent with our assumption that at most one site on a nucleosome can change state in one time step, our choice of parameters satisfies r XY ,p XY %1. Note that f i A (t) and f i R (t) depend on the lattice state in a neighborhood of site i within the range of recruitment specified in the first bullet above.
In section 4.3, where we treat localization of AR states, we allow the exchange transitions probabilities p i XY to vary from site to site, but everywhere else we consider p i XY to be the same at each site, p i XY~p XY :

Simulation Parameters
To assign roughly reasonable values to the parameters r XY and p XY , we first consider that our model time step, t?tz1, corresponds to a real time step Dt~2 min. We have numerically verified that our simulation results are independent of our choice of Dt so long as Dt is sufficiently small. To estimate a rough range for the parameters r XY and p XY , we set p XY , r XY &(Dt=T), where T is the characteristic time scale of the relevant process (see Table 1), and, as required, the Dt that we have chosen is such that Dt=T is small compared to one for all such processes. We fix as many parameters (Table 1) as possible using experimental information (see Table 2). Because the authors are not aware of any experimental measurements of the characteristic time for recruitment demethylation and methylation via exchange, we will consider these probabilities as free parameters in our numerical simulations below. Previous work [32] suggests that the loss of active marks is faster than the loss of repressive marks. In particular, it has been shown that nucleosome turnover is faster in regions bound by trithorax-group proteins. Therefore, we selected the model parameters so that all rates associated with active mark are faster than those associated with the repressive mark.
Specifically, we assume that r UR =r UA~pUR =p UA~rRU =r AUp RU =p AU~0 :5 in the simulation (when nonzero). Regarding the cell cycle, for embryonic stem cells the cell cycle length is about 12 hours, which, with our Dt~2 min, corresponds to 360 time steps of our discrete time model per cell cycle. Finally, motivated by Ref. [31], we take l~2, corresponding to a fairly short range of recruitment.

Results
We now illustrate the utility of our model by employing it to investigate dynamic changes of histone modification patterns. As described in the Introduction, both nucleation sites and recruitment of methylation may be involved in the establishment of bivalent domains. As noted above, we suggest that certain nucleosomes act as nucleation sites during the early stages of development. These nucleation sites may be instrumental in the formation of bivalent domains. We incorporate nucleation sites into our model by assigning them a higher value of p UA and p UR than other sites, and we model the absence of nucleation sites by lowering its value of p UA and p UR .
In Sections 4.1 and 4.2, we discuss the formation and decay of AR states with different initial conditions in the absence of nucleation sites. In Section 4.3, we study the effect of nucleation sites on dynamics of the formation of AR states. Finally, in Section 4.4, we consider how varying the cell-cycle length affects AR states. Taken together, these analyses demonstrate the utility of our model for systematic investigation of the dynamic properties of bivalent domains.

Formation of AR States
The formation of bivalent domains has been experimentally observed in studies of the early stages of embryogenesis [33] and in studies of cell reprogramming [34]. In particular, studies of cell reprogramming observe this formation process to be gradual [35].
In this section we use our model to simulate the formation of regions that are dense with AR states, and we identify such regions with bivalent domains. In the simulations, we take p UA~pUR~0 for all nucleosomes and fix r UA = 0.046 (corresponding to an H3K4me3 methylation timescale of 30 mins) and p AU~0 :005. Also, r AU and r RU are considered to be very small (for simplicity, we set r AU = r RU = 0), so that the AR states can be established and  persist for a long time. For the initial state of the lattice in the simulations, we consider a situation where there are a relatively small number of nucleosomes in AR states, with all other nucleosomes initially in the UU state. In particular, we choose the initial number of AR nucleosomes to be five (out of the 80 nucleosomes on the lattice), and we study how AR states spread to other nucleosomes on the lattice. To investigate the effect of the initial spatial distribution of AR nucleosomes, we consider two extreme cases: the localized case in which all five initial AR state nucleosomes are located at five consecutive nucleosome sites in the center of the lattice, and the delocalized case in which the five initial AR state nucleosomes are located at equally spaced sites spanning the entire lattice (at sites 1, 20, 40, 60, 80). Fig. 4 shows results for the space-time evolution of the distribution of nucleosomes for both localized (left column of figure panels ) and delocalized (right column of figure panels) initial states. For the localized case, there appears to be two fronts, a fast UU?AU front (corresponding to the blue to yellow transition in the left panel of Fig. 4C), followed by a AU?AR front (yellow to red transition in the left panel of Fig. 4C) that propagates at a slower speed than the UU?AU front. The slow AU?AR front is clearly seen in the left panels of Figs. 4A and B, while the UU?AU front is evident in the left panel of Fig. 4B. These two fronts propagate symmetrically in space in the average space-time plots (Fig. 4A and 4B) but, due to fluctuations, more asymmetrically in space in the single run plot (see left panel of Fig. 4C). Examining a range of parameters, we find that the fastest front corresponds to either a UU?AU transition (as in Fig. 4B) or a UU?RU transition (not shown). For the delocalized case, we also observe that the spreading of the active marks is faster than that of the repressive marks. This can be easily seen from the typical single run plot in the right panel of Fig. 4C.
Finally, we also studied the effects of varying the number of AR nucleosomes in the initial condition on the above simulations. Using the same parameters values as above, we plot (Fig. 5) the final average fraction of AR nucleosome at the end of the final simulated cell cycle (10 cell cycles) as function of the initial number m of AR nucleosomes which are taken to occupy the m nucleosome sites in the center of the lattice. As shown in Fig. 5, the average fraction of final AR nucleosomes initially increases with increasing m. We observe that past m §4 (i.e., m=80 §0:05) the value is essentially constant up to m~80 with AR nucleosomes spanning the whole lattice. For a given m, each simulation can be categorized into two groups, (1) the final spatial average level of AR nucleosomes is approximately equal to the corresponding large m limiting value, or (2) all AR nucleosomes vanish. Thus at low m, the value plotted on the vertical axis of Fig. 5 can be thought of as the limiting larger-m value (basically the value at m~4) multiplied by the fraction of runs in category (1). In the early stage of a simulation, the spreading of histone marks compete with the loss of histone marks via histone turnover. If either type of mark is lost totally, it cannot recover (i.e., the run is in category 2). On the other hand, we find that histone marks do not die out if there are enough of them on the lattice (the run is then in category 1). As a result, the average fraction of AR nucleosomes is larger with larger m, and with smaller p AU and p RU (compare the red and blue plots in Fig. 5). The above simulations suggest that in order for AR states to form when p UA and p UR are small, a sufficient number of initial AR nucleosomes is required. Taken together, these results have shown that the formation of bivalent domains undergoes two distinct phases: expansion and stabilization. In the expansion phase, the border of bivalent domains expands to neighboring nucleosomes. The expansion process is relatively fast (w10 nucleosomes per cell-cycle in our simulation) but quite noisy. As a result, only a sparse subset of nucleosomes are marked with the AR state. During the stabilization phase, the nucleosome state configuration is further refined and eventually reaches an equilibrium. Even then, the state of individual nucleosomes is still highly dynamic and equilibrium is only reached in the statistical sense.

Decay of AR States
In this section we use our model to simulate the decay of AR states. All parameters are the same as in section 4.1 except that r AU and r RU are taken to be non-zero. This is motivated by experimental findings that recruitment of demethylases is important for the decay of bivalent domains [22,23], and occurs during cell differentiation. Also, we consider an initial condition in which all nucleosomes are in AR states. Results are shown in Figs. 6, 7, and 8 for different values of r AU and r RU keeping their ratio fixed at r AU =r RU~2 . Fig. 6 shows results for the space-time evolution of the distribution of all four nucleosome states AR, UR, AU, and UU for three values of r AU :2r RU . In Fig. 6A, for the case r AU~0 :004, the initial level of AR nucleosomes rapidly (in about one cell-cycle) drops to a lower level of AR nucleosomes, but there still remains a substantial presence of AR nucleosomes which persists to the end of the run. In contrast, for both r AU = 0.016 and r AU = 0.034, where there is again similar very rapid decreases of the level of AR nucleosomes, now the final level is essentially zero. In addition, it is seen that the level of AR nucleosomes takes longer to fully decay for r AU~0 :016 than for r AU~0 :034. The latter case is consistent with the experimental observations [10,35] that an essentially complete loss of bivalent domain can occur very rapidly. To further explore how the decay of AR states depends on the recruitment demethylation rates, we plot the fraction of simulation runs that have at least one AR nucleosome on the lattice as a function of time in Fig. 7, and the final average fraction of AR nucleosomes (averaged over 1000 runs) as a function of r AU in Fig.  8. Comparing Fig. 6A to Fig. 7, we observe that the fraction of runs with at least one AR nucleosome plotted in Fig. 7 shows a slower decay compared to the decay of AR levels in Fig. 6A. This suggests that lineage-control genes in bivalent domains may become active without the full destruction of repressive marks. In Fig. 8, as might be anticipated, we observe that, in general, smaller histone turnover (p AU ) and smaller recruitment demethylation rate give a higher final average fraction of AR nucleosomes. Also, the value of r AU at which the average fraction of AR nucleosomes drops to zero is lower for larger p AU . Our results suggest that a large recruitment demethylation rate in a cell is important for cell differentiation. This is consistent with experimental findings [22,23].
In a real situation, a change from low to high values of the recruitment demethylation rates during cell differentiation will take place by processes not included in our model, and these processes may take some time. Thus our simulation use of constant non-zero initial r AU and r RU results in a determination of the characteristic decay time associated only with processes that are included in our model, and the true decay rate of AR state nucleosomes may be longer than this time due to the finite time for r AU and r RU to change. Overall, we observe that the decay determined from our model of AR state nucleosomes in response to high initial value of recruitment demethylation rate is relatively fast, as compared to the time that it takes to establish AR states spanning the lattice in Section 4.1. We conclude from this that processes included in our model do not prevent rapid decay of AR state nucleosomes, and that rapid decay, as seen in experiments [10], can occur in response to rapid increase of r AU and r RU .
In addition, it is interesting to emphasize the probabilistic nature of these results. For example, Fig. 6D shows results of typical single realizations. This figure also shows that the final state for r AU~0 :016 is different from that for r AU = 0.034. For the case r AU = 0.016, we observe that AU nucleosomes are dominant in the lattice at the end of the simulation (see also the second panels of Figs. 6B and 6C). However, for the case of r AU = 0.034 at long time, green regions of UR nucleosomes form at the upper edge (see third panel of Fig. 6D), while the AU nucleosomes are at the lower edges. This is because the AU and UR states can both be stable for this combination of parameters (see third panels of Figs. 6B and 6C). Our results suggest that the strength of the recruitment demethylation (i.e., the values of r AU and r RU ) is not only important for the decay of bivalent domains, but also strongly influences the possible final state following decay.

The localization of AR States
The next issue that we discuss is the effect of nucleation sites (i.e., in our model, p UA~pUR w0 at these sites). The existence of such sites is suggested by the finding [4,25] that DNA specific sequences can recruit protein binding factors like transcription factor which in turn recruit histone marks to the DNA. In section 4.1, we took p UA~pUR~0 , and we found that AR nucleosomes either span the whole lattice or disappear. Although similar broad bivalent domains are observed, narrow bivalent domains are also detected in some experiments [3,24]. A recent model [11] has previously been used to simulate the dynamics of localized histone modification domains, but that model allowed only a single type of histone modification, and therefore it cannot address the dynamics of bivalent domains. Using our model, we will be able to analyze interactions among the placements of active and repressive histone marks, histone turnover rate, and crosstalk between active and repressive histone marks. We consider p UA and p UR w0 for the central nucleosome (corresponding to the case that the central nucleosome is a nucleation site). For the initial condition, we consider that there are five AR nucleosomes located at the five consecutive nucleosome sites in the center of the lattice, with all other nucleosomes initially in UU state. Using our previous parameter ratios (i.e., r UR =r UA~pUR =p UA~rRU =r AU~pRU =p AU~0 :5), we explore the parameter space regions for which our model reproduces narrow and broad distributions of AR nucleosomes.
We consider cases of both relatively small and relatively large recruitment demethylation rates (r AU and r RU ). The former and latter choices are meant to simulate cell environments far before, and during, cell differentiation, respectively. We run the simulations for four cell-cycles such that the averaged nucleosome state configuration reaches an equilibrium. In particular, a steady spatial distribution of AR nucleosomes seems to be reached within the first cell-cycle, and change very little thereafter. Therefore, the time for establishment of a highly localized AR distribution (v 1 cell-cycle) is much shorter compared to that of establishing a very broad and uniform AR distribution (about 5 cell-cycles and see Figure 4A). For the case of small recruitment demethylation rate,  distributions of AR nucleosomes. The widths of these bounded distributions reflect the balance between the continuous placement of histone marks on the nucleation site, the spreading of histone marks by the recruitment process, and the destruction of histone marks via exchange [11]. From the simulations, we find that the width of the distributions of AR nucleosomes depends more on p AU and p RU , which they are inversely related to the width of the AR distribution. On the other hand, the amplitude of the distributions depends more on p UA and p UR (i.e., the continuous placements of histone marks on the center nucleosome) (Figs. 9A-B).
Next, we did simulations using the same parameters as in Fig.  9A but with larger recruitment demethylation rates (r AU and r RU ). The results are shown in Fig. 9C. Both of the corresponding distributions in Fig. 9A become narrower in Fig. 9C. In particular, the changes in the broad distribution (right panel) is particularly dramatic. This suggests that it may be easier to see changes in the broad bivalent domain than the narrow one during cell differentiation in experiments. Overall, our results demonstrate that nucleation sites can be responsible for the onset of bounded domains of AR nucleosomes. Also, narrow distributions can be obtained via either enhanced histone demethylation via exchange or via enhanced recruitment. We have also studied the distributions of AR nucleosomes, active marks, and repressive marks, using other reasonable parameter choices (see Figure S2 and the corresponding texts in Text S1). Taken together, these results suggest that highly localized bivalent domain patterns can be established surrounding nucleation sites, similar to the one-mark scenario described in previous study [11]. However, the local dynamics is more complex because multiple states are involved in the competition. The end configuration is an equilibrium resulting from the balance of multiple molecular forces.

States
During DNA replication, the nucleosomes, along with their associated histone marks, must be dissociated from the mother strand. How these marks are reassembled to the newly synthesized strands remains poorly understood. Recent studies suggest that the nucleosome, along with their associated marks, are randomly distributed to daughter strands [31]. In this section, we use our model to study the impact of DNA replication on the level of AR nucleosomes.
We choose parameters which correspond to cell environments during the formation of bivalent domains (see Fig. 10). Also, we assume that nucleation sites lose their properties at the very beginning of the simulations, so that there are no nucleation sites. We then vary the cell cycle lengths from 6 hours to 24 hours, which corresponds to varying the cell cycle length from that in stem cell to that in differentiated cells. We run the simulations for 10 cell cycles such that the average level of AR nucleosomes over a cell cycle reaches a stable value. Fig. 10 shows the average level of AR nucleosomes as a function of cell cycle length, where these levels are computed by averaging the number of AR nucleosomes at the end of the simulations over the lattice and over all simulation runs. Fig. 10 shows that the average level of AR nucleosomes is, in general, larger for longer cell-cycle. This result is expected, since there is more time for the lattice to recover from the loss of AR nucleosomes, caused by DNA replication, when the cell-cycle is longer. But the significance of cell-cycle length seems to be weaker for strong bivalent domains (blue curve in Fig. 10). This result is consistent with the experimental finding that higher levels of histone marking are observed when the length of the cell cycle increases [36].

Discussion
Development of computational models of bivalent domain dynamics can help to elucidate the mechanism of chromatin domain formation, and give insight for formulating and analyzing experimental studies. In this paper we introduce a model that incorporates multiple histone marks on a nucleosome and the interactions among these marks. We have illustrated the potential use of our model by employing it to investigate the dynamics of bivalent domains, with the following results.
Our main conclusion is that the formation of bivalent domains are highly stochastic at individual nucleosomes, but reproducible patterns can been obtained by averaging a large number of simulations. Dynamic changes of these patterns are maintained by the subtle balance of the multiple factors including the exchange rate, recruitment, distribution of nucleation sites, and cell-cycle length, resulting high degree of plasticity which might be advantageous for facilitating smooth transitions between cell-states during development.
Our analysis suggests that the formation of bivalent domains is in general a slow, two-step process, which can be divided into an   expansion and a stabilization phase. In contrast, the decay of bivalent domains, induced by demethylase activities, is much faster. This asymmetry between formation and decay dynamics may be an important feature for development control and perhaps needs to be taken into consideration into development epigeneticbased therapeutic approaches.
Specific epigenetic patterns can be established through targeted recruitment of chromatin regulators to specific genomic sequences. The effect of such nucleation sites on the establishment of highly localized epigenetic patterns has been studied via computational models in a number of previous studies [11,37]. We have extended these investigations by considering multiple histone marks in our model. As expected, we found that the strength of a nucleation site plays an important role in maintenance of localized bivalent domains. In the absence of nucleation sites, the bivalent domains either expands to the whole nucleosome array or disappears entirely. Our analysis is consistent with numerous experiment studies, which show that GC-rich DNA sequences are required for establishment of bivalent domains [26].
One limitation of our current model is that many kinetic parameters remain unknown, preventing us from making more quantitative predictions. Nevertheless, the major conclusions described above are robust with respect to parameter value changes therefore may reflect true biological principles. It will be interesting to test these principles by conducting quantitative experimental measurements. Figure S1 Illustration for the explanation of the states of the 6-state model. (TIFF) Figure S2 An example of the distribution of AR nucleosomes, active, and repressive marks. This plot illustrates that the 4-state model described in the main text can simulate bivalent domains (blue) in which the active mark (green) is less extensive than the repressive mark (red) (i.e., the bivalent domains (blue) are buried in the repressive domains (red)). The details of simulation can be referred to Section 4.3 in the main text. Here, distributions of AR nucleosomes (blue), ARzUR nucleosomes (red), and ARzAU nucleosomes (green) are plotted at the end of the simulation runs (time = 1800). The average levels of nucleosomes are averaged over 1000 simulation runs. In the simulation, p i~40 UA~0 :03 and p i~40 UR~0 :015. The other parameters are r UA~0 :029, r UR~0 :021, p AU~0 :025, p RU~0 :015, r AU~0 :004 and r RU~0 :002.

(TIFF)
Text S1 The supplementary text provides details regarding the formulation of our model and another example of localization of AR states related to our results in Section 4.3. (PDF) Figure 10. Average AR nucleosome level vs. cell-cycle length. The average AR nucleosome level is plotted as a function of cell-cycle length. These levels of AR nucleosomes are computed by averaging the final number of AR nucleosomes in the simulations over all the runs and the whole lattice. r UA~2 r UR~0 :046, p UA~pUR~0 , and r AU~rRU~0 . doi:10.1371/journal.pone.0077944.g010