Bottom-up modeling of chromatin segregation due to epigenetic modifications.

We use a chromosome-scale simulation to show that the preferential binding of heterochromatin protein 1 (HP1) to regions high in histone methylation (specifically H3K9me3) results in phase segregation and reproduces features of the observed Hi-C contact map. Specifically, we perform Monte Carlo simulations with one computational bead per nucleosome and an H3K9me3 pattern based on published ChIP-seq signals. We implement a binding model in which HP1 preferentially binds to trimethylated histone tails and then oligomerizes to bridge together nucleosomes. We observe a phase reminiscent of heterochromatin-dense and high in H3K9me3-and another reminiscent of euchromatin-less dense and lacking H3K9me3. This segregation results in a plaid contact probability map that matches the general shape and position of published Hi-C data. Analysis suggests that a roughly 20-kb segment of H3K9me3 enrichment is required to drive segregation into the heterochromatic phase.

Chemical modifications to a chromosome correlate with its 3D organization in the nucleus (1), and this 3D organization is integral to controlling gene expression (2). Thus, identifying the physical mechanisms by which epigenetic modifications control chromosomal organization would provide a predictive framework for bridging epigenetics and expression.

The fundamental organizational unit of chromosomal DNA is the nucleosome, which consists of a 147-bp segment of DNA wrapped around eight histone proteins. Among the host of epigenetic modifications, methylation of the ninth lysine of the C-terminal tail of histone protein 3 (i.e., H3K9) is a prototypical example where enrichment of such a modification leads to large-scale changes in physical state and biological function. Enrichment in H3K9 methylation correlates with the DNA and associated proteins composing chromatin existing in a heterochromatic state that is densely packed and generally less expressed than more loosely packed euchromatin (2). Roughly micrometer-size domains of euchromatin and heterochromatin are observed in cell nuclei by electron (3) and fluorescence microscopy (4). Predicting how epigenetic markers affect the physical segregation of the chromosome is critical to interpreting the biological consequences of the epigenetic code (5).

The determination of the 3D organization of DNA from nucleosomes to the full chromosome remains a challenge. Interphase chromatin organization appears quite random (6), with cell-to-cell variability in relative loci positions (7). The most prominent nonrandom feature of chromatin organization is the physical connectivity of the DNA backbone itself, and physical models of chromatin must begin with this connectivity (8–21).

Experimental advances in chromosome-capture techniques (e.g., Hi-C) have revealed many additional features of chromatin organization (1, 22). These techniques use chemical cross-linking of loci that are in spatial proximity, followed by sequencing. The most prominent feature of the resulting contact maps is a plaid-like pattern of high contact probability indicative of two or more types of chromatin that self-associate.

The observed Hi-C patterns have inspired a number of block copolymer models of chromatin where the blocks represent genomic regions with various epigenetic marks (17–21). Simulations of alternating block copolymers produce plaid contact patterns reminiscent of the Hi-C data, but do not incorporate genomic data (17, 19, 23). Jost et al. (19) demonstrate that diagonal Hi-C features are predicted when the blocks of the polymer correspond to epigenetic domains and that a few off-diagonal features are retained when the simulations are seeded with initial configurations containing them.

Further efforts aim to reproduce the Hi-C data with in silico block copolymers, using a top–down approach where model parameters are trained for agreement with Hi-C contact maps (17, 20, 21). DiPierro et al. (20, 21) classify the chromatin into six block types, using 11 epigenetic marks. After training the interaction constants between these block types as well as nonspecific binding energy as a function of genomic distance between loci, the model is able to predict Hi-C contact patterns, conclusively showing that the information in the epigenetic marks is sufficient to predict much of the observed Hi-C pattern. However, the DiPierro et al. (20, 21) model does not explain the physical mechanisms behind the fitted interactions, and the coarse-grained nature of the model makes it difficult to predict biological mechanisms for how chromosomal segregation affects accessibility and how epigenetic marks are established and maintained.

In this paper, we use a bottom–up approach to modeling chromosomal DNA. We develop a physical polymer model of chromatin and explore the consequences of introducing a binding model where heterochromatin protein 1 (HP1) is more likely to bind to trimethylated nucleosomes, which we identify from H3K9me3-ChIP. In contrast to the top–down approach of fitting Hi-C, we use only experimentally measured or physically motivated parameters in our model and observe to what extent this untrained model reproduces the Hi-C data. Our untrained model is not susceptible to reproducing experimental biases in the Hi-C data and is advantageous for establishing the causal effect of the introduced interactions on the organization. While we do not completely reproduce the Hi-C contact maps, we observe that this binding model alone explains much of the plaid pattern. Visualizing our modeled microstructure, we observe a phase segregation reminiscent of heterochromatin and euchromatin. Additionally, the model aids in the interpretation of the ChIP-seq data by providing physically motivated guidelines for predicting whether a chromosomal segment is segregated.

The Model

Chromosomal DNA.

We simulate each of the nucleosomes of human chromosome 16 as an individual bead. We use the shearable, stretchable wormlike chain model described in (24–26), which captures the energetics of a semiflexible polymer chain while allowing for the greater bead spacing needed for large-scale simulations. For simplicity, we model chromatin as having the same stiffness as bare DNA with nucleosomes spaced a backbone path length nm apart, corresponding to a 50-bp linker between nucleosomes. See for the effect of reducing this stiffness. A more detailed model of chromatin would include the effects of geometrical constraints from bound proteins, chain topology, linker histones, the entry and exit angles of the DNA from each nucleosome, geometric effects that tend to align adjacent fibers, and variation in linker lengths. Here, we present a bottom–up model on which detailed effects can be built.

Confinement.

The DNA is highly confined by the nucleus and, in the absence of a confining boundary, it would expand to occupy an unconfined radius of gyration much larger than the radius of the nucleus. We simulate only a single chromosome (of 46), so we choose a confinement m in diameter, roughly the size of a chromosome territory (27). We argue that since each chromosome is primarily found within its chromosomal territory in the nucleus (28), this smaller single-chromosome confinement is physically reasonable.

Nonspecific Interaction.

Interactions between chromosomal segments include electrostatic interactions, van der Waals forces, steric interactions, and interactions with a complex solvent containing a zoo of interacting proteins. We take a coarse-grained approach and define the interaction free energy from the sum of these effects to bewhere is the volume fraction of chromatin calculated inside discrete bins of width . To calculate the total interaction free energy, we divide the system into cells of volume and use linear density interpolation to calculate for each cell as in ref. 29. We choose to capture density fluctuations on the size scales of interest. Each bead is given a volume of , roughly corresponding to the volume of a nucleosome. We choose the parameter to distribute chromatin without leaving large chromatin-less regions, but not so large to prevent significant variation in chromatin density. See for details on the choice of . While this binned density approach does not guarantee individual nucleosomes do not overlap, it is appropriate for capturing equilibrium density fluctuations of disordered chromatin on length scales greater than .

Methylation Profile.

We focus on how chromatin organization depends on H3K9me3, which is associated with pericentric heterochromatin (30). We assume a fixed H3K9me3 profile based on ChIP-seq data from GM12878 cells (31, 32). We divide the ChIP-seq signal into 200-bp bins, corresponding to our nucleosome spacing. As there are two H3 tails per nucleosome, we apply cutoffs to the binned signal to classify each nucleosome into both tails trimethylated, one tail trimethylated, or neither tail trimethylated. For simplicity, we used evenly spaced cutoffs chosen so that roughly half of all tails are classified as trimethylated. See for details. Our trimethylated state is meant to capture the effects of in vivo trimethylation [20–30% of histone tails (33, 34)] on HP1 binding, and to a lesser extent, the effects of dimethylation [30–40% of tails (33, 34)].

H3K9me3-Dependent HP1 Binding.

The chromatin regulatory protein HP1 binds histone tails through its chromodomain. Its chromo-shadow domain oligomerizes with other HP1s from different histone tails, which may be from different nucleosomes in close spatial proximity (35). HP1 preferentially binds to trimethylated tails, condensing regions of the chromatin marked with H3K9me3. In this work, we do not have the detail to differentiate the paralogs HP1- and HP1-, so we refer only to the combined effect of the two.

We follow Mulligan et al. (30), using data from Canzio et al. (36) to find methylated or unmethylated histone-HP1 binding affinities and . When nucleosomes with bound HP1 come within an interaction radius, they experience a free-energy benefit from oligomerization (30). In this work, we choose the interaction distance as based on an estimate of the histone tail length (37). We use a two-state model (bound by an HP1 or not bound) for each H3 tail, biased by the chemical potential induced by the concentration of free HP1. For computational purposes, we use a coarse-grained approach for calculating the interaction between separate nucleosomes based on , the local number density of bound HP1 interpolated from a linear grid with spacing as in ref. 29. In summary, we include binding energies (i) between HP1 and histone tails, (ii) for the oligomerization of HP1, and (iii) for explicit intranucleosomal HP1 interaction,where , is the number of nucleosomes, is 1 if the jth H3 tail of the ith nucleosomes is bound by HP1 and 0 otherwise, and is or , depending on the methylation state of the corresponding histone tail. Note that to prevent double counting of the interaction between two tails of the same nucleosome.

Monte Carlo Algorithm.

We use Metropolis Monte Carlo sampling to draw configurations from the equilibrium ensemble of the above system, which has total energy from polymer-chain deformation energy, nonspecific repulsion, and HP1 binding: . Our simulations use moves which rotate a segment of polymer about the axis which runs through its ends, change the binding state of the histone tails, rotate a single bead, translate beads, and pivot the end of the chain. Moves were performed repeatedly to bring the system to equilibrium (). Our model does not include any topological constraints on the DNA. Therefore, our individual configurations are representative of thermodynamic equilibrium in the case where topoisomerase activity is sufficiently high to relax topological constraints at equilibrium. Our Fortran source code was developed based on previously published code (38); the version used in this paper corresponds to the branch MBS2018_PNAS of the GitHub repository (https://github.com/SpakowitzLab/wlcsim/tree/MBS2018_PNAS).

Results

H3K9me3–HP1 Binding Leads to Phase Segregation.

Fig. 1 shows a 115-nm–thick slice of a snapshot of our simulation. Phase segregation is observed between a dense phase enriched in H3K9me3 (purple), reminiscent of heterochromatin, and a more loosely packed region enriched in nucleosomes with neither tail trimethylated (cyan). A similar phase segregation driven by HP1 has been observed in vivo (39).

Fig. 1.

(Top) Schematic representation of our model of chromosomal DNA. Each bead represents a single nucleosome and is colored cyan, tan, or purple if neither, one, or both of its histone-3 tails are trimethylated at lysine 9, respectively. Small black beads represent the number of HP1 proteins bound to a nucleosome. The chromatin density and the number density of bound HP1 are calculated based on the displayed grid. (Bottom) A 115-nm–thick slice through the center of a simulation snapshot with an HP1 concentration of M.

The nucleosome density difference between the two phases is controlled by the nonspecific interaction function . The choice of shown in Fig. 1 is sufficiently small that the density of the heterochromatic phase is limited by the volume fraction cutoff of . See for determination of .

Heterochromatin is typically found in the periphery of the nucleus, although examples exist of inverted nuclei such as the rod photoreceptors in nocturnal mammals (40). Fig. 1 exhibits a structure with heterochromatin in the interior of the chromatin territory because no interaction between the chromatin and the confinement is included other than keeping the beads within the confinement. An interaction between the chromatin and nuclear lamina could be added to cause heterochromatic regions to move to the periphery, such as that used by ref. 41. In this paper, we focus on the thermodynamics of the binding model with as few free parameters as possible. We note that despite their striking visual difference, regular and inverted nuclei have been shown to have quite similar Hi-C profiles (41).

HP1 Concentration Dependence.

We probe the effects of HP1 on segregation by varying its total concentration. The fraction of histone-3 tails bound by HP1 is plotted in Fig. 2 for nucleosomes with neither, one, or both tails methylated. At low HP1 concentration the probability that two bound HP1 molecules from different nucleosomes encounter each other is small. Without the oligomerization of HP1, the chromatin is a free polymer which distributes itself randomly throughout the nucleus as shown in Fig. 2. With increasing HP1 concentration the average fraction of tails bound by HP1 increases first for nucleosomes with fully trimethylated tails and later for nucleosomes with only one tail trimethylated. At intermediate concentrations (e.g., 285 M, as displayed in Figs. 1, 2, and 3) there are enough bound HP1s in regions high in H3K9me3 to cause them to condense but not enough in regions low in the mark. This leads to the segregation of the regions. At sufficiently high HP1 concentrations nearly all nucleosomes are bound and condensed.

Fig. 2.

(A–D) Simulation snapshots at progressively higher HP1 concentrations. The plot shows, for each type of nucleosome, the fraction of tails bound by HP1 at different total HP1 concentrations. Cyan corresponds to nucleosomes with neither tail trimethylated, tan to one tail trimethylated, and purple to both tails trimethylated. The larger points correspond to the snapshots displayed in A–D with C being at M, the value displayed in Figs. 1 and 3.

Fig. 3.

(Top Left) Simulation slice with colors corresponding to the number of tails trimethylated. Pie charts show proportions of each color in the interior of each phase. (Top Right) Same as Top Left except the beads are recolored by the average methylation of their genomic neighbors, using a sliding window of 101 beads. Cyan corresponds to an average fraction trimethylated 0.53, followed by tan, and finally purple for average trimethylation 0.705. (Bottom) Accuracy with which the sliding-window average correctly classifies bead regions of high and low density as measured by number of neighboring beads. Optimal window size is 100 nucleosomes.

We constructed a specific binding model with experimentally derived parameters. However, any model with an attractive potential (in our model HP1 coupling ) along with a specificity for heterochromatin (in our model differential binding energy ) would show similar behavior. For example, we model the variability in ChIP-seq signal from H3K9me3 as coming from differing methylation of histones. At least some of this variability is caused by nucleosome density, which has a similar specificity in that it would also increase the bound HP1 content and H3K9me3 signal. We present HP1 concentration as modulating HP1 binding and, hence, the degree of compaction. In this regard, the HP1 level in the cell acts as a global regulator of chromatic state. Alternatively, phosphorylation of serine 10, by reducing the binding affinity of HP1 to H3K9 (42), could regulate chromatin architecture in the same way as we presented for HP1 concentration.

Predicting Chromatin State from H3K9me3.

Our simulations provide a means of predicting whether each nucleosome is segregated into the densely packed heterochromatic phase or the more lightly packed euchromatic phase based on the H3K9me3 profile from ChIP-seq. We aim to extract general rules from our simulations for interpreting arbitrary ChIP-seq profiles for different chromosomes, cell types, or species without the need to perform simulations for each case. While it is true that the methylation state of a nucleosome is correlated with its chromatin state, the snapshot in Fig. 3, Top Left indicates that there are a number of unmethylated nucleosomes (cyan) in the heterochromatic region (dense, mostly purple phase) and vice versa. The phase in which a nucleosomes resides is not determined just by the methylation state of that nucleosome, but also by the methylation state of its neighbors on the chromatin fiber. As in other copolymer systems (43), the strength of segregation is determined both by the strength of the interactions (defined by in our model) and by the length of the segments with common chemical identity (the size of methylated regions in our model). Thus, when predicting in which phase a particular nucleosome is found, one must assess the average methylation of its genomic region.

It is not a priori clear how many adjacent nucleosomes over which to average the degree of methylation when predicting the chromatin state of a nucleosome. Longer methylated regions are more likely to be segregated into the heterochromatic phase due to the additive binding energy of more HP1 molecules. The chain connectivity ensures that proximal genomic segments are in close spatial proximity, leading to more prevalent contacts and higher cooperativity in HP1 binding. However, genomic loci that are sufficiently separated do not cooperatively affect each other’s chromatin state, since the chain connectivity does not ensure that they are in close spatial proximity. Therefore, the thermodynamic advantage of cooperative HP1 binding starts to break down once the region size becomes larger than the characteristic length of chromatin in the heterochromatic phase.

We define a simulated nucleosome as being in the heterochromatic phase if it is surrounded by a high density of chromatin (). We then predict which phase each nucleosome will be in, based on the average methylation state nucleosomes to either side. Fig. 3 displays the accuracy of this prediction for a range of . Classification based on the optimal width, , is displayed in Fig. 3, Top Right. For simplicity, we present a sliding-window average; similar results are obtained by replacing the window with a Gaussian or power-law filter ().

The above observations provide guidance for classifying the genome into regions of heterochromatin and euchromatin based on ChIP-seq data for H3K9me2/3. For this classification, our model recommends an average of the H3K9me3 signal over a window kb wide. Regions with low averaged signal will be euchromatic, regions with high ChIP-seq signal will be heterochromatic, and intermediate levels will be at the boundary or stochastically incorporated into either phase.

Beyond its usefulness for interpreting ChIP-seq data, the optimal window size implies that modulating the methylation level of a roughly 20-kb–long strand of chromatin would cause it to move between heterochromatic and euchromatic phases in the nucleus. To emphasize this, we perform simulations with a modified input methylation sequence containing blocks of fully methylated (purple in Fig. 4) or fully unmethylated (cyan) nucleosomes interspersed among the ChIP-seq–based methylation profile. As shown in Fig. 4, the blocks that contain only five (un)methylated nucleosomes are scattered throughout both phases. Their positions are determined by their surrounding methylation profile. Progressively larger blocks are more likely to be incorporated into their respective phases, with blocks of at least 80 nucleosomes being well segregated. We note that this analysis differs from that presented in Fig. 3 because the latter incorporates the autocorrelation in methylation state present in the original ChIP-seq. The decaying number of neighbors to either side of the inserted block in Fig. 4 suggests that segregating blocks drag roughly 10-kb regions with them into their respective phases. The decay at the edge of unmethylated (cyan) blocks is generally outside the block while the decay around methylated blocks (purple) starts inside the block. This indicates that loci near epigenetic boundaries tend to sway euchromatic as measured by local density.

Fig. 4.

Simulations with blocks of 5, 10, 40, and 80 nucleosomes in the input methylation sequence set either fully methylated (purple) or fully unmethylated (cyan). The two black curves show the average fraction of histone tails trimethylated in the vicinity of the inserted blocks. The average number of neighbors—a measure of chromatin compaction—shows that 80 nucleosome blocks segregate into their respective phases (as seen in the snapshots) and drag roughly 10-kb adjacent regions with them. In the snapshots the blocks are shown in color and surrounding chromatin with wild-type methylation is shown in gray.

Phase Segregation Reflected in Contact Maps.

In Fig. 5, we show a comparison of the simulated contact map (Lower Left triangle) with a Hi-C contact map for GM12878 cells (1) (Upper Right triangle). Bright red corresponds to a high likelihood that the two loci are in contact while white corresponds to a low probability. The simulated contact map was calculated by counting beads within as “in contact” and averaging over 18 snapshots. Fig. 5 was generated at an intermediate HP1 concentration of M (see for Hi-C maps at other HP1 concentrations). We emphasize that none of the simulation parameters (apart from the color scale) are adjusted to fit the Hi-C data. Our goal is to understand which features of the Hi-C data are predicted by the model as a way of both verifying the model and interpreting the Hi-C data.

Fig. 5.

Comparison of simulated contact probabilities (Lower Left triangle) and those from Hi-C data from Rao et al. (1), GEO accession no. GSM 733664 (Upper Right triangle). The numbering is in megabase pairs along chromosome 16. The region shown is the majority of the arms of the chromosome. The model predicts multiple features of the Hi-C contact map.

A prominent feature of the Hi-C map is the plaid pattern with rectangles of 1–4 Mb. Many of the rectangles found in the simulated contact map agree in both size and genomic location with the experimental data. In the simulation, this pattern is caused by regions high in H3K9me3 condensing together into a dense phase with elevated contact probability. While the corresponding rectangles in the Hi-C data could be caused by something else that is correlated only with H3K9me3, the affinity between HP1 and H3K9me3 provides a plausible mechanism. There are also rectangles in the Hi-C data which are not matched by the simulation. Some of these may be the result of interactions with epigenetic marks other than H3K9me3 that are not included in the model. Biochemical measurement of the interactions selective to other marks would allow for a more complete physics-based model of chromatin. We note that while the Hi-C pattern appears like a checkerboard (i.e., corners touching), the simulation pattern has squares surrounded by continuous white bands. This difference may suggest that the heterochromatic and euchromatic regions are closer in density than in the model. Alternatively, the Hi-C data collection method may not result in a higher read count from a dense area in the way that counting neighbors does.

The Hi-C data also contain a number of white lines with widths on the order of 100 kb. These features are reliably matched by the simulation. However, this good agreement could also be explained by the existence of genomic regions that Hi-C and ChIP-seq are biased against and therefore show up as white in both.

Another prominent feature of the Hi-C contact map is the bright region near the diagonal, representing a higher contact probability for regions which are close genomically. While the connectivity of the chromatin polymer in our model does result in an increased contact probability for regions very close to each other genomically, this effect influences the contact probability only at short length scales not visible in Fig. 5. On the scale of megabases, the simulation shows little preferential contact for genomically proximal loci. This difference in contact probability scaling has been explained in the past by arguing that the chromatin is trapped in a nonequilibrium fractal globule state (15) or by the action of loop extrusion proteins (13). Our simulations confirm that extra physics are required to explain this scaling.

Summary

The 3D chromatin architecture and its regulation by epigenetics are a complex system involving many interactions, most of which have unknown energetics. By introducing a simple model with as few unknown parameters as possible, we can understand the implications of these few interactions. Our model consists of a self-repulsive polymer, discretized on the single-nucleosome level and placed in a confinement commensurate to a single-chromosome territory. The only distinguishing feature we give the beads is a methylation state based on the ChIP-seq signal for H3K9me3. Introducing an HP1-binding model based on measured energetics (30, 36), we show that these interactions are capable of causing a phase segregation into a dense, H3K9me3-rich phase and a less dense, H3K9me3-lean phase. This segregation is capable of reproducing some of the features of the Hi-C map on the megabase size scale. We leverage this model to find that translocation of a nucleosome into heterochromatin is governed by the methylation state of the 20 kb surrounding it. Our approach is amenable to further exploration of the physical phenomena that underlie the organization and function of chromatin in a living cell.