Mathematical modeling of histone modifications reveals the formation mechanism and function of bivalent chromatin.

Bivalent chromatin is characterized by occupation of both activating and repressive histone modifications. Here, we develop a mathematical model that involves antagonistic histone modifications H3K4me3 and H3K27me3 to capture the key features of bivalent chromatin. Three necessary conditions for the emergence of bivalent chromatin are identified, including advantageous methylating activity over demethylating activity, frequent noise conversions of modifications, and sufficient nonlinearity. The first condition is further confirmed by analyzing the existing experimental data. Investigation of the composition of bivalent chromatin reveals that bivalent nucleosomes carrying both H3K4me3 and H3K27me3 account for no more than half of nucleosomes at the bivalent chromatin domain. We identify that bivalent chromatin not only allows transitions to multiple states but also serves as a stepping stone to facilitate a stepwise transition between repressive chromatin state and activating chromatin state and thus elucidate crucial roles of bivalent chromatin in mediating phenotypical plasticity during cell fate determination.

Introduction

A variety of covalent modifications on histones, such as methylation and acetylation, play important roles in governing the chromatin state. Such diversity of histone modifications confers additional complexity of transcription regulation (Kouzarides, 2007; Lawrence et al., 2016; Li et al., 2007). An appealing biological phenomenon is the co-occurrence of activating histone modifications and repressive histone modifications at the same chromatin region, namely, bivalent chromatin (Bernstein et al., 2006). Bivalent chromatin is initially reported to mark development-related genes in embryonic stem cells, and further studies validate its existence in lineage-committed progenitor cells, tumor cells, and other cell types (Azuara et al., 2006; Bernstein et al., 2006; Chaffer et al., 2013; Kinkley et al., 2016; Matsumura et al., 2015). A prominent bivalent chromatin form is composed of both the activating modification H3K4me3 and the repressive modification H3K27me3 (trimethylation of lysine 4 or lysine 27 on histone H3), which are catalyzed by the Trithorax group (TrxG) proteins and Polycomb group (PcG) proteins, respectively (Schuettengruber and Cavalli, 2009; Schuettengruber et al., 2011). Other forms of bivalent chromatin such as H3K4me3/H3K9me3 and H3K9me3/H3K36me3 are also observed (Matsumura et al., 2015; Mauser et al., 2017). Such bivalent chromatin is found to be associated with cellular phenotypical plasticity during cell differentiation and oncogenesis (Bernstein et al., 2006; Chaffer et al., 2013; Matsumura et al., 2015; Rugg-Gunn et al., 2010). For example, the bivalent chromatin is viewed as a poised state with bidirectional potential by maintaining key developmental genes at low expression levels and waiting for further activation or repression upon appropriate signals (Bernstein et al., 2006).

However, several key issues of the bivalent chromatin remain controversial. First, the pattern formation mechanism of bivalent chromatin is still unclear. Intuitively, bivalent chromatin results from the balance between TrxG proteins and PcG proteins (Voigt et al., 2013), as suggested by the observations that mutations or experimental perturbations tilting such balance dramatically reduce bivalent histone modifications (Béguelin et al., 2013; Shema et al., 2016). How the bivalent chromatin emerges from the interaction between TrxG and PcG still remains unknown. Second, a quantitative analysis of bivalent chromatin composition at the nucleosome level is lacking. Most of the bivalent chromatin domains are identified by using separate single-antibody chromatin immunoprecipitation (ChIP) experiments for each of the modifications. Consequently, it is not sufficient to depict the modification state of each nucleosome at the chromatin region, and it is unclear whether the two marks reside simultaneously on the same nucleosome. Although genome-wide sequential ChIP or other similar techniques like co-ChIP and re-ChIP can detect bivalent nucleosomes carrying both two marks (Kinkley et al., 2016; Sen et al., 2016; Weiner et al., 2016), the extent to which the bivalent nucleosomes account for the bivalent chromatin remains to be investigated. Third, the causal role of bivalent chromatin in mediating phenotypical plasticity remains largely controversial. Several experimental studies suggest that many bivalent chromatin domains are resolved upon differentiation and the expression levels of the corresponding genes change as well (Bernstein et al., 2006). Although such coincidence establishes the association between the bivalent chromatin and its proposed function, the causality is difficult to confirm because perturbations of bivalent chromatin inevitably lead to the change of a broad range of genes, making it difficult to interpret the observations (Sneppen and Ringrose, 2019; Voigt et al., 2013).

Previous mathematical models of histone modifications mainly focus on the emergence of epigenetic bistability and spatial patterns as well as heritability (Dodd et al., 2007; Zhang et al., 2014). The theoretical works on modeling the bivalent chromatin state are still limited and oversimplified. For example, when modeling the H3K4me3/H3K27me3 bivalent chromatin, the intermediate, yet critical, modification steps such mono- and dimethylations are usually neglected (Huang and Lei, 2017; Ku et al., 2013; Sneppen and Ringrose, 2019), which have been implicated to dramatically affect the system behavior by increasing nonlinearity (Berry et al., 2017). A comprehensive model has been developed to incorporate four kinds of histone modifications; however, details like the mutual exclusiveness of H3K4me3 and H3K27me3 on the same histone tail as well as different nucleosome interaction ranges are ignored, which may be critical to the bivalent dynamics (Sneppen and Ringrose, 2019).

In this article, we develop a minimal yet fully equipped mathematical model to imitate the stochastic dynamics of two antagonistic histone modifications H3K4me3 and H3K27me3. Our model incorporates the essential elements, including the intermediate steps of histone modifications, the mutual exclusiveness of H3K4me3 and H3K27me3 on the same histone tail, and different nucleosome interaction ranges. To decipher the formation mechanism of bivalency, we investigate the parameter constraints that yield the bivalent chromatin state and analyze the biological importance of those constraints. We demonstrate that advantageous methylating activity over demethylating activity, noisy conversions of modifications, and nonlinearity are required for the emergence of bivalent chromatin. To probe the relationship between bivalent nucleosome and bivalent chromatin, we analyze the composition of bivalent chromatin at the nucleosome level under different nucleosome interaction regimes and predict that no more than 50% of nucleosomes at the bivalent chromatin domain are bivalent nucleosomes. Finally, by analyzing the transitions among chromatin states in response to external stimulus or noise, we find that bivalent chromatin allows the transition to multiple states and can serve as a stepping stone to facilitate a stepwise transition between repressive chromatin state and activating chromatin state, showing the potential role of bivalent chromatin in mediating phenotypical plasticity during cell fate decision.

Results

Modeling H3K4me3/H3K27me3 bivalent chromatin

To capture the details of histone modification patterns, we propose a minimal model while incorporating all the essential interaction elements. Inspired by a previous work (Berry et al., 2017), which focuses on H3K27 methylations for modeling the histone modification dynamics, our model is constructed to study the interplay between two distinct kinds of histone modifications, i.e., H3K4 methylations and H3K27 methylations. The basic settings of our model are as follows:

Histone modifications. In our model, a linear array of nucleosomes as well as their interactions, representing a chromatin region in the cell, is considered (Figure 1A). Although each nucleosome actually is an octamer composed of two copies of four core histones (H2A, H2B, H3, and H4), we focus on the modifications on H3 histone and assume that a nucleosome consists of two copies of H3 histones by neglecting the other three types of histones for simplification. Each histone H3 is assumed to have one modification site and can be in one of the seven modification states: K27me3, K27me2, K27me1, Empty, K4me1, K4me2, and K4me3 (Figure 1B). Here “Empty” means no methylation on histone H3. “Implicit” in the above assumption is that H3K27 methylations and H3K4 methylations are mutually exclusive on the same histone tail, owing to the experimental observation that H3K27 methylations and H3K4 methylations can hardly exist on the same histone tail (Shema et al., 2016; Voigt et al., 2012).

Figure 1

A model of histone modification at bivalent chromatin domain

(A) Schematic illustration of the model for histone modification dynamics.

(B) The modification state of a histone H3 and feedbacks mediated by H3K4me3 and by H3K27me3. Black and gray arrows indicate methylation and demethylation reactions, respectively.

(C) An illustration of histone “neighbors” for local interaction regime and global interaction regime. For the local interaction regime, the neighbors of the fourth histone are histone numbers 1 2, 3, 5, and 6 (indicated by the red dashed box); for the global interaction regime, the neighbors include all the other N-1 histones (indicated by the gray dashed box).

(D) Classification of the chromatin states. Four chromatin states are defined according to both the number of H3K4me3 and the number of H3K27me3. In the left plot, the region for the unmodified state is bounded by the curve and the axes in the first quadrant. The boundary between the regions for the bivalent state and the activating state is y = x/2, and the boundary between regions for the bivalent state and the repressive state is y = 2x. The right panel illustrates the possible composition for four chromatin states at the nucleosome level.

Reactions. For each histone H3, the modification state can change toward neighboring states through one-step reactions (Figure 1B). Although there are 12 possible reactions in total, at most two possible reactions can occur for a histone with a specific modification state. For example, if the current modification state of the ith histone (i = 1, 2, …, N, where N is the total number of histones) is H3K4me2, then it can be either methylated to H3K4me3 or demethylated to H3K4me1. The reaction rate depends on not only the background catalytic activity of histone modification enzymes but also the modifications on neighboring nucleosomes. For instance, existing H3K27me3 can facilitate H3K27 methylations on neighboring nucleosomes by recruiting the enzymatic complex PRC2 (polycomb repressive complex 2) (Hansen et al., 2008; Margueron et al., 2009), and also enhance H3K4 demethylation by interacting with KDM5A through PRC2 (Pasini et al., 2008), thus forming a positive feedback for H3K27me3 deposition (Figure 1A). Similarly, a recruitment-mediated positive feedback for H3K4me3 deposition is also suggested by previous studies (Agger et al., 2007; Dou et al., 2006; Issaeva et al., 2007; Lee et al., 2007; Ruthenburg et al., 2007; Voigt et al., 2013; Wysocka et al., 2005). Thus, our model incorporates the feedbacks mediated by H3K4me3 and H3K27me3, by assuming that H3K4me3 (or H3K27me3) promotes all the six reactions toward H3K4me3 (or H3K27me3) on neighboring nucleosomes (Figure 1B).

Range of feedback interactions. Different nucleosome interaction ranges are considered when modeling the feedbacks from neighboring nucleosomes. In our model, by specifying the range of “neighbors” of each histone along with the assumption that only neighbors can interact with each other, we can establish different nucleosome interaction regimes. Two typical interaction regimes are studied: the local interaction regime where the neighbors of one histone include the other histone on the same nucleosome and histones in adjacent nucleosomes, and the global interaction regime where the neighbors include all the other histones regardless of their distance in the one-dimensional chromatin fiber (Figure 1C).

Reaction rate formulas. Based on assumptions (1)–(3), the formulas used to model all 12 possible reactions for the ith histone are presented in Table 1. Although the reaction rates for different degrees of lysine methylation (i.e., me0→me1, me1→me2, and me2→me3) may be naturally different (McCabe et al., 2012), to reduce the number of parameters used in the model, the reaction rates for the three methylation reactions on H3K4 (or H3K27) are assumed to be the same. Similarly, the reaction rates for the three demethylation reactions (i.e., me3→me2, me2→me1, and me1→me0) on H3K4 (or H3K27) are set to be the same. Each reaction rate contains a feedback term that reflects recruitment-mediated (de)methylation and a noise term representing noisy (de)methylation resulting from background enzyme catalytic activity. For example, for the ith histone carrying H3K4me2, the methylation rate and demethylation rate can be formulated as and , respectively. Here, E and F are defined as fractions of H3K4me3 and H3K27me3 among the neighbors of the ith histone, respectively. The methylation rate consists of a noise term and a feedback term , where is the recruitment-mediated H3K4 methylation rate constant and reflects the positive feedback mediated by neighboring H3K4me3 on the deposition of H3K4 methylations. Likewise, the demethylation rate contains a noise term and a feedback term representing the inhibitory effect of neighboring H3K27me3 on H3K4me3 by enhancing the removal of H3K4 methylations. Functions f and g are chosen as Hill functions to consider the nonlinear feedback effect. The rate formulas for all histones are constructed in a similar way.

Table 1

Rate formulas of all 12 reactions for the ith histone

Reactions	Reaction rate
Empty → H3K4me1, H3K4me1 → H3K4me2, H3K4me2 → H3K4me3	k_K4mef(E_i)+α_K4me
H3K4me3 → H3K4me2, H3K4me2 → H3K4me1, H3K4me1 → Empty	d_K4meg(F_i)+β_K4me
Empty → H3K27me1, H3K27me1 → H3K27me2, H3K27me2 → H3K27me3	k_K27meg(F_i)+α_K27me
H3K27me3 → H3K27me2, H3K27me2 → H3K27me1, H3K27me1 → Empty	d_K27mef(E_i)+β_K27me

See also Table S1 for parameters.

DNA replication. The random loss of histone modifications during DNA replication is also incorporated in our model. We assume that DNA replication occurs once per cell cycle. When DNA replication occurs, the parental nucleosomes are randomly partitioned and redistributed into two daughter cells and mixed with newly synthesized, unmodified nucleosomes. In simulation we track one daughter cell, and upon DNA replication, each nucleosome is replaced by an unmodified nucleosome with a probability of 0.5.

The stochastic simulations of our model are performed by using the Gillespie algorithm (Gillespie, 1977) (see STAR Methods for details). Besides, we develop a deterministic model based on the reaction rate formulas (see STAR Methods). This model will be used to obtain steady-state solutions and provide a guide for the stochastic simulations. All parameters in our model are shown in Table S1. There are four free parameters in total, including , , , and . Since we fix recruitment-mediated demethylation rate constants as ones, i.e., , parameters kK4me and kK27me can be interpreted as the ratio of recruitment-mediated methylation rate and demethylation rate for H3K4 and H3K27, respectively, also called methylation-to-demethylation ratios. Parameters and are noise-to-feedback ratios, defined as the ratio of noisy (de)methylation rate to recruitment-mediated (de)methylation rate constant for H3K4 and H3K27, i.e., , , respectively.

Emergence of the bivalent chromatin state

By directly simulating the system with the Gillespie algorithm (Gillespie, 1977), we obtain the temporospatial dynamics of histone modifications. The chromatin states can be evaluated based on the abundance of H3K4me3 and H3K27me3 over the chromatin region. We classify the chromatin states into four categories, including activating chromatin state, repressive chromatin state, bivalent chromatin state, and unmodified state (Figure 1D). If both H3K4me3 number and H3K27me3 number are at low levels, i.e., most of the histones are mono/dimethylated or unmodified, such chromatin state is named as the unmodified state. Otherwise, the chromatin is regarded as the activating state if the number of H3K4me3 is twice larger than that of H3K27me3, the repressive state if the number of H3K27me3 is twice larger than that of H3K4me3, or the bivalent state if the numbers of H3K4me3 and H3K4me27 are balanced (Figure 1D).

Depending on the parameter choice, the four chromatin states can either exist individually or coexist with others. In particular, our model can exhibit three distinct modes: monostability, bistability, and tristability. For example, there exists only the repressive chromatin state when is much greater than (Figures 2A–2C), and the system can achieve bistability switching between the activating state and the repressive state when is close to (Figures 2D–2F). Besides, within a certain range of parameters, the bivalent chromatin state can emerge in a tristable mode, coexisting with the activating state and the repressive state. Simulation of the stochastic model and analysis of the corresponding deterministic model both confirm such tristability (Figures 2G–2I). It is worth noting that, unlike the case under the global interaction regime (Figure 2), stochastic simulations show that there is no clear bistable or tristable behavior under the local interaction regime because the boundaries between H3K4me3 patches and H3K27me3 patches would randomly walk along the chromatin region (Figure S1).

Figure 2

Typical simulations of histone modification dynamics and emergence of the bivalent chromatin state

(A–C) The dynamics (A), landscape (B), and phase plane (C) of a system exhibiting the monostable repressive chromatin state. Parameters: , , .

(D–F) The dynamics (D), landscape (E), and phase plane (F) of a system operating in a bistable mode where activating state and repressive state coexist. Parameters: , , .

(G–I) The dynamics (G), landscape (H), and phase plane (I) of a tristable system where activating state, repressive state, and bivalent state coexist. Parameters:  = , . The two-dimensional chromatin state landscape (B, E, and H) is obtained by the density of stochastic trajectories, and the two-dimensional phase plane (C, F, and I) is obtained by the deterministic trajectories (gray lines), where red solid dots and red circles represent stable steady states and unstable steady states, respectively. The arrows indicate the evolving directions of trajectories. The global interaction regime is adopted in all stochastic simulations.

See also Figures S1 and S2.

We also perform a simulation of the deterministic model with 104 sets of randomly assigned parameters and identify a total of 16 types of combination patterns of the four chromatin states (Figure S2). In addition to the tristable mode, the bivalent state can be observed in two bistable modes, coexisting with either the activating state or the repressive state; however, the parameter space for either bistable mode is smaller than that for the tristable mode (Figure S2B).

Formation mechanisms of the bivalent chromatin state

High methylation-to-demethylation ratio and high noise-to-feedback ratio

We next investigate the parameter space to answer what conditions are required for generating the bivalent chromatin state. First, we set the noise-to-feedback ratios for H3K4 and H3K27 as the same, i.e., . For a given level of noise-to-feedback ratio γ, we next calculate the time-averaged probabilities (i.e., the fraction of time spanned by a chromatin state) of four chromatin states under different combinations of and (Figures 3A and 3B). We find that, under the global interaction regime, for most of the parameter combinations, the system exhibits monostable mode or bistable mode (Figure 3A), a result consistent with previous studies (Berry et al., 2017; Dodd et al., 2007; Sneppen and Ringrose, 2019). Nevertheless, we still find that a relatively small parameter space exists for the bivalent state. When the methylation-to-demethylation ratios for H3K4 and H3K27 are well matched and both large enough, the bivalent state is usually observed in a tristable mode where the activating state and the repressive state as well as the bivalent state coexist (Figures 3A and S3A–S3C). The requirement of high methylation-to-demethylation ratio (i.e., advantageous methylating activity over demethylating activity) for the bivalent state is further confirmed by analyzing the deterministic model (Figure S3D). Although there is no clear bistable or tristable behavior under the local interaction regime, we can still calculate the time-averaged probability of each chromatin state (Figure 3B) and observe a similar pattern of parameter space for the four chromatin states to that under the global interaction regime. Compared with the global interaction regime, the local one shows a lower probability of the bivalent state but allows the bivalent state under smaller methylation-to-demethylation ratios (Figures 3A and 3B). In addition to the methylation-to-demethylation ratio, the noise-to-feedback ratio plays a role in the emergence of the bivalent state. A high noise-to-feedback ratio, which reflects frequent noise conversions of histone modifications, leads to the large probability of the bivalent state; with the increase of noise-to-feedback ratio, the monostable bivalent state may emerge (see also Figure S3C). These results also hold for a small system size N = 20 (Figures S4A and S4B). We also investigate the effect of cell cycle length on bivalent state (Figures S4C–S4F). Despite the fact that frequent cell division could induce random loss of histone modifications and may destabilize chromatin states, it contributes to a high probability of bivalent chromatin state (Figures S4C and S4D). The reason is that the activating state and repressive state experience a more dramatic decrease of the mean dwell time than the bivalent state, as the cell cycle length becomes shorter (Figures S4E and S4F).

Figure 3

The effects of methylation-to-demethylation ratio and noise-to-feedback ratio on the emergence of bivalent chromatin state

(A and B) Parameter space for four different chromatin states under global interaction regime (A) and local interaction regime (B). For each interaction regime, time-averaged probabilities of four chromatin states under different combinations of kK4me and kK27me are calculated. For each parameter set, 100 simulations with randomly assigned initial values are used for averaging; for these initial values, numbers of H3K4me3 and H3K27me3 obey a discrete uniform distribution on the region and the other histones are set to be Empty.

(C) The bifurcation diagram for the steady-state H3K4me3 number as a function of the methylation-to-demethylation ratio. Solid lines and dashed lines represent stable steady states and unstable steady states, respectively.

(D) The maximum of the eigenvalue real parts of the Jacobean matrix for the intermediate steady state (indicated by the black line in (C)).

(E and F) The same plots as (C-D) except that the noise-to-feedback ratio γ is set to be the control parameter. Parameters for (C-D) are set to be , whereas those for (E-F) are set to be .

(G) The number of bivalent genes (upper panel, extracted from Extended Data Figure 8B of a previous study [Liu et al., 2016]) and the relative protein expression level of H3K4me3 and H3K27me3 writers (lower panel, extracted from Extended Data Figure 2D of a previous study [Liu et al., 2016]) in different embryo stages. The number of bivalent genes for the blastocyst stage embryo is estimated by the average of the numbers for the inner cell mass and trophectoderm. The median values of the relative protein expression levels of H3K4me3 and H3K27me3 writers are used and then normalized to the maximum among different stages.

(H) The fold change of the number of bivalent genes and the relative protein expression levels of H3K4me3 and H3K27me3 writers between adjacent stages.

See also Figures S3 and S4.

By utilizing the deterministic model, we analyze how the methylation-to-demethylation ratio and the noise-to-feedback ratio affect the bivalent chromatin state. As shown in the bifurcation diagram in Figure 3C, the bivalent state derives from a branch between the activating state branch and the repressive state branch, i.e., an intermediate state. When the methylation-to-demethylation ratio () is small, such intermediate state behaves as an unmodified state since the number of H3K4me3 (and H3K27me3) is close to zero. As k increases, the bivalent state emerges owing to the increased amounts of H3K4me3 and H3K27me3. We also notice the non-monotonic effect of the methylation-to-demethylation ratio on the stability of the bivalent state: the intermediate state becomes stable once k exceeds some threshold, but the stability will reach its maximum and then decline as k continues to increase (Figure 3D). When the methylation-to-demethylation ratio is fixed, increasing the value of noise-to-feedback ratio γ has no effect on the H3K4me3 (or H3K27me3) number of the bivalent state but can stabilize the bivalent state (Figures 3E and 3F). Thus, we conclude that both the high methylation-to-demethylation ratio and the high noise-to-feedback ratio collectively contribute to the emergence of the bivalent chromatin state: the high methylation-to-demethylation ratio pushes the system away from the unmodified state but toward the bivalent state, whereas the high noise-to-feedback ratio stabilizes the bivalent state.

A recent experimental study shows the positive role of a high methylation-to-demethylation ratio in generating the bivalent chromatin state (Liu et al., 2016). Liu et al. performed a genome-wide analysis of H3K4me3 and H3K27me3 in mouse pre-implantation embryos and observed that the number of bivalent genes changes dynamically over different embryo developmental stages (the upper panel of Figure 3G). The relative protein expression levels of H3K4me3 and H3K27me3 writers are measured (the lower panel of Figure 3G). By analyzing how the number of bivalent genes and expression levels of H3K4me3 and H3K27me3 writers change over adjacent stages, we find that simultaneous elevation of expression levels of antagonistic modification writers is associated with the increased number of bivalent genes (i.e., four-cell stage to eight-cell stage and Morula stage to Blastocyst stage in Figure 3H). This is consistent with our finding that the higher methylation-to-demethylation ratios lead to the more probable bivalent chromatin state.

Nonlinearity

As mentioned above, the bivalent chromatin state usually emerges along with the tristability (Figure S2B), which typically results from the nonlinearity. The nonlinearity in our model arises not only from the Hill functions (versus linear functions) to describe the feedback but also by incorporating multiple positive feedback loops in the reactions converting between H3K4me3 and H3K27me3. For example, the overall rate for consecutive conversions from H3K27me3 to H3K4me3 is proportional to , a nonlinear function naturally generated by the multiple positive feedback loops mediated by H3K4me3.

Here, we show that sufficient amount of nonlinearity is necessary to generate the chromatin tristability, leading to the bivalent state. To elucidate the effect of nonlinearity, we constructed 11 variants of the standard model with weakened nonlinearity by reducing intermediate histone modification states, eliminating the feedback on either methylation reactions or demethylation reactions, or using linear feedback (Figure 4A), then investigated their robustness of generating the desired tristability (i.e., coexistence of activating, repressive and bivalent states; see also the inset of Figure 4B). For each of the 12 models (one standard model plus 11 variants), the robustness is measured by Q value (Qiao et al., 2019), defined as the number of parameter sets that yield the specified tristability among 10,000 sets of randomized parameters by simulating the deterministic dynamics. Compared with the standard model, all these variants show decreased Q values, i.e., reduced robustness of generating tristability with a bivalent state, and only the three seven-state models with nonlinear feedback achieve non-zero Q values (Figure 4B). Furthermore, the system with feedback on only demethylation reactions shows less robust tristability than that with feedback on only methylation reactions (Figure 4B), which is consistent with their robustness of bistability (Dodd et al., 2007). Besides, we notice that some variants with linear feedback can also generate tristability; however, the intermediate state of the generated tristability is not the bivalent state but an unmodified state (Figures 4C and S5), which indicates a special need for nonlinear feedback for the emergence of the bivalent state. In addition, enhancing nonlinearity by increasing the Hill coefficient within a certain range dramatically improves robustness of generating the desired tristability for the standard model (Figure 4D). Since high Hill coefficients usually reflect strong cooperativity, our result also indicates that high cooperativity is required for robust tristability. Here, cooperativity of histone modifications means that the production rate of one kind of modification (i.e., H3K4me3) exhibits a sigmoid dependence on the abundance of existing modifications. Taken together, these results illustrate the key role of nonlinearity in generating tristability with a bivalent state.

Figure 4

The role of nonlinearity in generating the bivalent state

(A) The standard model and its variants with reduced nonlinearity. These models include the (a) seven-state model with feedbacks on both methylation and demethylation, (b) seven-state model with feedbacks on methylation only, (c) seven-state model with feedbacks on demethylation only, (d) three-state model with feedbacks on both methylation and demethylation, (e) three-state model with feedbacks on methylation only, and (f) three-state model with feedbacks on demethylation only. Each of (a-f) can use either nonlinear feedback or linear feedback, thus forming 12 models in total. For the linear feedback, we choose .

(B and C) Q values of 12 models (including the standard model and 11 variants) for generating tristability with a bivalent state (B) or with an unmodified state (C). Each inset indicates the pattern of stable steady states on the two-dimensional phase plane for the specified tristability mode. For each model, we use the Latin hypercube method to sample parameters from log-uniform distributions with the parameter ranges [0.1,10] for and and [0.01 100] for and .

(D) The Q value as a function of the Hill coefficient for the standard model. Here, Hill functions f and g have the form of and , respectively, where n is the Hill coefficient.

See also Figure S5.

Composition of bivalent chromatin at the nucleosome level

Although the chromatin region exclusively occupied by bivalent nucleosomes that contain H3K4me3 and H3K27me3 on the opposing H3 histone tails is obviously bivalent, the inverse is not always true. To clarify the relationship between the bivalent chromatin and the bivalent nucleosome, we investigate the composition of the bivalent chromatin state at the nucleosome level. We find that the proportion of bivalent nucleosomes at the bivalent chromatin depends on what kind of interaction regime is applied (Figure 5). Under the global interaction regime, the bivalent state shows the pattern of randomly arranged H3K4me3 and H3K27me3 along the chromatin (Figure 5A), and the numbers of nucleosomes carrying H3K4me3/H3K4me3, H3K4me3/H3K27me3, H3K27me3/H3K27me3 at the bivalent chromatin are approximately in the ratio (Figures 5B and 5C). Furthermore, such ratio could be explained by a trinomial distribution with two parameters and , which represent probabilities for a histone to be H3K4me3 and H3K27me3, respectively (Figure S6). However, for the local interaction regime, the composition of the bivalent state is different: the chromatin region is composed of balanced H3K4me3 patches and H3K27me3 patches (Figure 5D), and the ratio of the three kinds of nucleosomes becomes , where (Figures 5E and 5F), which reflects a reduction of bivalent nucleosome proportion compared with the case under the global interaction regime. Intuitively, under the local interaction regime, existing K4me3 (or K27me3) promotes deposition of new K4me3 (or K27me3) on only adjacent histones, and thus, the formation of K4me3 (or K27me3) patches, which hinder the formation of bivalent nucleosomes; however, under the global interaction regime, newly added K4me3 (or K27me3) does not necessarily reside around the existing K4me3 (or K27me3), which benefits the formation of bivalent nucleosomes. Furthermore, the proportion of bivalent nucleosomes is showed to be an increasing function of the noise-to-feedback ratio γ, and the upper bound 50% can be approached when extremely large γ is used (see the inset of Figure 5F).

Figure 5

The composition of bivalent chromatin state at the nucleosome level under different interaction regimes

(A) The typical dynamics under global interaction regime. Upper panel: the temporospatial dynamics of histone modification. Lower panel: the dynamics of H3K4me3 number and H3K27me3 number.

(B) The time-averaged percentages of different nucleosome types at the bivalent chromatin state. Only percentages of nucleosomes with H3K4me3/H3K4me3, H3K4me3/H3K27me3, and H3K27me3/H3K27me3 are shown. The data are shown as mean ± SD for 100 repeated and independent simulations.

(C) The relationship between H3K4me3/H3K27me3 nucleosome percentage and (the square root of) the product of H3K4me3/H3K4me3 and H3K27me3/H3K27me3 nucleosome percentages.

(D-F) The same plots as (A-C) except that local interaction regime is applied. Inset of (F): the percentage of bivalent nucleosomes at the bivalent chromatin under local interaction regime as a function of noise-to-feedback ratio. For (A-B) and (D-E), parameters are set to be and . For (C) and (F), each data point represents a parameter set with randomly assigned and as well as a fixed γ, which yields the bivalent state. For the inset of (F), we set and vary from −2 to 2 with an increment 0.1.

See also Figures S6 and S7.

Taken together, an important finding is that the maximum proportion of bivalent nucleosomes at the bivalent chromatin is (or near considering stochasticity due to the finite system size N) regardless of nucleosome interaction regimes. Moreover, can reach its maximum value 50% when methylation-to-demethylation ratios are far greater than one and symmetric parameters are adopted. In particular, when methylation-to-demethylation ratios are far greater than one, histone modifications at bivalent chromatin are almost either H3K4me3 or H3K27me3, i.e., ; symmetric parameters further lead to . Thus, can reach a maximum of 50%.

It should be noted that the fact that the proportion of bivalent nucleosomes cannot exceed 50% in our model results from mutual exclusiveness of K4 methylations and K27 methylations on the same H3 histone tail. Because of such mutual exclusiveness, a nucleosome can be in one of the only four modification states, i.e., H3K4me3/H3K4me3, H3K4me3/H3K27me3, H3K27me3/H3K27me3, and Empty/Empty, if we neglect the mono- and dimethylations. Then high methylation-to-demethylation ratio, required for the emergence of bivalent chromatin, will increase the proportions of all kinds of nucleosomes except Empty/Empty nucleosomes. Thus, H3K4me3/H3K27me3 nucleosomes must coexist with H3K4me3/H3K4me3 nucleosomes and H3K27me3/H3K27me3 nucleosomes. However, the case without such mutual exclusiveness could yield a bivalent state with the percentage of bivalent nucleosomes more than 50% and even near 100%. We remove the mutual exclusiveness in a modified model by assuming that each H3 histone has two separate modification sites, one for K4 methylations and the other for K27 methylations (Figures S7A and S7B). In this case, the high methylation-to-demethylation ratio leads to a bivalent chromatin state, where near all H3K4 and H3K27 modification sites are occupied by H3K4me3 and H3K27me3, respectively; in other words, almost all nucleosomes are bivalent (Figure S7C).

Bivalent chromatin as an intermediate state to mediate cell plasticity

Finally, we explored the functional roles of bivalent chromatin in mediating phenotypical plasticity during cell fate determination. To this end, we first investigated whether the bivalent chromatin state allows transitions to multiple states in response to either external stimulus or noise. To analyze how the bivalent state responds to the external stimulus, we performed a bifurcation analysis of the deterministic model and found that the bivalent state can make transitions to other states as the control parameter (i.e., the external stimulus) changes (Figure 6A). When the methylation strength of H3K4 surpasses that of H3K27, the bivalent state switches to the activating state; when H3K4 methylation strength declines, it can switch to the repressive state. In the parameter regime of tristability without external stimulus, noise-driven transitions among the three chromatin states may spontaneously arise. For example, the noise caused by stochastic chemical reactions, and the random loss of histone modifications during DNA replication, is found to be sufficient in triggering the transition from the bivalent state to either the activating state or the repressive state (Figure 5A). Furthermore, such plasticity of the bivalent state is found to be enhanced by cell division. As division occurs more frequently, the dwell time of the bivalent state becomes shorter, benefiting the switching to other states (Figure 6B). This result is significant because the division-enhanced plasticity only takes place for a subset of the intermediate states. For example, in the seven-state model with feedbacks on methylation only and with the linear feedback, there exists a tristable mode where the unmodified state serves as the intermediate state (Figures 4C and S5); however, the dwell time of the unmodified state increases as the cell division becomes faster (Figure 6C).

Figure 6

Bivalent chromatin as an intermediate state to mediate cell plasticity

(A) The bifurcation diagram of H3K4me3 number as a function of parameter kK4me. Parameters: and . The green arrows illustrate the transition direction of the bivalent state as kK4me changes.

(B and C) The mean dwell time of the bivalent chromatin state (B) and the unmodified chromatin state (C) under the condition of no cell division, period = 40, period = 20, and period = 10. The data are shown as mean ± SD for multiple independent simulations, and the SD is obtained by using the bootstrap method. Parameters for (B) are set to be and ; for (C), the seven-state model with feedbacks on methylation only and linear feedback (see Figures 4Ab and S5) is used, and parameters are set to be and .

(D and E) The MAP from the repressive chromatin state to the activating chromatin state for the parameter regime of tristability (D) and bistability (E). Inset: the accumulated action from the repressive chromatin state to the activating chromatin state along the relative coordinate ([0,1]) of the MAP. The MAP is calculated by using the geometric minimum action method (Heymann and Vanden-Eijnden, 2008), and 101 points are used in the discretization of the path. The dots in red, blue, and black represent activating, repressive, and bivalent states, respectively; circles represent unstable steady states. Parameters: , for (D) and for (E).

In addition to the multi-directional transition potential, the bivalent chromatin can serve as an intermediate state, facilitating the noise-driven transitions between the repressive state and the activating state. We used the minimum action method (E et al., 2004; Heymann and Vanden-Eijnden, 2008; Yu et al., 2018) to compute the most probable path, i.e., minimum action path (MAP), for the transition. The “action” used here refers to Wentzell-Freidlin action, which measures the level of “difficulty” of a transition; the most probable transition path could be obtained by minimizing such action: the less the action, the more probable the path. Of interest, the MAP from the repressive state to the activating state passes through the bivalent state, indicating that the transition is a two-step process (Figure 6D). As a result, the transition to the activating state from the bivalent state is easier than that from the repressive state, since the required action for the former is less than the latter (see the inset of Figure 6D). Moreover, by comparing the entire action of the MAP for the cases with and without the bivalent state (i.e., tristability versus bistability), we found that the existence of the bivalent state can decrease the least action required for the transition between the two extreme chromatin states and, thus, makes the transition more probable (Figures 6D and 6E). These results indicate that the bivalent state serves as a stepping stone to facilitate the stepwise transition between the two extreme chromatin states.

Discussion

Bivalent chromatin has been associated with many important biological processes including cell differentiation and oncogenesis. However, the underlying mechanisms of the bivalent chromatin remain unclear. For instance, how does the bivalent chromatin emerge from the interaction between PcG proteins and TrxG proteins? To what extent do the bivalent nucleosomes account for the bivalent chromatin? In addition, could bivalent chromatin have a causal role in mediating phenotypical plasticity?

To address these questions, we developed a minimal but fully equipped model to perform a quantitative study of two antagonistic histone modifications H3K4me3 and H3K27me3. We identified that the high methylation-to-demethylation ratio and high noise-to-feedback ratio collectively contribute to the emergence of bivalent chromatin: high methylation-to-demethylation ratio pushes the system away from the unmodified state but toward the bivalent state, while high noise-to-feedback ratio stabilizes the bivalent state. Our finding is supported by a recent experimental observation that the increased number of bivalent genes is associated with simultaneous elevation of the expression levels of H3K4me3 and H3K27me3 writers (Liu et al., 2016). Dodd et al. have noticed that a high noise-to-feedback ratio can increase the probability of intermediate states when studying the requirements for chromatin bistability (Dodd et al., 2007), which is consistent with the enhanced bivalent state stability induced by a high noise-to-feedback ratio. Besides, sufficient nonlinearity involved in our model, which has proved to be important for generating bistability (Berry et al., 2017; Dodd et al., 2007), could also extend the parameter space for generating multi-stability containing the bivalent state.

We quantitatively studied the composition of bivalent chromatin at the nucleosome level under different nucleosome interaction regimes: under the global interaction regime, proportions of different types of nucleosomes obey a trinomial distribution; under the local interaction regime, the proportion of bivalent nucleosomes shrinks and is an increasing function of noise-to-feedback ratio. Of most importance, we predicted that bivalent nucleosomes account for no more than 50% of nucleosomes at the bivalent chromatin domain regardless of nucleosome interaction regime, which results from mutual exclusiveness of K4 methylations and K27 methylations on the same H3 histone tail. This result agrees with the scattered distribution of H3K4me3/H3K27me3 bivalent nucleosomes around the transcription start site of some bivalent genes (Sen et al., 2016). This could also explain why some of the previous models find it difficult to produce a chromatin state with dominant bivalent nucleosomes (Huang and Lei, 2017; Ku et al., 2013; Sneppen and Ringrose, 2019). However, for the model without such mutual exclusiveness, bivalent nucleosomes can be dominant with some parameter settings, which is consistent with a recent study (Sneppen and Ringrose, 2019). Although removing the mutual exclusiveness is unrealistic for modeling H3K4me3/H3K27me3 bivalent chromatin, it could be feasible for other forms of bivalent chromatin.

Finally, by analyzing the transitions among chromatin states in response to external stimulus or noise, we demonstrated that bivalent chromatin can serve as a poised state, which has the potential to switch to multiple states and facilitates a stepwise transition between repressive chromatin state and activating chromatin state, to mediate phenotypical plasticity. Our results demonstrate why bivalent chromatin is often involved in many developmental processes where cell fate decision is required (Bernstein et al., 2006; Matsumura et al., 2015; Rugg-Gunn et al., 2010; Sachs et al., 2013). Such plasticity of bivalent chromatin also agrees with an experimental study wherein the bivalent gene Zeb1 can readily respond to microenvironmental signals and thus enables the transition to cancer stem cells (CSCs) from non-CSCs (Chaffer et al., 2013). However, the bivalent state is unlikely to be the only form of the poised chromatin. Several studies have suggested that bistability is the nature of the poised chromatin (Hong et al., 2011; Lorzadeh et al., 2016; Sneppen and Ringrose, 2019). This is not contradictory to our results. In fact, both the bivalent state and bistability reflect a balance between TrxG proteins and PcG proteins; from a mathematical point of view, under the condition with symmetric parameters, both the bivalent state and bistability can emerge, largely depending on the model and parameter setting. Taken together, the combination of the bivalent state and bistability could provide a full picture of the nature of the poised chromatin.

Limitations of the study

In the present model, we do not take into account the transcription processes, which have been implicated in the formation of bivalent chromatin (Voigt et al., 2013). Recent theoretical analyses have shown that incorporating the feedback from transcriptional products into the histone modification regulation can induce intermediate states (Huang and Lei, 2019; Zhang et al., 2019). So, it will be interesting to investigate the effect of transcription processes on the formation of bivalent chromatin in a model coupling both transcription processes and histone modification regulation.

Besides, we do not consider different timescales for H3K4me3 and H3K27me3. In our model, we assume the same demethylation timescales for both H3K4 and H3K27 (i.e., ). Under this assumption, we identify that comparable methylation rate constants for H3K4 and H3K27 (i.e., ) are required for the emergence of the bivalent state. However, demethylation reactions for H3K4 and H3K27 are carried out by different enzymes (e.g., KDM5A for H3K4me3 and UTX for H3K27me3) (Lee et al., 2007; Pasini et al., 2008) and thus may exhibit different timescales; previous experimental studies have also indicated a difference in the dynamics of H3K4me3 and H3K27me3, suggesting slow dynamics for restoration of H3K27me3 after replication and fast dynamics for H3K4me3 (Petruk et al., 2017; Reverón-Gómez et al., 2018). Although changing timescales usually does not affect the steady-state solutions of a deterministic system, it may reshape the probability landscape of a stochastic system. In this case, our model could be readily used to explore how timescale difference affects the emergence of the bivalent state. On the other hand, according to our current results, similar timescales for H3K4me3 and H3K27me3 can contribute to the emergence of the bivalent state, and thus might be the case at bivalent loci. We expect that some new technology could be developed and experimentally test the predictions of our model. Although our model only focuses on H3K4me3/H3K27me3 bivalent chromatin, we believe that the underlying mechanisms are general and could be naturally applied to other combinations of histone modifications.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Lei Zhang (zhangl@math.pku.edu.cn).

Materials availability

This study did not generate new unique materials.

Data and code availability

Simulations in this paper are performed in software MATLAB (https://www.mathworks.com/). The code is accessible here: http://bicmr.pku.edu.cn/∼zhanglei/publication.htm and https://github.com/Wei-BioMath/bivalency.

Method details

Stochastic model

We consider a linear array of nucleosomes and each nucleosome has two copies of H3 histones. There are in total N histones, indexed by 1 to N. We also assume the (2k-1)-th and 2k-th histones belong to the same nucleosome, where k=1, …, N/2. Let denotes the modification state of the i-th histone, and . In practice, we use seven numbers -3, -2, -1, 0, 1, 2, 3 to represent respectively.

The reaction rate, also termed as the propensity function, has been listed in Table 1. Despite 12 possible reactions in total, each histone i with state can take part in at most two reactions, i.e., and . The propensity function for is formulated as

and the one for iswhere and represent fractions of K4me3 and K27me3 among the neighbors of the ith histone respectively; is the number of histone neighbors; the function is the Kronecker delta. Functions f and g are chosen as Hill functions and when modeling the nonlinear feedback, and chosen as when modeling the linear feedback. The range of histone neighbors is defined according to the nucleosome interaction regime. For the global interaction regime, the neighbors are all the other N-1 histones, and . For the local interaction regime, except four histones (number 1, 2 and N-1 and N) on the boundaries, there are 5 neighbors for each histone, i.e.,

We set for all histones for the local interaction regime.

We use the Gillespie algorithm (Gillespie, 1977) to simulate the stochastic dynamics of the whole system. There are N species (i.e., histones) and 2N reaction channels in total. Denote the state of each histone at time t, and denote vector the state of the whole system. We denote and represent the two propensity functions for the i-th histone: one for and one for respectively. The state change of the whole system induced by the j-th reaction channel can be represent by a vector with the same length as (t): if j is odd, then the -th element of is 1 and the other N-1 elements are zeros; if j is even, then the -th element of is −1, and the other N-1 elements are zeros. The Gillespie algorithm is as follows:

Initialization. Set the time variable t=0, and assign desired values to (0).

Calculate and store the values of 2N propensity functions at current time t. Also calculate the sum of the 2N values, denoted as .

Generate two independent random numbers and , both of which obeys a uniform distribution over the interval [0,1], and calculate the waiting time for the next reaction to occur. Determine the integer μ according to . Here, μ is the index of reaction channel chosen to occur.

Update the time variable as , and update the system state as .

Repeat the step 2)-4) until the terminating condition is satisfied.

Besides, DNA replication during the cell cycle is also simulated in our model. Denote L the cell cycle length. If the current time plus the waiting time for the next reaction exceeded the time when DNA replication should have occurred, i.e., for some positive interger n, then we stop the Gillespie algorithm and execute DNA replication. In particular, time variable t is updated as nL rather than , and each nucleosome will be replaced by an unmodified nucleosome (i.e., both and are set to be zeros, k=1, …, N/2) with a probability of 0.5. After such DNA replication, resume the Gillespie algorithm by repeating the step 2)-4). For all stochastic simulations, DNA replication is applied unless otherwise specified.

Deterministic model

Under the assumption of the homogeneous system and global interaction regime, the mean behavior of the stochastic model can be approximated by a set of ODEs. Denote , , , , , and the number of K27me3, K27me2, K27me1, Empty, H3K4me1, K4me2 and K4me3, respectively. Then the deterministic model can be described by ODEs as follows:where

Similar ODEs can be constructed for the model variants shown in Figure 4.

Calculation of time-averaged probability and mean dwell time

Time-averaged probability. Let denotes a simulated stochastic trajectory obtained by the Gillespie algorithm, where and . At each time point , let and denote the number of elements in vector that are equal to “K4me3” and “K27me3”, respectively. As seen in Figure 1D, we define four regions , , and to represent territories of the activating state, repressive state, bivalent state and unmodified state, respectively, with the following forms:

An indicator function is used to classify the four chromatin states at each time point:

Finally, the time-averaged probability for each chromatin state Z is calculated based on the formula:

Mean dwell time. For each chromatin state Z, we first extract the start time points when the system enters this state and the end points when the system leaves. In particular, is defined as a subsequence of , which satisfies and for ; likewise, is a subsequence of , which satisfies and for . Then the mean dwell time (MDT) for each chromatin state Z can be calculated bywhere is the number of elements in . In practice, to avoid extreme short dwell time due to noisy transitions at state boundaries, we apply moving average (see MATLAB function “smooth”) to smooth the trajectory to obtain a robust estimation of the MDT.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Data used in Figure 3G	(Liu et al., 2016)	https://doi.org/10.1038/nature19362
Software and algorithms
MATLAB code to perform stochastic and deterministic simulations	This paper	http://bicmr.pku.edu.cn/∼zhanglei/publication.htm & https://github.com/Wei-BioMath/bivalency
The Geometric Minimum Action Method	(Heymann and Vanden-Eijnden, 2008)	https://doi.org/10.1002/cpa.20238
MATLAB	MathWorks	https://www.mathworks.com/