A simulation model of heterochromatin formation at submolecular detail.

Heterochromatin is a physical state of the chromatin fiber that maintains gene repression during cell development. Although evidence exists on molecular mechanisms involved in heterochromatin formation, a detailed structural mechanism of heterochromatin formation needs a better understanding. We made use of a simple Monte Carlo simulation model with explicit representation of key molecular events to observe molecular self-organization leading to heterochromatin formation. Our simulations provide a structural interpretation of several important traits of the heterochromatinization process. In particular, this study provides a depiction of how small amounts of HP1 are able to induce a highly condensed chromatin state through HP1 dimerization and bridging of sequence-remote nucleosomes. It also elucidates structural roots of a yet poorly understood phenomenon of a nondeterministic nature of heterochromatin formation and subsequent gene repression. Experimental chromatin in vivo assay provides an unbiased estimate of time scale of repressive response to a heterochromatin-triggering event.

Introduction

Heterochromatin gene repression plays a central role in cell development and differentiation through intricately timed targeted gene repression (Allshire and Madhani, 2017; Grewal and Moazed, 2003). In particular, heterochromatin formed on regulatory regions of transcription factors (TF) Oct4, Sox2, and Klf4 locks newly differentiated cells in their specialized states. Conversely, in regenerative medicine applications derepression of Oct4 transcription is a key trigger to cell reprogramming into induced pluripotent stem cells (iPSC) (De Los Angeles et al., 2015; Jerabek et al., 2014). Because of the importance of heterochromatin for normal development and cell state maintenance, disruption of this pathway represents a key role in multiple human diseases including hard to treat cancers (Carone and Lawrence, 2013; Dialynas et al., 2008; Morgan and Shilatifard, 2015).

Multiple facets of heterochromatin gene repression have been revealed. In the late 1930s, heterochromatin, which was then seen as densely stained stripes in chromosome bodies, was linked to the regulation of gene activity (Brown, 1966). In the early 2000s, heterochromatin was comprehensively characterized by molecular biologists as a chromatin region enriched in dimethylation and trimethylation of histone H3 at lysine 9 (H3K9me2/3), DNA methylation marks, and the binding of heterochromatin protein 1 (HP1) (Fischer et al., 2009). Most recently, several groups have hypothesized that phase-separated droplets displaying reduced diffusion could play a major role in driving heterochromatin assembly (Larson and Narlikar, 2018; Tatarakis et al., 2017) and that HP1 has an intrinsic propensity to form phase-separated droplets (Sanulli and Narlikar, 2021). Molecular mechanisms of this liquid-liquid or polymer-polymer phase separation are yet to be understood (Erdel and Rippe, 2018). On the one hand, formation of ordered and collapsed chromatin globules may be driven by HP1-mediated chromatin bridging (Hiragami-Hamada et al., 2016; Kilic et al., 2018; Machida et al., 2018). On the other, heterochromatin formation and related phase separation may arise from the inherent propensity of HP1 to form phase-separated condensates (Larson et al., 2017; Strom et al., 2017). Despite this overwhelming experimental evidence, the structural mechanisms underlying the heterochromatin formation process are yet to be fully understood. It was long believed that heterochromatin exists, at least partly, in the form of a 30 nm fiber (Felsenfeld and McGhee, 1986), a periodic arrangement of nucleosomes. However, recent electron microscopy (EM) studies in live cells did not detect significant presence of 30 nm fiber structures (Nishino et al., 2012). Unfortunately, the kilobase-scale 3D structure of heterochromatin at molecular resolution is currently in a blind spot of experimental structure determination techniques. It is too small for large-scale scale techniques, such as, ChromEMT (Ou et al., 2017), ChIA-PET (Tang et al., 2015), Hi-C (Stevens et al., 2017), or 3D-FISH (Markaki et al., 2012), but too large, flexible, and heterogeneous for atom-scale or molecule-scale NMR, cryo-EM, or x-ray crystallography.

There is an ongoing effort to computationally and theoretically comprehend molecular structure, function, and properties of the chromatin (D. Bascom and Schlick, 2018; Korolev et al., 2016; Laghmach et al., 2021, 2020; Moller and de Pablo, 2020; Ozer et al., 2015). A number of simulation studies demonstrated that high concentrations of multivalent DNA-binding particles, representing generic transcription factors or epigenetic effectors, drive chromatin condensation, and phase separation (Barbieri et al., 2012; Brackley et al., 2016; Jost et al., 2014; MacPherson et al., 2018; Michieletto et al., 2017a; Nuebler et al., 2018). Several other studies focused on the "glassiness" of the condensed chromatin fiber (Michieletto et al., 2017b; Shi et al., 2018). More generally, ad hoc particle-based simulations proved successful in modeling various aspects of the chromatin structural transformations (Buckle et al., 2018; Falk et al., 2019; Gürsoy and Liang, 2016; Laghmach et al., 2020, 2021; MacPherson et al., 2018; Vasquez et al., 2016; Verdaasdonk and Bloom, 2011). Worm-like chain is a typical model used to represent the chromatin fiber (Kang et al., 2015; MacPherson et al., 2018; Michieletto et al., 2016, 2017a; Shi et al., 2018). It enables a straightforward implementation using standard algorithms and force fields. Most of the reported models (Barbieri et al., 2012; MacPherson et al., 2018; Michieletto et al., 2016, 2017a; Shi et al., 2018) were tuned to correctly reproduce the chromatin behavior at multi-megabase scale. Most models of protein-mediated heterochromatin formation represented chromatin-binding proteins implicitly as a nucleosome state modifier, turning such modified nucleosomes into multivalent attractors of other modified nucleosomes (Barbieri et al., 2012; MacPherson et al., 2018; Michieletto et al., 2016, 2017a; Shi et al., 2018). To enable chromosome-scale simulations, the fiber is often composed of particles embedding multiple nucleosomes (Michieletto et al., 2016; Nuebler et al., 2018; Shi et al., 2018) with persistence length longer than nucleosome spacing (Brackley et al., 2013; MacPherson et al., 2018).

We developed a model to simulate the HP1-mediated chromatin compaction to a greater detail than previous models to incorporate previously unaccounted factors affecting the dynamics of the process. To more accurately reflect the entropic burden of the HP1-mediated compaction, the model explicitly represents molecular and submolecular entities involved in the process. For instance, HP1 is represented by chromo-domains and chromo-shadow-domains of HP1 separated by a disordered linker or disordered histone tails. The model parameters were derived from experimental data for protein-protein binding and particle diffusion in the cellular environment, enabling higher spatial and temporal resolution. Although our model displays a comparable to previously reported models’ degree of coarse-graining (Barbieri et al., 2012; MacPherson et al., 2018; Michieletto et al., 2016, 2017a; Shi et al., 2018), we introduced four independent reversible events – binding of two nucleosomes by HP1 chromodomains, a close encounter of HP1-bound nucleosomes, and a CSD-mediated dimerization of the nucleosome-bound HP1 proteins – to better reflect the entropic cost of the fiber condensation. Furthermore, compared to previous work, our fiber model allows a higher degree of compaction consistent with the evidence of kilobase-scale heterochromatinization on gene promoter regions (Hathaway et al., 2012) and compact nucleosome clusters visible in single-cell ChromEMT images (Ou et al., 2017). On the whole, this simple model provides a surprisingly complete interpretation of what we know about heterochromatin.

The new model was applied to better understand the structural mechanism of heterochromatin formation on short regulatory regions, such as an approximately 10kb-long Oct4 promoter. In particular, we investigated whether it would allow a bridge-mediated mechanism of fiber condensation. It is still an open question whether the formation of compact heterochromatin domains is driven by a liquid-liquid phase separation (LLPS) because of an intrinsic propensity of HP1 to form separated droplets or by HP1-mediated bridging of sequence-remote nucleosomes (Erdel and Rippe, 2018; Hiragami-Hamada et al., 2016; Larson et al., 2017; Machida et al., 2018; Strom et al., 2017). Since previous computational efforts have already explored the LLPS-driven compaction hypothesis (MacPherson et al., 2018), in this study we focus on the possibility of bridging mechanism in a system with parameters derived from experimental data on interaction and diffusion of HP1 and H3K9me3-marked nucleosomes. The bridging hypothesis was recently supported by experimental data suggesting that heterochromatin can adopt compact states without showing hallmarks of HP1-driven LLPS (Erdel et al., 2020). Another question we examined here is a nondeterministic response of a 10kb region to conditions inducing heterochromatin formation, which is recruitment of HP1 and histone methyltransferases to the locus of interest. As we show in the experimental section using the Chromatin in vivo Assay (CiA), the repressive response to a triggering event may take hours to days. Owing to the limitations of single cell experimental detection of chromatin fibers, the structural roots of such a nondeterministic response were unclear. Here, our computational simulations provide insights into both the aforementioned questions.

Model

The simulated system consists of a pre-methylated chromatin fiber and multiple copies of HP1 protein at a putative cellular concentration level. Nucleosomes, the constitutive units of chromatin, and HP1 proteins are composed of particles and feature submolecular details essential for heterochromatin formation, as described in the following lines. Particles may represent objects as large as protein complexes and as small as individual residues. The submolecular features enable the system to have targeted pairwise interactions. Separation in this manner creates the ability to mimic specific multivalent interactions of multi-domain proteins and multi-protein epigenetic proteins. Particles may be tethered to one another to form a body or a chain of particles. A body is a set of particles in which particles do not move relative to each other. A chain is a set of particles in which each particle is constrained to be within a defined distance range from its chain neighbors. Each particle is characterized by its mass (a readily accessible and experimentally measurable property) (see Table 1 for the particles’ masses and radii). We assume that the particle/body volume and radius can be calculated from its mass through a collective estimate of the protein density (Fischer et al., 2004). The model was simulated by a Monte Carlo method where the displacement of a particle on each Cartesian coordinate was drawn from a respective zero-centered Gaussian distribution G(Δx), G(Δy) or G(Δz) with a variance s derived from Equations 1, 2, and 3. From Fick’s 2nd law, the root mean square displacement of a freely diffusing Brownian particle of type A can be represented aswhere, Δt is the time step and D is the diffusion coefficient, which, according to Stokes-Einstein relationship, can be written aswhere k is the Boltzmann constant, T is temperature, η is viscosity, and r is the radius of the particle of type A, which, for a spherical particle, can be written aswhere m is mass and d is density. Density was assumed to be the average protein density estimated as 1.35 g/cm3 (Fischer et al., 2004). The viscosity of the cell lysate was previously estimated to ∼3 cP (Arosio et al., 2016; Ye et al., 2013).

Table 1

Individual particles’ parameters utilized in the simulations

Particle	Radius (nm)	Mass (kDa)
Nucleosome_core	2.5	200
H3_core	1	20
H3_K9	1	10
H3_K9me3	1	10
HP1_cd	1	15
HP1_csd	1	15

kilodaltons.

The generated new position of a particle must satisfy acceptance criteria intended to prevent the violation of distance constraints (see below) that tether members of a chain to each other. Whenever the new position causes a deviation from a canonical distance range, the displacement vector is damped exponentially as exp(-kΔd), where k is a spring-like constant for a given tether and Δd is the magnitude of the distance past the allowed tether length range. The tether-specific constants k were chosen to keep average tether lengths within the allowed range (see Table 2 for a full set of tether parameters) and avoid frequent overstretching (as can be seen from the sampled tether lengths distributions in Figure S5A).

Table 2

Parameters of particle-particle tethers utilized in the simulations

Particle 1	Particle 2	Dmin (nm)a	Dmax (nM)	kb
Nucleosome_core	Nucleosome_core	13	17	0.5
HP1_cd	HP1_csd	2	4	0.3
H3_core	H3_K9	1	4	0.3
H3_core	H3_K9me3	1	4	0.3

Minimal and maximal particle-particle distances allowed without a penalty applied.

as well as the “spring-like” constant (k).

The time step in Equation 1 was set to 1 μs to maximize the sampling rate while keeping the particle displacements within the system’s resolution, which is by avoiding overstretching the particle-to-particle tethers. The dependence of the displacement magnitude on time step enables, in theory, assessment of a time scale for the events occurring in our simulations. We realize though, especially in the light of recent research on phase separation (Larson et al., 2017; Strom et al., 2017), that the diffusion rates in a cellular compartment of interest may differ from bulk measurements (Chow and Skolnick, 2015; Young et al., 1980) used to parameterize the diffusion in our simulations. Hence, the time intervals given throughout this report should be considered as reference values to compare the course of events between simulation runs (though orders of magnitude are consistent with the experimental time course visualized with CiA results). The time step of 1 μs was validated by calculating diffusion parameters of the system’s molecules from the simulation data using Equations (Equation 1), (Equation 2), (Equation 3). For instance, a posterior estimate of the time step based on the mean square displacement of nucleosome particles (nucleosome-related experimental data were not used to parameterize the particles’ displacement), we found that the effective time step is 0.941 μs. Moreover, the HP1 diffusion coefficient is 0.550 μm2/s, which is in line with 0.6–0.7 μm2/s determined by Schmiedeberg et al. (2004) and close to 1.4 μm2/s by Muller et al. (Müller et al., 2009).

Chromatin fiber

The chromatin fiber was simulated as a chain of 51 nucleosomes. A nucleosome was represented as a body composed of three particles. Two of these particles are histone H3 core domains and the third represents the remaining portion of the nucleosome core (see next section for the description of H3 and its interactions and Tables 1, 2, and 3 for the particles’ parameters). Nucleosomes were allowed to freely move within a distance range between 13 and 17 nanometers (nm) from their adjacent neighbors with an exponential damping beyond these limits. The damping multiplier, applied to the distance Δd by which the allowed distance range would be violated in case of an undamped displacement was equal to exp(-kΔd) with k = 0.5 for nucleosome-nucleosome tethers. The distance constraints above were informed by a large body of experimental evidence on nucleosome spacing. Linker DNA varies broadly in length depending on species or cell types and can be as short as 20 base pairs (bp) (6 nm) and as long as 90 bp (27 nm) (Baldi et al., 2018; Singh and Mueller-Planitz, 2021). However, in most human cells, the linker lengths are closely distributed around an average of ∼55 bp (∼15 nm) (Kornberg and Lorch, 2007; Kornberg and Stryer, 1988; Schones et al., 2008). Linker lengths may evolve in time because of nucleosome sliding, either spontaneous or mediated by ATP-dependent chromatin remodeling complexes (Bowman and Deindl, 2019; Brandani et al., 2018; Lequieu et al., 2017; Mueller-Planitz et al., 2013). Furthermore, at a constant linker length, distances between adjacent nucleosomes may vary because of a number of factors including bending, twisting, and other double-strand DNA fluctuations (Strick et al., 2000), as well as possible kinks that putatively occur on the interface of DNA-protein interaction (Connolly et al., 2018; Schindler et al., 2018). Overall, the fiber configurations sampled by this Monte Carlo simulation model are consistent with the most recent experimental evidence suggesting a high degree of disorder of the chromatin fiber in live cells (Baldi et al., 2020; Fussner et al., 2011; Maeshima et al., 2019; Nozaki et al., 2017), as well as with available EM images of 10-nm chromatin fiber somewhat reminiscent of a freely jointed chain (Baldi et al., 2020; Grigoryev et al., 2009; Woodcock and Horowitz, 1997). We verified that the nucleosome-to-nucleosome tether parameters are consistent with the time step of 1 μs in Equations 1 and that overstretched or overcompressed tethers constitute only a minor population. In particular, the probability density function (PDF) for nucleosome-to-nucleosome distances (Figure S5A) shows that over 90% of nucleosome particle displacements resulted in distances within the range of 13–17 nm with a median value of 15.4 nm. We also verified that our fiber model displays standard characteristics of a self-avoiding walk. First, the mean radius of gyration (R) over all trajectories of the non-methylated fiber (43.2 nm) exactly matches the mean R calculated from the formula  = (N∗b)1/2/6 = 43.3 nm, where N = 50 and b = 15 nm. Second, to test the dynamic properties of the fiber, we compared relaxation times () for the slowest mode calculated either directly from the simulation autocorrelation function or from the Rouse model, both ways resulting in comparable hundred-millisecond-scale relaxation times (i.e., equal to respectively 0.240 and 0.265 s; see details in supplemental information).

Table 3

Association/dissociation probabilities of the interacting particles utilized in the simulations

Particle 1	Particle 2	Reaction	Probability
H3_K9me3	HP1_cd	Association	0.000001
HP1_csd	HP1_csd	Association	0.000001
H3_K9me3	HP1_cd	Dissociation	0.00000001
HP1_csd	HP1_csd	Dissociation	0.00000001

An important example of a kilobase-scale heterochromatin represented by our model is the Oct4 promoter region, which, along with a group of similarly behaving regulatory elements, is subject to a so-called self-bounded H3K9 methylation (Hathaway et al., 2012). In particular, it features a nucleation site from which methylation spreads toward the boundaries where it decays because of a competition of opposing histone methyltransferase and demethylase activities. In such regions, the H3K9me3 PDF has a bell-like shape and is centered on the nucleation site (Hathaway et al., 2012). The H3K9 methylation process can stochastically occur rapidly but may take up to 5–6 days to reach an equilibrium H3K9me3 PDF (Hathaway et al., 2012). Hence, to model the Oct4 promoter in the course of the methylation spread, three ensembles of methylated fiber in random 3D conformations were generated using a gaussian PDF centered on the central nucleosome: 1) an ensemble of 25 non-methylated fibers (corresponding to day zero of the H3K9me3 spread), 2) an ensemble of 25 approximately half-methylated fibers (corresponding to an average day-2 H3K9me3 PDF), and 3) an ensemble of 25 almost fully methylated fiber (corresponding to an average day-5 H3K9me3 PDF) (Figure 1A)(Hathaway et al., 2012). Here and in the following sections, we use the terms "methylation" and "H3K9me3" interchangeably, although a methylated fiber may also feature populations of mono- and di-methylation marks to which HP1 chromodomain is able to bind (Al-Sady et al., 2013). Hence, in our model, only one collective methylation state is used; for the sake of convenience, we label it as H3K9me3.

Figure 1

HP1-mediated chromatin fiber condensation

(A) Methylation patterns used in simulations corresponding to days 0, 2, and 5 of the heterochromatinization process; shades of gray reflect the probability for a nucleosome at this position to have an H3K9me3 mark.

(B and C) Metrics used monitor the physical state of the chromatin fiber: radius of gyration () and its moving average ().

(D) Time charts for R from a sample trajectory (top) and for N, N, N, and L averaged over an ensemble of 25 trajectories (bottom) of a nearly fully methylated chromatin fiber; Structures shown exemplify typical conformations of the chromatin fiber at different stages of the condensation process; See Figure S1 for the remaining 74 time charts at 3 different methylation level.

(E) Magnified examples of fiber conformations (graytube) in a transition and fully condensed states; H3 histone tails are shown as blue sticks and HP1 chromo and chromo-shadow domains, as red spheres. See also Figure S1.

Proteins and protein-protein interactions

Two types of proteins – HP1 and histone H3 – are explicitly represented in the model (see Tables 1, 2, and 3 for the protein particles’ parameters). There are 102 copies of HP1, each modeled as a chain of two particles – chromodomain (CD) and chromo-shadow domain (CSD). This number of copies in a 400 nm box corresponds to a concentration of 2.65 μM, which is within the experimentally determined range of 1–10 μM (Canzio et al., 2011; Müller et al., 2009; Müller-Ott et al., 2014). The CD and CSD particles of HP1 were allowed to freely move within a distance range between R+R and D, where R and R (equal to 1 nm each) are radii of CD and CSD domains, respectively, and D (equal to 4 nm) is the standard distance maintained by a disordered linker connecting CD and CSD with an exponential damping beyond D. The damping multiplier k in exp(-kΔd) was set to 0.3. We made sure that the time step of 1 μs would not result in a systematic overstretching of the CD-CSD tether. PDF for CD-CSD distances over the ensemble of simulations (Figure S5B) is typical of a random polymer (as is the case for the actual CD-CSD linker). Histone H3 is modeled as a two-particle chain consisting of an H3 core (a part of the nucleosome body) and an H3K9 residue (that can be in either a non-methylated or a methylated state). The H3 core and H3K9 are maintained within a distance range of 1–4 nm with exponential damping outside that region exp(-kΔd), where k = 0.3.

The protein-protein interactions in the simulated system include HP1 CSD homodimerization and binding of HP1 CD to H3K9me3. The association (p) and dissociation (p) probabilities were used to parameterize stochastic processes “carried” by interacting particles (see Tables 1, 2, and 3 for a full list of interaction parameters). Association of two particles into a body occurs when (i) one associable particle (e.g., HP1 CSD) finds itself in the same 10 nm cubic cell with another and (ii) stochastic process on one of the respective particles produces a value smaller than p. After a new body is formed, in each time step, it may split into the two constituent particles if the stochastic process on one of them produces a value smaller than p. Here, the p parameters were determined as inverse residence times of HP1 on heterochromatin and HP1-dimer lifetimes that were experimentally determined to be on the order of minutes (Canzio et al., 2011; Müller et al., 2009; Müller-Ott et al., 2014; Teif et al., 2015). The p values were estimated as , where [F] is the concentration of free binding partner in the 10 nm reaction cube (∼1,500 μM) and k/k is the ratio of experimentally determined association and dissociation rates (e.g., in one study (Buecker et al., 2015) k and k for HP1α on heterochromatin were measured as 0.41 and 0.15, respectively). Eventually, given significant experimental error associated with these measurements, p and p values for both merge reactions were rounded to 10−5 and 10−8, respectively. To test the choice of the p and p values, we ran a series of simulations with only HP1 proteins STAR Methods and measured the equilibrium binding constant as K = [HP1d]/[HP1m/2]2, where [HP1d] is concentration of HP1 dimers and [HP1m] is concentration of free HP1 monomers. K values of 5.93 ± 2.08 μM were obtained, i.e., close to the respective experimental values of 1–5 μM (Mendez et al., 2011; Müller-Ott et al., 2014; Sanulli et al., 2019a).

Results

HP1 induces highly condensed chromatin fiber state

We first sought to investigate how HP1 binding can influence chromatin fiber structure. Because the fiber composition is uniform, it putatively behaves as an ideal polymer chain. Hence, radius of gyration (R) (Figure 1B), a standard polymer metric, was used to monitor the fiber behavior. We also deemed plausible that CSD dimerization and HP1-H3K9me3 interaction may produce four distinct species of molecular complexes (Figure 1C): HP1 dimers (the number of such complexes will be referred to as N), HP1 monomers bound to H3K9me3 (N), HP1 dimers bound to H3K9me3 (N), and HP1 dimers bridging two nucleosomes (N). Therefore, metrics reflecting the numbers of HP1 proteins involved in all of the aforementioned complexes were analyzed to approximate medium and long-range chromatin fiber interactions. Of note, although R was monitored in all three ensembles of simulation runs (for unmethylated, partially methylated and fully methylated fiber), HP1-chromatin complexes cannot be formed by design in simulations with unmethylated fiber. Figure 1D shows an R time chart for one of the 25 simulations with fully methylated fiber displaying, where the fiber displays a steady condensation toward a compact state with R between 16 and 21 nm. In this particular example, HP1 molecules start bridging adjacent nucleosomes after ∼30 s with the first bridge between sequence-remote nucleosomes formed shortly after 120 s. This first sequence-remote bridge brings all nucleosomes on the fiber closer together and hence facilitates encounters in 3D space for all other sequence-remote nucleosomes, thus further increasing the odds of forming more HP1-mediated bridges. As bridging progresses, we observe a clear trend on decreasing R. In the beginning of the simulation (Figure 1D), R shows a moving average of ∼45 nm, with individual R varying between 23 to 85 nm. These R values reflect with a high precision the behavior of an ideal polymer (a theoretically predicted mean R value is 44 nm (Rubinstein and Colby, 2003)). At its minimum, by the end of the simulation, the R moving average decreases to ∼18 nm corresponding to an approximately 15-fold decrease in nucleosome density (varying as the cube of R). At this point, sixteen HP1-mediated nucleosome-nucleosome bridges are formed with an average loop size varying between 17 to 19. In addition, Figure 1 shows sampled conformations of the fiber along the trajectory. In conformations that do not feature loops longer than 15–20 nucleosomes and having R over 40 nm, every single nucleosome is clearly accessible to interactions with any external entity such as transcription machinery. In contrast, fiber conformations displaying ten or more bridges, such as those in Figure 1E, do have significant numbers of sterically occluded nucleosomes.

We also monitored the collective trends in molecular interactions leading to the fiber condensation. As can be seen in the time chart reflecting quantities of molecular complexes formed in the system with fully methylated fiber, HP1-mediated nucleosome-nucleosome bridges are largely the dominant type of molecular complexes. An average of 15 HP1-mediated bridges by the end of a typical trajectory; however, at the same time, up to 5 single HP1 molecules are bound to H3K9me3 and up to 3 HP1 free dimers are formed. On the whole, up to 35 HP1 molecules (or 34% of HP1 available in the simulation box) are immobilized on heterochromatin (see the Discussion section for interpretation of these results in the context of putative mechanisms of heterochromatin formation).

HP1-induced chromatin condensation obstructs the conformational mobility of the fiber

Next, we studied how HP1-mediated nucleosome bridging affects the fiber’s conformational mobility. In particular, we examined how stable are the fiber conformations produced throughout the simulations in different fiber states (condensed or uncondensed). To this end, we calculated autocorrelation functions (AF) for simulations of the fully methylated, partially methylated, and unmethylated fiber. Pearson’s correlation coefficient (r) values were calculated for sets of the end-to-end distances R representing each trajectory (one snapshot out of 100,000 was exported to this end) shifted by lag times multiple of 0.01 s up to 100 s. The r values obtained were then averaged over the respective simulation ensembles. As can be seen in Figure 2A, the unmethylated fiber collectively behaves as an ideal chain showing zero autocorrelation after a lag time of approximately 2.5 s. The AF of the partially methylated fiber also falls sharply (within ∼2.5 s) to ∼0.09 correlation that suggests no significant relationships between the lagged trajectories. However, further decay to zero correlations takes ∼100 s (Figure 2A). This residual correlation must be because of a small number of HP1 bridges between the methylated nucleosomes adjacent to the center. Finally, the AF of the fully methylated fiber falls (within ∼2.5 s) to ∼0.25 correlation (i.e., significantly higher than for the partially methylated fiber), suggesting a stronger relationship even between largely lagged trajectories; it does not decay to zero within the lag time of 100 s. Two AFs for individual end-to-end trajectories show even higher levels of lasting correlations (Figure 2A). The latter result suggests that multiple lasting HP1 bridges between sequence-remoted nucleosomes result in a slower-evolving disordered state of the fiber also observed in previous studies (Michieletto et al., 2017b; Shi et al., 2018). A compelling way to see how conformational mobility of the fiber evolves throughout the simulation is through averaging distance maps over an extended period of time (0.1 s), so that only long-lasting inter-nucleosome contacts would appear on the map. For instance, the averaged distance matrix after approximately 85 s of simulation (Figure 2C) shows no persistent contacts between sequence-remote nucleosomes even though specific instantaneous conformations (Figure 2D) feature such contacts. As the fiber condenses, first off-diagonal elements appear in the averaged distance matrix (e.g., Figure 2E). Later, for more condensed fiber states, averaged distance matrices (e.g., Figure 2F) show a complex pattern of multiple sequence-remote inter-nucleosome contacts very similar to those seen for specific instantaneous conformations (e.g., Figure 2G). The distance matrices in Figures 2C–2G are reminiscent of contact maps resulting from Hi-C experiments. However, a typical bin size on a Hi-C map (>10 kpb) exceeds the size of our whole system (which is two orders of magnitude smaller than an average topologically associated domain (TAD)). Although there are similarities between the TAD and local heterochromatin formation processes (both are forms of polymer self-organization) there are also significant differences because of differing protein-mediated mechanisms involved.

Figure 2

Conformational mobility of the fully methylated chromatin fiber undergoing the HP1-induced condensation

(A) Autocorrelation functions (AF) for the end-to-end distances in simulations of the fully methylated, partially methylated, and unmethylated fiber with lag times ranging from 1 to 100 s (2.5 s in the insert chart); Three curves, rendered in thicker lines, represent the AFs averaged over the ensemble of 25 runs for the fiber of the same methylation level; Two more curves, rendered in thinner lines, represent the AFs for single runs showing the highest level of lasting autocorrelation.

(B) Time chart for R in a simulation featuring a broad range of compaction states.

(C–G) Examples of nucleosome-nucleosome distance matrices at times t indicated by dotted lines; (B)(E)(F) distance matrices with elements averaged over 0.1 second; (D)(G) instantaneous distance matrices; the dimensionality of the matrix (51x51) corresponds to the number of nucleosomes in the fiber, where each off-diagonal element is the distance between the i-th and j-th nucleosomes.See also Figures S1, S2, and S3.

Nondeterministic character of the fiber condensation

Relatively few simulations of the fully methylated fiber (e.g., those exemplified in Figures 1D and 3A) featured a steady sliding into the condensed state. Overall, a great variety of temporal R patterns was observed over the full ensemble of trajectories, most of which display three distinct chromatin states. First, the ideal-polymer behavior – with R mostly ranging between 40 – 45 nm but sometimes descending slightly below 35 nm – was observed in the beginning of every simulation of the methylated fiber (Figures 3, S1, and S2) and all along the trajectories of the non-methylated fiber (Figure S3). The ideal-polymer behavior was also seen along four entire trajectories of half-methylated fiber (Figures S2K and S2P). Besides obvious cases of non-methylated or non-HP1-bound fiber, ideal-polymer behavior can also be seen in HP1-bound fiber featuring short HP1-mediated loops. Such HP1-coated fiber has no or very few vacant non-dimerized chromatin-bound HP1 molecules available to bridge sequence-remote nucleosomes. Another clearly identifiable chromatin state corresponds to a highly condensed fiber with multiple (up to 20) HP1 bridges between sequence-remote nucleosomes and is characterized by low values (e.g., Figures 3A, 3C, 3D, S1B, S1C, S1F, S1H, S1L, S1P, S1Q, and S1X). This condensed state – that is, randomly folded, dense fiber (e.g., Figure 1E) showing slow conformational diffusion – is likely to cause gene repression, because most of its nucleosomes are inaccessible by transcription factors. In most simulations, condensation of the highly methylated fiber was apparently irreversible (Figure 3), although there are cases of a reverse transition from a relatively long-lived condensed state back to the ideal polymer state (Figures 3D, S1H, S1P, and S1S). The third observed fiber state corresponds to a transition from the ideal polymer to a highly condensed one. This state is characterized by fewer (less than ten) HP1 bridges between sequence-remote nucleosomes, between 25 and 35nm and standard error of ∼10 nm. It might be lasting, up to 500 s (e.g., Figure 3A) or shorter than 100 s (e.g., Figures 3C and 3D), probably depending on the fiber configuration at the moment when first HP1 bridges between sequence-remote nucleosomes are formed. In addition, as can be seen in Figure 3B, the transition state can be easily reverted to that of an ideal polymer.

Figure 3

Stochastic character of the fiber condensation

(A–D) Time charts for and for sample trajectories of a fully methylated chromatin fiber; The trajectories exemplify varying patterns of behavior observed for the fully methylated fiber: early progressive condensation with a high density state reached after ∼450 s of simulation (A), no full condensation throughout the simulation time (B), late rapid condensation after ∼700 s (C), and intermittent condensation in the intervals 120–180 s, 320–380 s, and 540–650 s (D).

(E) Ensemble average charts for a and in simulations of the non-methylated, half-methylated, and fully methylated fiber (25 simulations/each).

(F) Ensemble average charts for a in simulations of the half-methylated and fully methylated fiber (25 simulations/each); no bridges were formed in simulations of non-methylated fiber.See also Figures S1, S2, and S3.

The observed variability in temporal patterns of the methylated fiber behavior might be one of the factors helping to rationalize an intriguing and yet ununderstood aspect of the H3K9me3/HP1-induced heterochromatin – a broad cell-to-cell variability in lag times preceding heterochromatin-induced gene silencing after triggering the heterochromatin formation process (Bintu et al., 2016; Hathaway et al., 2012).

Another outcome of this study is that the fiber condensation is clearly H3K9me-dependent. First, as expected, none of the simulations of the non-methylated fiber shows any signs of condensation. Second, the highly condensed state (with  < 22) was observed only in simulations of the fully methylated fiber (in 19 out of 25 simulations for longer than 120 s). Third, half-methylated fiber, in every simulation, was switching back and forth between a semi-condensed and a random-polymer state. A clear role of the methylation level can be seen in ensemble-averaged time charts in Figures 3E and 3F. In particular, the fully methylated fiber continues to condense to  = 25 nm, whereas the half-methylated fiber plateaus at 34 nm after ∼300 s of simulation (Figure 3E). Similarly, the ensemble-averaged number of HP1 bridges steadily grows to ∼16 during the whole simulation time in simulations of the fully methylated fiber, whereas it varied between 6 and 10 for the half-methylated fiber (Figure 3F). The stability of the fiber condensation is likely because of its 3D conformation and the topology of the HP1 bridging. When the fiber is highly entangled, a few bridges between sequence-remote nucleosomes would “freeze” it in a stable densely packed state, so that dissociation of a single HP1 bridge would not destabilize the structure until a new bridge is formed. Thus, the described results suggest that a near-full methylation of the fiber is a prerequisite for its condensation.

Experimental repression of Oct4 by HP1

To validate how closely computational simulation resembled the properties of heterochromatin formation in cells, we also used the chemical induced proximity (CIP)-based Chromatin in vivo Assay (CiA) to visualize the kinetics of heterochromatin formation (Figure 4A)(Hathaway et al., 2012). The CiA system uses rapamycin to reversibly tether HP1α to an engineered Oct4 loci, which contains a Gal4 DNA binding array followed by an eGFP reporter replacing one allele of Oct4. Once the HP1α is docked, a flurry of activity transforms this allele in a start of euchromatin chromatin to a heterochromatin structure with all the properties of physiologic heterochromatin. Following initial chemical tethering of HP1α neighboring histones adjacent to the recruitment site, the CiA:Oct4 promoters have been demonstrated to lose active chromatin posttranslational modification (PTM, e.g., H3K4me3 and H3ac) while gaining signature heterochromatin marks like H3K9me3 (Hathaway et al., 2012). Here we evaluate the eGFP signal over an intricate time course to help validate the computational model, especially the timescale at which chromatin becomes repressed heterochromatin in individual cells. We measure this by both fluorescence microscopy and flow cytometry; thus, we can assess the heterochromatin formation process with single cell resolution reporting the number of GFP (−) cells over time.

Figure 4

Time course of experimental generation of heterochromatin by chemically tethering HP1 to CiA:Oct4 allele

(A) CIP-Rapamycin mediated tethering of HP1α to CiA:Oct4 promoter stimulated heterochromatin gene repression.

(B) Representative images of mESC during HP1 stimulated repression.

(C) Average GFP levels over time following HP1 recruitment∗.

(D) Number of GFP (−) cells∗. See also Figure S4.

∗Error bars reflect standard deviations calculated from triplicate sample analyses.

We treated the CiA cells with CIP-rapamycin over a time course of 5 days, measuring GFP signal at fixed intervals. During the first 12h of rapamycin treatment, cells remained mostly in a euchromatin state with high levels of GFP expression (Figure 4C). Cells with repressed Oct4 appear on day 1 after the rapamycin treatment, and the GFP signal in cells was reduced (Figure 4B,C). This trend of GFP repression continues as CIP-HP1 recruitment continued to day 5 (Figure 4C). The reporter measurement demonstrated a dramatic reduction between 36 h and 2 days, which represented a reflection point where the bulk cell population shifted to lower GFP signal (Figures 4C and 4D). By day 3, the majority of the cells reached the GFP silenced state and as the rapamycin treatment continues the silenced states are more uniform (Figures 4C and 4D and Figure S4A). During each of the rapamycin treatment times, cell colonies behaved similarly, with the exception of a few colonies that remained GFP positive even after 5 days of CIP-HP1 treatment (Figure 4B).

Discussion

Although there are decades worth of experimental data regarding the molecular markers of heterochromatin, computational effort to gain structural insight into heterochromatinization process is relatively recent (Barbieri et al., 2012; Brackley et al., 2016; Jost et al., 2014; MacPherson et al., 2018; Michieletto et al., 2017a, 2017b; Nuebler et al., 2018; Shi et al., 2018). Ideally, a computational model of heterochromatin is expected to rationalize its known traits (Allshire and Madhani, 2017) including (i) a structural mechanism by which heterochromatinization of regulatory regions may repress gene activity, (ii) stochastic nature of gene repression (Bintu et al., 2016; Elgin and Reuter, 2013), (iii) the respective timescale (Hathaway et al., 2012), (iv) condensed state of heterochromatin and phase separation (Larson et al., 2017; Strom et al., 2017), and (v) roles for all molecular players in the process. Here, molecular simulations of the chromatin fiber reveal structural transformations undergone by the fiber upon H3K9 methylation and HP1-mediated bridging in short, kilobase-scale chromatin regions. Although our model is too simple to account for all the aforementioned criteria and provide an accurate timescale, it does offer a structural rationale for several key features of the heterochromatinization process. In particular, multiple HP1-CSD-mediated bridges between sequence-remote nucleosomes induce the fiber entanglement into stable highly condensed structures. In those "frozen" random structures, a significant proportion of nucleosomes are fully buried inside the fiber clew, likely becoming inaccessible by the transcription machinery. Currently available experimental evidence suggested two possible mechanisms of HP1-driven heterochromatinization and phase separation – (i) through intrinsic propensity of HP1 to form phase-separated droplets (Larson et al., 2017; Sanulli et al., 2019b; Shakya et al., 2020) and (ii) HP1-mediated bridging. A recent study by Erdel et al. (2020) brought compelling data in favor of the latter hypothesis. In particular, the study (Erdel et al., 2020) has shown that the HP1 concentration in heterochromatin (of ∼3 μM (Müller-Ott et al., 2014)) is well below the half-saturation concentration for in vitro droplet formation of ∼40–45 μM. Moreover, an experimental study by (Nickels et al., 2021) suggests that heterochromatin propagates in sudden bursts at different time intervals with HP1 H3K9me reader responsible for this bursting behavior. In agreement with this recent evidence, our simulations demonstrate that a stable highly condensed chromatin state for kilobase-size chromatin regions can be achieved even at a typical cellular HP1 concentration of 2.65 μM, after broadly varying lag times.

There is another notable aspect of the heterochromatin formation process that might be, at least partially, explained by our model. As can be seen in Figure 4C, cells with repressed Oct4 start appearing on day 1 after HP1 tethering to the Oct4 promoter region is triggered in CiA cells. The share of the repressed (i.e., GFP(−)) cells rapidly grows until day 5 when it approaches a plateau, although a small fraction of GFP(+) cells persists even on day 5. In addition, it has previously been shown that the H3K9 methylation density of the Oct4 promoter region increases with time between days 0 and 5 of the heterochromatinization process (Hathaway et al., 2012). Accordingly, in simulations, we show that because of a high entropic penalty, HP1-mediated bridging of sequence-remote nucleosomes is an extremely rare event occurring on a minute scale. A state where several long-range HP1 bridges coexist may take multiple minutes to occur, depending on the initial fiber conformation. Even longer time may be required to achieve a stable condensed state in which every bridge is redundant, hence allowing a conformation to survive breaking of one or more bridges. It is remarkable that such a simple model that neglects important aspects of the systems energetics, such as explicit crowding or nonspecific protein-protein and protein-DNA interactions, was still able to provide a comprehensive spatial and temporal view of the heterochromatin formation process. This suggests that the process is largely determined by entropic costs of adopting system’s configurations favorable for formation of quaternary H3-HP1 complexes. Therefore, a Monte Carlo simulation model combining an accurate parameterization of protein-protein binding free-energies (from experimental data) with realistic configurational sampling provides a comprehensive structural interpretation for the currently accumulated evidence on heterochromatin.

To sum up, this study proposes a structural depiction of HP1-mediated heterochromatin formation on short regulatory regions that provides a sound interpretation of several open questions regarding the mechanisms of heterochromatin formation. The HP-mediated heterochromatin is a major "barrier to cell fate changes" (Becker et al., 2016) and a better understanding of its physical roots may have significant implications for basic and applied chromatin research.

Limitations of the study

It is of note that in our designed system with pre-methylated nucleosomes the chromatin condensation is a minute-scale process. In live cells, there are two other important processes that would further slowdown heterochromatin formation. First, a permanent competition of methyltransferases and demethylases permanently redesigns the methylation pattern (Hodges and Crabtree, 2012). Hence, broken HP1 bridges cannot always be readily re-formed when one of the bridged nucleosomes gets demethylated. Second, although the diffusion of single particles in our simulations already implicitly accounts for molecular crowding, the presence of adjacent and sequence-remote chromatin fiber remains unaccounted for. Although difficult to quantify, both presence of writer/eraser enzymes and fiber crowding would likely be able to further slowdown the fiber dynamics to the experimentally observed day-scale times. Furthermore, one interesting hypothesis can be inferred by combining our computational and experimental results. According to experimental data (Figures 4B–4D), there is a significant surge in the number of cells with the repressed target gene on day 2 of the experiment, whereas virtually no reduction in GFP signals occurred during day 1. At the same time, our computational simulations showed that the chromatin fiber has to be fully methylated to adopt a highly condensed state putatively associated with a repressed downstream gene. Therefore, it might be deduced that a minimum time required for the full promoter methylation is about a day.

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, Dmitri Kireev (dmitri.kireev@unc.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Cell culture

Mouse embryonic stem (ES) cells CiA:Oct4 infected with N118 and N163 plasmids (N118, nLV EF-1α-Gal-FKBPx1-HA-PGK-Blast; N163, nLV EF-1α-HP1α (CS)-Frbx2(Frb+FrbWobb)-V5-PGK-Puro) were grown in cell culture media containing high-glucose DMEM (Corning, 10013CV), 20% FBS(Gibco 26140–079), 10 mM HEPES pH 7.5(Corning, 25060Cl), NEAA(Gibco 11140050), Pen/Strep (Corning, 20002CL), 2-Mercaptoethanol(Gibco, 21985023), and 1:500 LIF conditioned media produced from Lif 1Ca (COS) cells. ES cells were cultured essentially as previously described (Hathaway et al., 2012). CiA:Oct4 with N118 and N163 were selected with puromycin (invivogen, ant-pr) and blasticidin (invivogen, ant-bl) for three days, and remove selection one day prior to rapamycin (LC laboratories, R-5000) treatment. For the time course experiment, CiA:Oct4 with N118 and N163 were plated in 80,000 cells per well in a 6-well plate for the time course experiment. Rapamycin from a 10 μM stock dissolved in ethanol was added to media at 3 nM concentration for the time course experiment. CiA:Oct4 cell lines were generated in a previous study (Hathaway et al., 2012). No animals were used in this work.

Method details

Materials and methods

Flow cytometry

Analysis of GFP expression was conducted with Attune Nxt machines (Thermo Fisher). All samples were analyzed in biological triplicate and standard devation was used to report error. Each individual sample was grow in one well of a six well plate. Cells were resuspended in 300ul of FACs buffer (1xPBS, 1% FBS) then and analyzed on the Attune NXT machine. Flow analysis was conducted using FlowJo cells were gated based on Ssc vs. Fsc, then Fsc-H vs. Fsc-A, then a histogram of eGFP fluorescence, and finally reported as a GFP(+) based on a threshold of GFP(+) cells.

Computational methods

Simulations

The simulated systems (as described in the Model section of this report) consisting of a 51-nucleosome chromatin fiber (featuring three different methylation patterns) and 102 copies of the HP1 CD-CSD chains were simulated in a 400 nm cubic box. The three systems differing by their methylation states were as follows: (i) no methylation, (ii) approximately half of 102 H3K9 particles marked as trimethylated, and (iii) almost all H3K9 particles marked as trimethylated. These three methylation states correspond to averaged methylation profiles experimentally observed on days 0, 2 and 5 of the heterochromatin formation process (Hathaway et al., 2012). The methylation state of each H3K9 particle was calculated based on the Gaussian probability density distributions centered on the nucleosome #26. Each simulation system was sampled 25 times to create ensembles of simulations. The 102 HP1 CD-CSD chains were uniformly and randomly distributed over the simulation box. To maintain a central location for the chromatin fiber in the box, all particle locations were adjusted relative to the central histone as the origin. Particles reflected off of the walls of the box. The time step utilized was 1 μs. Frames were recorded every 1,000 steps.

Tables 1, 2, and 3 provide a full set of parameters utilized in the simulations (see the Model section for detailed explanation on how the parameters used are linked to experimental data with respective literature citations).

Structural metrics

The radius of gyration (R) for the chromatin fiber was calculated usingwhere R is the distance between each element to the center of mass, m is the mass, and n is the total number of elements. R, N, N, N, and L were calculated from the exported trajectory snapshots using the Pipeline Pilot software (3dsbiovia.com). On average, ∼700,000 frames per trajectory were processed. Since the latter number is far beyond the resolution of the time charts in Figures 1 and 3 and Figures S1, S2, and S3, the individual points on the charts represent average values over 100 frames. R for each 0.1 s segment (10,000 frames) was calculated with the time average of the previous 120 s.

Relaxation of the non-methylated chromatin fiber

The relaxation time of the slowest mode for the Rouse model of an ideal chain in solution can be written as (Avdoshenko et al., 2017; de Gennes and Witten, 1980; Doi et al., 1988):where is the mean-square end-to-end distance usedas a measure of the polymer size and is the translational diffusion coefficient. For single nucleosomes or short chromatin regions (of ∼3 kbp), was experimentally estimated as ∼0.1μm2s(Amitai, 2018; Amitai et al., 2017; Wachsmuth et al., 2016). Given the mean end-to-end distance observed in our simulations of the non-methylated fiber of ∼100 nm, after applying Equation 4 we obtain an estimate of .

Furthermore, we assessed from the autocorrelation functions for the end-to-end distances. The mean autocorrelation function with a time lag step of 1 millisecond (ms) was calculated as an average of autocorrelation functions for 25 individual simulations of the non-methylated fiber (see supplemental files for calculated end-to-end distances and autocorrelation functions). The decline in autocorrelation by a factor of e (i.e., ∼2.72) occurred after a lag time of ∼240 ms.

Equilibrium constants of HP1 dimerization

Equilibrium binding constants for HP1 chromo-shadow domains were estimated through simulations of HP1 only in the simulation box (i.e., without the fiber). A total of 5 simulations were run at varying concentrations of HP1 (see below table). The binding constants were calculated as , where and are concentrations of respectively HP1 dimer and free monomer in a 400 × 400 × 400 nm box (see Table S1 for resulting KD values) at 1 s intervals over a minute-scale period (after 70 to 130 s equilibration times). The resulting mean K over all simulations and standard deviation are 5.93 ± 2.08 μM. As expected, the K values do not show clear correlation with the respective concentrations of HP1 although the constants at two highest concentrations are somewhat lower than others.

Quantification and statistical analysis

For flow cytometry analysis, three biological replicates of each sample point of the CiA-Oct4 experiments were conducted. Fluorescence was measured using the Attune Flow cytometry machine (ThermoFisher, Inc). The Flow cytometry data were analyzed using FlowJo software (BD, Inc). In each experiment, the live and single cell population was separated then cells were analyzed based on GFP expression, the average as well as standard deviation were calculated using Prism (GraphPad Software, Inc.). Graphs were generated with Prism (GraphPad Software, Inc.).

RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and proteins
Blasticidin	Invivogen
Puromycin	Invivogen
Rapamycin	LC laboratories	R-5000
Deposited data
Properties of the Chromatin Fiber sampled from MB simulations	This study	https://doi.org/10.17632/t4dvkhy4vy.2
Software
Matlab R2019a	MathWorks	www.mathworks.com
Molecular Biosystems	This study	N/A
Maestro molecular modeling suite 2017-1	Schrödinger	www.schrodinger.com
PyMol 2.1.1	Schrödinger	www.schrodinger.com
Pipeline Pilot (data processing software)	BIOVIA	www.3ds.com
FlowJo	FlowJo	www.flowjo.com

Effective binding constants for HP1 dimerization

Total HP1		Free monomers		Dimers		KD, μM
Count	[μM]	Count	[μM]	Count	[μM]
51	1.32	46.97	1.22	2.01	0.052	7.1
75	1.95	67.67	1.76	3.66	0.095	8.1
102	2.65	87.78	2.28	7.11	0.184	7.03
151	3.92	106.86	2.77	22.07	0.573	3.36
204	5.29	140.79	3.65	31.61	0.82	4.07

Counts and concentrations of HP1 molecules (total, free monomers, and dimers) and equilibrium binding constants for HP1 chromo-shadow domains.