- 1College of Life Sciences, Zhejiang University, Hangzhou, China
- 2Key Laboratory of Structural Biology of Zhejiang Province, School of Life Sciences, Westlake University, Hangzhou, China
- 3Westlake Laboratory of Life Sciences and Biomedicine, Hangzhou, China
- 4Institute of Biology, Westlake Institute for Advanced Study, Hangzhou, China
The organization of genomes in space and time dimension plays an important role in gene expression and regulation. Chromatin folding occurs in a dynamic, structured way that is subject to biophysical rules and biological processes. Nucleosomes are the basic unit of chromatin in living cells, and here we report on the effective interactions between two nucleosomes in physiological conditions using explicit-solvent all-atom simulations. Free energy landscapes derived from umbrella sampling simulations agree well with recent experimental and simulation results. Our simulations reveal the atomistic details of the interactions between nucleosomes in solution and can be used for constructing the coarse-grained model for chromatin in a bottom-up manner.
Introduction
Chromatin is highly compacted and condensed into the small space of the nucleus (Goloborodko et al., 2016b; Finn and Misteli, 2019). Human diploid genome contains about six billion DNA base pairs, and it will be approximately 2 m long if fully extended in a double helix (Fraser et al., 2015; Dans et al., 2016; Saurabh et al., 2016). In contrast, the size of the nucleus is only a few micrometers, and much remains to be understood that how chromatin fit into the nucleus in such a compacted and condensed way. The organization and dynamics of chromatin in the nucleus are found to be neither random nor stochastic, instead it is well defined and regulates the gene expression intricately.
The organization and folding processes of chromatin can be mainly divided into four layers: 1) The antiparallel double helical and right-handed B-DNA structure with each base pair rising up about 3.4 Å along the helical axis, which makes DNA very stable while it provides potential for the binding of proteins. 2) About 146 base pairs wrap around a histone octamer (also called histone core), including two copies each of the four histone core proteins (H2A, H2B, H3, and H4) in a left-handed superhelix, which constitute a nucleosome. Nucleosomes are the basic units of chromatin and provide controlled accessibility for DNA-binding proteins such as transcription machines and structural maintenance of chromosome (SMC) complexes. The height of a nucleosome is about 55 Å, so the formation of nucleosome compacts DNA by about 9 times as 146 base pairs × 3.4Å ÷ 55Å ≈ 9. Even though the DNA sequences binding to histone core are nonspecific, the binding affinity between some DNA sequences and histone core are higher than others (Field et al., 2008; Teif et al., 2012; Teif and Clarkson, 2019). About 75–90% of the genome are organized in the form of nucleosomes (Field et al., 2008). 3) Under the view of electron microscopes, chromatin appears as “beads on a string,” in which beads correspond to nucleosomes and the string between beads corresponds to the double helical DNA called linker. The string-like chromatin is shaped by loops (Alipour and Marko, 2012; Rao et al., 2014; Goloborodko et al., 2016a, Goloborodko et al., 2016b), topologically associating domains (TADs) (Dixon et al., 2012; Schwarzer et al., 2017; Krietenstein et al., 2020), and A/B compartments (Dekker et al., 2002; Zhao et al., 2006; Lieberman-aiden et al., 2009; Dekker et al., 2013). Because chromatins are fluid and dynamic, it is appropriate to use structural ensembles to describe the chromatin structure. While TADs and compartments fluctuate at a single-cell level, population-averaged TADs and compartments are tissue specific, meaning that their patterns are similar and conserved in one cell line and have high variability across different cell lines (Cheng et al., 2020; Contessoto et al., 2021). The variability of the organization and dynamics in this layer suggests direct impacts on the expression and regulation of genes (Finn and Misteli, 2019). 4) Each chromatin has its own territory in the nucleus, in which locus from the same chromatin have higher probability to localize together and form exclusive subregions (Stam et al., 2019). In the normal cell nucleus, euchromatins, which are gene-rich and transcriptionally active segments, cluster together in the interior, whereas heterochromatins, which are gene-poor and silenced, cluster together in the nuclear envelope. Visualization in real time and simulation in silicon provides a comprehensive understanding of the phase separation and chromosome territory (Liu et al., 2020; Oliveira Junior et al., 2020; Su et al., 2020). Overall, chromatin is organized in a highly hierarchical architecture. Each of these four layers is also highly dynamic, which provide potential to regulate important biological processes, in particular the gene expression.
The human genome is the blueprint of life consisting of more than 20,000 genes and millions of regulatory candidate elements (Dunham et al., 2012). Despite intensive efforts, it is far from complete to understand how these elements function and interact with each other in the spatial and temporal dimensions to regulate gene expression, as cells with identical DNA sequence can function differently. The three-dimensional structure and dynamics of chromatin play critical roles in bringing into physical proximity the regulatory elements with target genes across hundreds of kilobases or even megabase distance. The abnormal chromatin organization leads to the occurrence of diseases, in particular cancer (Goes et al., 2011; Meaburn et al., 2016). 3C-based methods such as Hi-C (Dekker et al., 2002; Lieberman-aiden et al., 2009) and FISH (Fraser et al., 2015; Fudenberg and Imakaev, 2017) are two mainstream experimental techniques to probe the organization and dynamics of chromatin. However, how to integrate and interpret the Hi-C and the FISH measurements remains challenging, and in some cases, they can even lead to contradictory results (Fudenberg and Imakaev, 2017).
Complementary to experiments, computer simulations provide unprecedented resolution to investigate the folding of chromatin. Polymer model theory can be used to study the organization and dynamics of chromatin fibers at genome scale. Coarse-grained (CG) models of nucleosome can be used to study the interaction and dynamics of nucleosomes array which is the local subregion of the chromatin fiber. All-atom molecular dynamics (MD) simulations of nucleosome can provide further detailed information at the atomistic level. Combining these models for multiscale simulations would be useful to enhance our understanding of chromatin folding processes. In the polymer model (Philip et al., 1993; Münkel and Langowski, 1998), chromatin is represented as polymers consisting of different monomers connected by harmonic bonds, and different persistence lengths are chosen according to the compaction ratio, for example, how many base pairs are coarse-grained into one monomer. The models also account for particular interactions related to the biological activities of chromatin, for example, loop extrusion and compartmental segregation. Recently, two sophisticated computational models with slightly different potential energy function formulas have been developed (Di Pierro et al., 2016; Fudenberg et al., 2016), and shed light on the mechanisms underlying the folding and organization of chromatin (Gibcus et al., 2018; Mirny et al., 2019; Stam et al., 2019; Banigan et al., 2020).
The quality of polymer model simulations depends critically on the accuracy of their potential energy functions, which are a summation of different pairwise interaction terms. Previous researches (Sanborn et al., 2015; Fudenberg et al., 2016; Goloborodko et al., 2016a; Goloborodko et al., 2016b) typically used grid search strategy to optimize the parameters in potential energy functions, which means trial-and-error until simulated properties match experiments. If there are four parameters in one energy function and each parameter have 10 grids, 10,000 sets of parameter combination are tried from which one will be selected, the one most consistent with the Hi-C contact map. In addition to grid search strategy, the maximum entropy principle (Zhang and Wolynes, 2015; Di Pierro et al., 2016) was also applied to derive the potential functions using the experimental contact map information as inputs. A prior knowledge is needed to derive models from the maximum entropy principle, where it is used by grid search strategies as the criterion to select the best parameter combination. An alternative way to construct coarse-grained models is using the “bottom-up” or “ab initio” approach. One can derive coarse-grained physical potentials from more detailed simulations, for example, explicit-solvent all-atom MD simulations (Noid et al., 2008; Li et al., 2016). Here, we present a first step toward constructing a multiscale polymer model for chromatin based on atomistic simulations of two nucleosomes.
Nucleosomes are fundamental units for chromatin folding and the main carrier of epigenetic marks, so we coarse-grain one nucleosome as one bead in the coarse-grained model. Since the high-resolution X-ray structure of the nucleosome was determined in 1997 (Luger et al., 1997), there have been many studies on the properties of nucleosomes (Portela and Esteller, 2010; Biswas et al., 2013; Feng and Li, 2017; Zhou et al., 2019). Recently, Funke et al. (2016) employed force spectrometer and single-particle electron microscopy to measure the forces and interaction profiles between nucleosomes by placing two nucleosomes close to each other in a variety of defined relative orientations and recording the frequency of their distances. Multiscale simulations revealed that increased secondary structure resulting from acetylation of H4 tail has an important effect on the rigidification and also impaired the interactions between stacked nucleosomes (Collepardo-Guevara et al., 2015). The effect of H4 tail on stabilizing the stacked nucleosome is also validated by a recent first atomistic simulation of stacked nucleosomes using conventional and steered MD simulations (Saurabh et al., 2016). Ishida and Kono (2017) explored the energy landscape of two stacked nucleosomes using umbrella sampling with nucleosomes restrained at a few distinctive orientations. Moller et al. (2019) evaluated the anisotropic energy landscape of stacked nucleosomes across a variety of parameters in configurational and environmental space using residue coarse-grained simulations. Lequieu et al. (2019) constructed a coarse-grained multiscale model of chromatin by mapping the energy landscape of stacked nucleosomes to a reduced coarse-grained topology. Interestingly, Spakowitz and co-workers investigated the effects of epigenetic modifications, especially methylation of lysine-9 of histone H3, on the organization and dynamics of chromatin using polymer model (MacPherson et al., 2018; Sandholtz et al., 2020).
Nucleosome positioning, referring to the location of the nucleosome dyads in linear DNA, regulates the accessibility of DNA to DNA-bound proteins (Portela and Esteller, 2010). Adjacent nucleosomes in sequence are connected by linkers, whose lengths are about tens of base pairs in eukaryotes. The average radius of folded protein is about a few nanometers (Milo et al., 2009), whereas the compartment structures are across megabases. So the process of compartmentalization is self-organized by the interactions between nucleosomes, rather than mediated by proteins. Modifications of nucleosome, such as DNA methylation and histone modification, change the properties of nucleosomes and thus alter their interaction landscapes. Theoretical modeling of the interactions between two nucleosomes, in particular its dependence on the distance and relative orientation of nucleosomes, are crucial to our understanding of chromatin folding.
In this work, we present all-atom MD simulations of two nucleosomes interacting with each other in physiologically relevant explicit-solvent environment and analyze their interaction landscapes in the context of the 30-nm chromatin model. We consider two simulation systems, one containing two linked nucleosomes, whereas the other one containing two unlinked nucleosomes proximal in space. Atomistic simulation results of these two systems could be used to determine the function forms and parameters in the model, such as the diameter of the beads, the strength of the harmonic bonds connecting beads, and the strength of the weak attraction between beads in different chromatin states. The manuscript is organized as follows: Details of MD simulations and trajectory analysis will be provided in the Methods section. In the Results section, unbiased MD simulation results will first be presented, followed by the potential of mean force (PMF) and coarse-graining calculations. The manuscript ends with a short discussion and conclusion.
Methods
System Setup
The interaction between two nucleosomes can be naturally divided into two types, one type representing the interactions between nucleosomes that are connected by the linker DNA and the other type representing the interactions between two spatially adjacent nucleosomes that are stacked together with no connecting linker. Accordingly, two simulation systems were set up as shown in Figure 1. The initial structures were taken from that of classical 30 nm fiber (pdb id: 6hkt) (Garcia-Saez et al., 2018). Although its resolution is relatively low (9.7 Å), it includes the information about the nucleosomal stacking and packing patterns which are important for building the simulated systems to model the interactions between nucleosomes. We note that the histone tails that are missing in the crystal structure are not modeled in the simulation systems.
FIGURE 1. (A) The template structure (pdb id: 6hkt) for the setup of simulation system includes six nucleosomes with a flat two-start helix with uniform nucleosomal stacking interfaces and an uncondensed nucleosome packing density. The LN system is marked by a red rectangle, whereas the ULN by a purple rectangle. (B) A snapshot of the LN system after 10 ns MD simulation. (C) A snapshot of the ULN system after 10 ns MD simulation (D) The distance, d, defined using the centers of nucleosomes and the angle, θ, defined using the superhelical axes of nucleosomes.
For the linked-nucleosomes (LN) system, two nucleosomes and their corresponding linker DNA base pairs were extracted from the experimental structure and solvated in an explicit TIP3P (Jorgensen et al., 1983) water box with dimensions of 463 Å × 142 Å × 110 Å (Figure 1B). CHARMM (Brooks et al., 1983; Brooks et al., 2009) was used to build the missing hydrogen atom coordinates and patch protein and nucleic acid terminals. 150 mM KCl was added with additional cations to neutralize the system. In total, the LN system contains 676,742 atoms which are composed of 16 protein chains, two DNA chains, 209,471 water molecules, 1,140 K+ ions, and 594 Cl− ions.
The unlinked-nucleosome (ULN) system contains two stacked nucleosomes plus 15 flanking base pairs extracted from the experimental 30 nm fiber structure. A similar procedure was used to solvate the system in a cubic water box with dimensions of 182 Å × 182 Å × 182 Å (Figure 1C). In total, the ULN system has 558,630 atoms composed of 16 protein chains, four DNA chains, 170,064 water molecules, 1,039 K+ ions, and 485 Cl− ions. Proteins and DNAs were modeled by the CHARMM36m protein force field (Huang et al., 2017) and the CHARMM36 nucleic acid force field (Hart et al., 2012), respectively.
Molecular Dynamics Simulations
MD simulations were performed using OpenMM (Eastman and Pande, 2010; Eastman et al., 2017) with the isothermal–isobaric (NPT) ensemble. Periodic boundary condition (PBC) was applied and particle mesh Ewald (PME) (Essmann et al., 1995) was used to compute all the nonbonded interactions with a real space cutoff at 9 Å. We noted that both electrostatics and van der Waals interactions are fully accounted for with no truncations, as the latter were treated by the recently developed LJ-PME method (Wennberg et al., 2015). All hydrogen-containing bonds were constrained by the SETTLE algorithm (Miyamoto and Kollman, 1992) and the Verlet integrator was used with a time step of 2 fs The temperature was maintained at 303.15 K using the Andersen thermostat (Andersen, 1980) with a damping coefficient of 1 ps−1. A Monte Carlo barostat (Aqvist et al., 2004) was used to maintain the pressure at 1 atm by attempting to change the box dimension every 25 steps. Both LN and ULN systems were simulated to 1 µs, and frames were saved every 5,000 steps (10 ps).
Umbrella Sampling and Potential of Mean Force Calculations
To investigate the free energy landscape between nucleosomes, advanced simulation techniques need to be employed. Here, we perform umbrella sampling simulations using the distance between the two nucleosomes as the reaction coordinate.
where the system Hamiltonian
Umbrella sampling simulations were carried out using OpenMM with a plugin for PLUMED (Tribello et al., 2014). The center of a nucleosome was defined as the center of mass of the phosphorus atoms of DNA wrapped on the nucleosome histone cores (LN system) or the center of mass of protein Cα atoms of the histone cores (ULN system). The reaction coordinate d was then defined as the distance between the centers of nucleosomes. For the LN system, 127 windows were used in total,
Results
Unbiased Molecular Dynamics Simulations
Unbiased MD simulations were carried out for 1,000 ns for both LN and ULN systems. For individual nucleosome, the binding between DNA and histone core was very stable as indicated by their root mean square deviations (RMSDs) being around 4 Å with respect to the initial structures (Supplementary Figure S3). The relative motion between nucleosomes, in contrast, was highly dynamic with respect to both their distance and orientation over the microsecond timescale (Figure 2). For nucleosomes connected by the linker, their distance d varied between 165 and 239 Å. The distance increased to 238 Å in the first 200 ns of the simulations and gradually decreased to 170 Å after 400 ns.
FIGURE 2. The distance and relative orientation between two nucleosomes in the LN (A, B) and ULN (C, D) systems from 1 μs unbiased MD simulations. Also shown are snapshots corresponding to different time points (E).
For unlinked nucleosomes, d increased almost linearly from the starting value of 60 to about 90 Å at the first 200 ns, and then fluctuated between 75 to 100 Å in the following 600 ns, and rose to 120 Å after 800 ns. This suggests that the equilibrium distance of two free nucleosomes in aqueous environment is much larger than the 62 Å in the crystal structure of 30 nm fiber, which are impacted by the crystal packing and the stacking of nucleosomes. Based on our simulations, the relative velocity between two nucleosomes can be estimated to be about 0.1 Å/ns or 0.01 m/s.
We also analyzed the relative orientation between the two nucleosomes, characterizing it using the angle between the superhelical axes of nucleosomes (Figure 1D). There were large variations of angle in 1,000 ns MD simulations for both LN and ULN systems. The length of the linker DNA in the linked-nucleosomes system is 40 base pairs, smaller than the persistence length of double helix DNA (about 150 base pairs). A weak correlation between orientation and distance d was observed in both the LN and ULN systems. In the LN system, the relative orientation angle has a tendency to decrease when d decreases. In the ULN system, the angle has a tendency to increase with the distance. The frequency of angle variation is significantly higher than that of distance. During 1,000 ns MD simulation, distance went down in the LN or up in the ULN system, whereas the angle vibrated up and down regularly.
Interaction Free Energy Landscapes
To study the effective interactions between nucleosomes in both systems, umbrella sampling simulations were performed to compute the PMFs as a function of distance between the nucleosomes. In general, the free energy profiles are shallow and flat, consistent with the hypothesis that nucleosomes are highly dynamic. The convergence of PMF is relatively good for both LN and ULN systems as indicated by the statistical uncertainties from the Monte Carlo bootstrapping calculations with WHAM (Grossfield, 2020–9) being low (Supplementary Figures S8–S11). For the LN system, the global minimum is found at
FIGURE 3. PMF calculated from umbrella sampling in the LN system along the distance between two nucleosomes (black solid line). The fitted harmonic potential function is shown as red dotted line.
The interaction between two unlinked nucleosomes in solution has a completely different free energy profile, featured by a strong repulsion wall at smaller distances, and a flat curve for larger d (Figure 4). It indicates that the interactions between ULN nucleosomes are repulsive. As histone tails are not included in our ULN simulation system, the repulsive interaction is consistent with previous simulations (Collepardo-Guevara et al., 2015; Saurabh et al., 2016; Ishida and Kono, 2017; Moller et al., 2019) and experiments (Dorigo et al., 2003; Gordon et al., 2005), which showed histone tails are crucial for holding stacked nucleosomes together. We note that the nucleosome has a disc-like shape with a height of 55 Å, so the sharp free energy rise below d = 60 Å is probably not caused by steric clash but instead by the unfavorable electrostatic interactions as each individual nucleosome is highly negatively charged. The asymptotic free energy was not fully resolved in our calculations, due to the limited size of the simulation box (182 Å × 182 Å × 182 Å). When the two nucleosomes are pulled away with d larger than 90 Å, interactions with their PBC images become possible such that the WHAM analysis is no longer valid.
FIGURE 4. PMF calculated from umbrella sampling in the ULN system along the distance between two nucleosomes (black solid lines). The fitted exponential and shifted Coulomb potential functions are shown as red dotted line and green solid line, respectively.
Similar to the LN system, we also analyzed the orientation of nucleosomes in each window. The relative orientation angle along the time series in each window is stable and vibrates regularly at a small interval (Supplementary Figure. S6) and correlates with the distance d (Supplementary Figure. S7). In the ULN system, orientation angle has a tendency to increase with d and has larger fluctuation at larger d, consistent with unbiased MD simulation in which angle increases with distance (Figures 2C,D). The computational PMF for the ULN system can be directly compared with a recent experimental measurement of the association of nucleosomes (Funke et al., 2016). By integrating two nucleosomes into a carefully designed and calibrated DNA origami-based force spectrometer, Funke et al derived the Boltzmann-weighted distance-dependent energy landscape for two nucleosomes interacting with each other. They observed three major features: a strong repulsion at distances smaller than 60 Å; a minimum located somewhere between 60 and 70 Å; and vanishing interactions at distances greater than 130 Å. The computational PMF shown in Figure 4 agrees well with strong repulsion at small distances and vanishing interactions at distal distance. No global minimum is found in our PMF due to the absence of histone tails in the ULN system, suggesting the importance of histone tails from an indirect perspective.
Constructing Coarse-Grained Potentials
With umbrella sampling, we are able to dissect the interactions between nucleosomes with atomistic details. On the other hand, the human chromatin contains about 28 million nucleosomes, so it is not feasible to use atomistic simulations to study chromatin folding in the foreseeable future. The free energy profiles we obtained can serve as the starting point to construct coarse-grained potentials to study the conformational dynamics of chromatin elements such as 30 nm fibers. Current polymer theory models for chromatin (Di Pierro et al., 2016; Fudenberg et al., 2016) are typically composed of five terms:
where the
If we construct a CG model in which each nucleosome is coarse-grained into one monomer, the
where k equals 0.01 kcal/mol/Å2 and
The free energy profile of the ULN system shows that the interactions between a pair of nucleosomes in distant distance are very weak (Figure 4). This suggests that the nucleosomes might not self-aggregate if there are no additional restraints such as histone tail effects, epigenetic modifications or binding of proteins such as cohesins or condensins.
we could determine the parameters to be
Another suitable potential energy function form to fit the PMF would be a shifted Coulomb potential.
However, fitting the PMF with such a shifted Coulomb potential leads to
Discussion and Conclusion
A spatially and temporally resolved understanding of chromatin organization is currently one of the central topics in molecular and cell biology. Here, we used classical MD simulations and enhanced sampling methods to study the basic element of chromatin, nucleosome. In particular, we studied the conformational dynamics and the free energy landscapes between two nucleosomes in aqueous environment. Not only the widely investigated stacked nucleosomes but also the linked nucleosomes are involved in our study. With umbrella sampling calculations, we obtained detailed free energy profiles for nucleosome pairs either free or connected by a linker DNA. Our simulation results compare favorably with a variety of experimental findings and MD simulations, in particular the PMF of the ULN system correlates well with a recent single-molecule force spectrometer measurement.
Considering that the configurations of nucleosomes are highly dynamic, the relative orientations of the stacked nucleosomes are often restrained to achieve quick convergence in previous simulation studies. In this work, the ULN (stacked) system was simulated freely without additional restrains at the cost of longer simulations, which makes the dynamics more natural and nucleosomes free to explore possible relative orientations. Both unbiased and biased simulations of unlinked nucleosomes show that electrostatic repulsion dominates distance distribution in the absence of histone tails. In addition to the ULN system, we also studied the dynamics and free energy profiles of linked nucleosomes. Results of LN and ULN systems could provide new insights into the organization and dynamics of larger chromatin elements such as 30 nm fibers.
Umbrella sampling is a useful technique to investigate the effective interaction between biological macromolecules. Recently, Lai et al (2020) studied the free energy profile between two identical DNA double helices using extensive umbrella sampling. In this study, we carried out 1 μs conventional MD simulations and found out that the relative motion between nucleosomes is very slow. We then performed accumulatively more than 3.8 μs umbrella sampling simulations on systems with more than half a million atoms to understand the slow- and large-scale dynamics between nucleosomes. This is only possible with recent advances in both hardware and software that utilize the graphics processing units (GPUs) for MD simulations. Simulation of the LN system (
We note that K+ ions were used as the cation in our simulation systems, whereas salt concentration and composition are more complicated in physiological environment and will impact the strength of nucleosome interactions. Mg2+ ions were known to play a crucial role in directly binding with and stabilizing nucleic acids. However, more accurate polarizable force fields might be needed to model Mg2+ ions (Huang et al., 2014; Lemkul and MacKerell, 2016; Walker et al., 2020). It would be interesting to investigate how the interaction landscapes between nucleosomes change with the different ion strength and composition. Another flaw of the current study concerns the limited simulation box size. For the ULN system, we are not able to obtain the asymptotic free energy over distance larger than 90 Å, as nucleosomes will interact through their images due to the PBC conditions. A larger and probably noncubic water box might be needed to overcome these limitations.
We carried out atomistic simulations of nucleosomes with the long-term goal to construct a coarse-grained potential for 30 nm chromatin structure to bridge the gap between all-atom model of nucleosomes and polymer model. This would constitute a “bottom-up” approach to determine the function forms for polymer models and to optimize their parameters. With the analytically fitted CG potentials presented in this work, one would already be able to simulate, for example, the classical 30 nm chromatin fibers to investigate the stability and dynamics of such nucleosome fiber arrays and then optimize the polymer model with this coarse-grained potential. For this purpose, we built a system consisting of 100 nucleosomes which had no open ends and was covalently bonded to itself through periodic boundary conditions. When we propagated the simulation system using derived potential functions (Eqs 3, 4) for
This highlights one of the current limitations of the preliminary CG model presented here. Conversion between different compacted states of 30 nm chromatin significantly depend on the histone tails, especially H4 tails, so simulations with histone tails need to be carried out to derive the corresponding free energy profiles. Epigenetic modifications often occur on the histone tails, changing the properties of nucleosomes such as charge, rigidification, and solvent exchange which make a big difference in the dynamics and organization of chromatin and gene expression. Considering the importance of histone tails for both nucleosome stacking interactions and higher level chromatin structures, we are currently performing similar simulations with full histone tails. However, this involves significantly more extensive MD simulations that are out of the scope of the current study, which represents a first step toward building up a CG model for chromatin with a bottom-up approach. The binding between the linker histone and linker DNA also plays important roles in chromatin compaction (Luque et al., 2014), especially the interactions between linked nucleosomes. Similar linked nucleosomes systems with linker histone should be studied in future work. On the other hand, the model is not applied to any practical biological problem, for example, loop extrusion and compartmentalization. In reality, the disk-like shape of nucleosomes and nonuniform distribution of their charges make interactions highly anisotropic, whereas interactions are assumed to be isotropic in our coarse-grained model.
As for the next step, we will perform similar umbrella sampling simulations to study how epigenetic modifications on nucleosomes change the interaction landscapes and derive the corresponding potential parameters. Modifications of interests include acetylation of lysines on histone H4, mono-, di-, and trimethylation of lysine 27 on histone H3 and monoubiquitination of lysine 119 on histone H2A. Parameterization for these epigenetic modifications would allow us to construct a transferable computational model for chromatin. We plan to construct a coarse-grained model at nucleosome resolution, aiming to make the model as simple as possible so that it can be used to simulate one whole chromosome even genome which consists of millions of nucleosomes. Four or more beads (monomers) will be needed to intimate the anisotropicity of one nucleosome and the highly negative charge will be included simultaneously in the future. Ultimately, the polymer model will be optimized using “bottom-up” strategy. Such a model is expected to be useful in understanding development-related and disease-related chromatin dynamics, for example, how the binding of Polycomb repressive complex 1 (PRC1) and 2 (PRC2) (Comet et al., 2016) induces the condensation of chromatin.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author Contributions
JH conceived the project, CZ performed simulations. JH and CZ analyzed the data and wrote the paper together.
Funding
The work is supported by Zhejiang Provincial Natural Science Foundation of China (Grant No. LR19B030001) and National Natural Science Foundation of China (Grant No. 21803057), Westlake Education Foundation and Tencent Foundation.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Westlake University Supercomputer Center for computational resource and related assistance.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.624679/full#supplementary-material.
References
Alipour, E., and Marko, J. F. (2012). Self-organization of domain structures by DNA-loop-extruding enzymes. Nucleic Acids Res. 40, 11202–11212. doi:10.1093/nar/gks925
Andersen, H. C. (1980). Molecular dynamics simulations at constant pressure and (or) temperature. J. Chem. Phys. 72. doi:10.1063/1.439486
Aqvist, J., Wennerström, P., Nervall, M., Bjelic, S., and Brandsdal, B. O. (2004). Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm. Chem. Phys. Lett. 384, 288–294. doi:10.1016/j.cplett.2003.12.039
Banigan, E. J., van den Berg, A. A., Brandão, H. B., Marko, J. F., and Mirny, L. A. (2020). Chromosome organization by one-sided and two-sided loop extrusion. eLife 9, e53558. doi:10.7554/eLife.53558
Biswas, M., Langowski, J., and Bishop, T. C. (2013). Atomistic simulations of nucleosomes. Wires Comput. Mol. Sci. 3, 378–392. doi:10.1002/wcms.1139
Brooks, B. R., Brooks, C. L., Mackerell, A. D., Nilsson, L., Petrella, R. J., Roux, B., et al. (2009). Charmm: the biomolecular simulation program. J. Comput. Chem. 30, 1545–1614. doi:10.1002/jcc.21287
Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S., and Karplus, M. (1983). Charmm: a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4, 187–217. doi:10.1002/jcc.540040211
Cheng, R. R., Contessoto, V. G., Lieberman Aiden, E., Wolynes, P. G., Di Pierro, M., and Onuchic, J. N. (2020). Exploring chromosomal structural heterogeneity across multiple cell lines. eLife 9, 1–21. doi:10.7554/eLife.60312
Collepardo-Guevara, R., Portella, G., Vendruscolo, M., Frenkel, D., Schlick, T., and Orozco, M. (2015). Chromatin unfolding by epigenetic modifications explained by dramatic impairment of internucleosome interactions: a multiscale computational study. J. Am. Chem. Soc. 137, 10205–10215. doi:10.1021/jacs.5b04086
Comet, I., Riising, E. M., Leblanc, B., and Helin, K. (2016). Maintaining cell identity: PRC2-mediated regulation of transcription and cancer. Nat. Rev. Cancer 16, 803–810. doi:10.1038/nrc.2016.83
Contessoto, V. G., Cheng, R. R., Hajitaheri, A., Dodero-Rojas, E., Mello, M. F., Lieberman-Aiden, E., et al. (2021). The Nucleome Data Bank: web-based resources to simulate and analyze the three-dimensional genome. Nucleic Acids Res. 49, D172–D182. doi:10.1093/nar/gkaa818
Dans, P. D., Walther, J., Gómez, H., and Orozco, M. (2016). Multiscale simulation of DNA. Curr. Opin. Struct. Biol. 37, 29–45. doi:10.1016/j.sbi.2015.11.011
Dekker, J., Marti-Renom, M. A., and Mirny, L. A. (2013). Exploring the three-dimensional organization of genomes: interpreting chromatin interaction data. Nat. Rev. Genet. 14, 390–403. doi:10.1038/nrg3454
Dekker, J., Rippe, K., Dekke, M., and Klechner, N. (2002). Capturing chromosome conformation. Science 295, 1306–1311. doi:10.1126/science.1067799
Di Pierro, M., Zhang, B., Aiden, E. L., Wolynes, P. G., and Onuchic, J. N. (2016). Transferable model for chromosome architecture. Proc. Natl. Acad. Sci. U.S.A. 113, 12168–12173. doi:10.1073/pnas.1613607113
Dixon, J. R., Selvaraj, S., Yue, F., Kim, A., Li, Y., Shen, Y., et al. (2012). Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature 485, 376–380. doi:10.1038/nature11082
Dorigo, B., Schalch, T., Bystricky, K., and Richmond, T. J. (2003). Chromatin fiber folding: requirement for the histone H4 N-terminal tail. J. Mol. Biol. 327, 85–96. doi:10.1016/S0022-2836(03)00025-1
Dunham, I., Kundaje, A., Aldred, S. F., Collins, P. J., Davis, C. A., Doyle, F., et al. (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57–74. doi:10.1038/nature11247
Eastman, P., and Pande, V. (2010). OpenMM: a hardware-independent framework for molecular simulations. Comput. Sci. Eng. 12, 34–39. doi:10.1109/MCSE.2010.27
Eastman, P., Swails, J., Chodera, J. D., McGibbon, R. T., Zhao, Y., Beauchamp, K. A., et al. (2017). OpenMM 7: rapid development of high performance algorithms for molecular dynamics. PLoS Comput. Biol. 13 (7), e1005659. doi:10.1371/journal.pcbi.1005659
Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. doi:10.1063/1.470117
Feng, Z.-X., and Li, Q.-Z. (2017). Recognition of long-range enhancer-promoter interactions by adding genomic signatures of segmented regulatory regions. Genomics 109, 341–352. doi:10.1016/j.ygeno.2017.05.009
Field, Y., Kaplan, N., Fondufe-Mittendorf, Y., Moore, I. K., Sharon, E., Lubling, Y., et al. (2008). Distinct modes of regulation by chromatin encoded through nucleosome positioning signals. PLoS Comput. Biol. 4, e1000216. doi:10.1371/journal.pcbi.1000216
Finn, E. H., and Misteli, T. (2019). Molecular basis and biological function of variability in spatial genome organization. Science 365 (6457), eaaw9498. doi:10.1126/science.aaw9498
Fraser, J., Williamson, I., Bickmore, W. A., and Dostie, J. (2015). An overview of genome organization and how we got there: from FISH to Hi-C. Microbiol. Mol. Biol. Rev. 79, 347–372. doi:10.1128/mmbr.00006-15
Fudenberg, G., and Imakaev, M. (2017). FISH-ing for captured contacts: towards reconciling FISH and 3C. Nat. Methods 14, 673–678. doi:10.1038/nmeth.4329
Fudenberg, G., Imakaev, M., Lu, C., Goloborodko, A., Abdennur, N., and Mirny, L. A. (2016). Formation of chromosomal domains by loop extrusion. Cell Rep. 15, 2038–2049. doi:10.1016/j.celrep.2016.04.085
Funke, J. J., Ketterer, P., Lieleg, C., Schunter, S., Korber, P., and Dietz, H. (2016). Uncovering the forces between nucleosomes using DNA origami. Sci. Adv. 2 (11), e1600974. doi:10.1126/sciadv.1600974
Garcia-Saez, I., Menoni, H., Boopathi, R., Shukla, M. S., Soueidan, L., Noirclerc-Savoye, M., et al. (2018). Structure of an H1-bound 6-nucleosome array reveals an untwisted two-start chromatin fiber conformation. Mol. Cell 72, 902–915. doi:10.1016/j.molcel.2018.09.027
Gibcus, J. H., Samejima, K., Goloborodko, A., Samejima, I., Naumova, N., Nuebler, J., et al. (2018). A pathway for mitotic chromosome formation. Science 359 (6376), eaao6135. doi:10.1126/science.aao6135
Goes, A. C. S., Cappellen, D., Santos, G. C., Pirozhkova, I., Lipinski, M., Vassetzky, Y., et al. (2011). Loop domain organization of the p53 locus in normal and breast cancer cells correlates with the transcriptional status of the TP53 and the neighboring genes. J. Cell. Biochem. 112, 2072–2081. doi:10.1002/jcb.23129
Goloborodko, A., Imakaev, M. V., Marko, J. F., and Mirny, L. (2016a). Compaction and segregation of sister chromatids via active loop extrusion. eLife 5, 1–16. doi:10.7554/elife.14864
Goloborodko, A., Marko, J. F., and Mirny, L. A. (2016b). Chromosome compaction by active loop extrusion. Biophys. J. 110, 2162–2168. doi:10.1016/j.bpj.2016.02.041
Gordon, F., Luger, K., and Hansen, J. C. (2005). The core histone N-terminal tail domains function independently and additively during salt-dependent oligomerization of nucleosomal arrays. J. Biol. Chem. 280, 33701–33706. doi:10.1074/jbc.M507048200
Grossfield, A. (2020). Data from: Wham: the weighted histogram analysis method. Version XXXX. http://membrane.urmc.rochester.edu/wordpress/?page_id=126.
Hart, K., Foloppe, N., Baker, C. M., Denning, E. J., Nilsson, L., and MacKerell, A. D. (2012). Optimization of the CHARMM additive force field for DNA: improved treatment of the BI/BII conformational equilibrium. J. Chem. Theor. Comput. 8, 348–362. doi:10.1021/ct200723y
Huang, J., Lopes, P. E. M., Roux, B., and MacKerell, A. D. (2014). Recent advances in polarizable force fields for macromolecules: microsecond simulations of proteins using the classical drude oscillator model. J. Phys. Chem. Lett. 5, 3144–3150. doi:10.1021/jz501315h
Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B. L., et al. (2017). Charmm36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods 14, 71–73. doi:10.1038/nmeth.4067
Ishida, H., and Kono, H. (2017). H4 tails potentially produce the diversity in the orientation of two nucleosomes. Biophys. J. 113, 978–990. doi:10.1016/j.bpj.2017.07.015
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi:10.1063/1.445869
Krietenstein, N., Abraham, S., Venev, S. V., Abdennur, N., Gibcus, J., Hsieh, T.-H. S., et al. (2020). Ultrastructural details of mammalian chromosome architecture. Mol. Cell 78, 554–565. doi:10.1016/j.molcel.2020.03.003
Lai, C.-L., Chen, C., Ou, S.-C., Prentiss, M., and Pettitt, B. M. (2020). Interactions between identical DNA double helices. Phys. Rev. E 101, 32414. doi:10.1103/PhysRevE.101.032414
Lemkul, J. A., and MacKerell, A. D. (2016). Balancing the interactions of mg2+ in aqueous solution and with nucleic acid moieties for a polarizable force field based on the classical drude oscillator model. J. Phys. Chem. B 120, 11436–11448. doi:10.1021/acs.jpcb.6b09262
Lequieu, J., Córdoba, A., Moller, J., and De Pablo, J. J. (2019). 1CPN: a coarse-grained multi-scale model of chromatin. J. Chem. Phys. 150, 215102. doi:10.1063/1.5092976
Li, M., Zhang, J. Z., and Xia, F. (2016). Constructing optimal coarse-grained sites of huge biomolecules by fluctuation maximization. J. Chem. Theor. Comput. 12, 2091–2100. doi:10.1021/acs.jctc.6b00016
Lieberman-aiden, E., van Berkum, N. L., Williams, L., Imakaev, M., Ragoczy, T., Telling, A., et al. (2009). Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 326, 289–293. doi:10.1126/science.1181369
Liu, M., Lu, Y., Yang, B., Chen, Y., Radda, J. S. D., Hu, M., et al. (2020). Multiplexed imaging of nucleome architectures in single cells of mammalian tissue. Nat. Commun. 11, 1–14. doi:10.1038/s41467-020-16732-5
Luger, K., Mäder, A. W., Richmond, R. K., Sargent, D. F., and Richmond, T. J. (1997). Crystal structure of the nucleosome core particle at 2.8 Å resolution. Nature 389, 251–260. doi:10.1038/38444
Luque, A., Collepardo-Guevara, R., Grigoryev, S., and Schlick, T. (2014). Dynamic condensation of linker histone C-terminal domain regulates chromatin structure. Nucleic Acids Res. 42, 7553–7560. doi:10.1093/nar/gku491
MacPherson, Q., Beltran, B., and Spakowitz, A. J. (2018). Bottom-up modeling of chromatin segregation due to epigenetic modifications. Proc. Natl. Acad. Sci. U.S.A. 115, 12739–12744. doi:10.1073/pnas.1812268115
Meaburn, K. J., Agunloye, O., Devine, M., Leshner, M., Roloff, G. W., True, L. D., et al. (2016). Tissue-of-origin-specific gene repositioning in breast and prostate cancer. Histochem. Cell Biol. 145, 433–446. doi:10.1007/s00418-015-1401-8
Milo, R., Jorgensen, P., Moran, U., Weber, G., and Springer, M. (2009). BioNumbers-the database of key numbers in molecular and cell biology. Nucleic Acids Res. 38, D750–D753. doi:10.1093/nar/gkp889
Mirny, L. A., Imakaev, M., and Abdennur, N. (2019). Two major mechanisms of chromosome organization. Curr. Opin. Cell Biol. 58, 142–152. doi:10.1016/j.ceb.2019.05.001
Miyamoto, S., and Kollman, P. A. (1992). Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 13, 952–962. doi:10.1002/jcc.540130805
Moller, J., Lequieu, J., and De Pablo, J. J. (2019). The free energy landscape of internucleosome interactions and its relation to chromatin fiber structure. ACS Cent. Sci. 5, 341–348. doi:10.1021/acscentsci.8b00836
Münkel, C., and Langowski, J. (1998). Chromosome structure predicted by a polymer model. Phys. Rev. E 57, 5888–5896. doi:10.1103/physreve.57.5888
Noid, W. G., Chu, J.-W., Ayton, G. S., Krishna, V., Izvekov, S., Voth, G. A., et al. (2008). The multiscale coarse-graining method. I. a rigorous bridge between atomistic and coarse-grained models. J. Chem. Phys. 128, 244114. doi:10.1063/1.2938860
Oliveira Junior, A. B., Contessoto, V. G., Mello, M. F., and Onuchic, J. N. (2020). A scalable computational approach for simulating complexes of multiple chromosomes. J. Mol. Biol. 433, 166700. doi:10.1016/j.jmb.2020.10.034
Philip, H., John, E. H., David, B. J., Rainer, K. S., and Lynn, R. H. (1993). Polymer models for interphase chromosomes. PNAS 90, 7854–7858. doi:10.1073/pnas.90.16.7854
Portela, A., and Esteller, M. (2010). Epigenetic modifications and human disease. Nat. Biotechnol. 28, 1057–1068. doi:10.1038/nbt.1685
Rao, S. S. P., Huntley, M. H., Durand, N. C., Stamenova, E. K., Bochkov, I. D., Robinson, J. T., et al. (2014). A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680. doi:10.1016/j.cell.2014.11.021
Sanborn, A. L., Rao, S. S. P., Huang, S.-C., Durand, N. C., Huntley, M. H., Jewett, A. I., et al. (2015). Chromatin extrusion explains key features of loop and domain formation in wild-type and engineered genomes. Proc. Natl. Acad. Sci. U.S.A. 112, E6456–E6465. doi:10.1073/pnas.1518552112
Sandholtz, S. H., MacPherson, Q., and Spakowitz, A. J. (2020). Physical modeling of the heritability and maintenance of epigenetic modifications. Proc. Natl. Acad. Sci. U.S.A. 117, 20423–20429. doi:10.1073/pnas.1920499117
Saurabh, S., Glaser, M. A., Lansac, Y., and Maiti, P. K. (2016). Atomistic simulation of stacked nucleosome core particles: tail bridging, the H4 tail, and effect of hydrophobic forces. J. Phys. Chem. B 120, 3048–3060. doi:10.1021/acs.jpcb.5b11863
Schwarzer, W., Abdennur, N., Goloborodko, A., Pekowska, A., Fudenberg, G., Loe-Mie, Y., et al. (2017). Two independent modes of chromatin organization revealed by cohesin removal. Nature 551, 51–56. doi:10.1038/nature24281
Stam, M., Tark-Dame, M., and Fransz, P. (2019). 3D genome organization: a role for phase separation and loop extrusion? Curr. Opin. Plant Biol. 48, 36–46. doi:10.1016/j.pbi.2019.03.008
Su, J.-H., Zheng, P., Kinrot, S. S., Bintu, B., and Zhuang, X. (2020). Genome-scale imaging of the 3D organization and transcriptional activity of chromatin. Cell 182, 1641–1659. doi:10.1016/j.cell.2020.07.032
Teif, V. B., and Clarkson, C. T. (2019). Nucleosome positioning. Ency. Bioinform. Comput. Biol. 1, 308–317. doi:10.1016/B978-0-12-809633-8.20242-2
Teif, V. B., Vainshtein, Y., Caudron-Herger, M., Mallm, J.-P., Marth, C., Höfer, T., et al. (2012). Genome-wide nucleosome positioning during embryonic stem cell development. Nat. Struct. Mol. Biol. 19, 1185–1192. doi:10.1038/nsmb.2419
Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., and Bussi, G. (2014). Plumed 2: new feathers for an old bird. Comput. Phys. Commun. 185, 604–613. doi:10.1016/j.cpc.2013.09.018
Walker, B., Jing, Z., and Ren, P. (2020). Molecular dynamics free energy simulations of atp:mg2+ and adp:mg2+ using the polarisable force field amoeba. Mol. Simul. 0, 1–10. doi:10.1080/08927022.2020.1725003
Wennberg, C. L., Murtola, T., Páll, S., Abraham, M. J., Hess, B., and Lindahl, E. (2015). Direct-space corrections enable fast and accurate lorentz-berthelot combination rule Lennard-Jones lattice summation. J. Chem. Theor. Comput. 11, 5737–5746. doi:10.1021/acs.jctc.5b00726
Zhang, B., and Wolynes, P. G. (2015). Topology, structures, and energy landscapes of human chromosomes. Proc. Natl. Acad. Sci. U.S.A. 112, 6062–6067. doi:10.1073/pnas.1506257112
Zhao, Z., Tavoosidana, G., Sjölinder, M., Göndör, A., Mariano, P., Wang, S., et al. (2006). Circular chromosome conformation capture (4C) uncovers extensive networks of epigenetically regulated intra and interchromosomal interactions. Nat. Genet. 38, 1341–1347. doi:10.1038/ng1891
Keywords: chromatin, nucleosome, molecular dynamics simulation, umbrella sampling, potential of mean force, coarse-grain
Citation: Zhang C and Huang J (2021) Interactions Between Nucleosomes: From Atomistic Simulation to Polymer Model. Front. Mol. Biosci. 8:624679. doi: 10.3389/fmolb.2021.624679
Received: 31 October 2020; Accepted: 09 February 2021;
Published: 12 April 2021.
Edited by:
Xiakun Chu, Stony Brook University, United StatesReviewed by:
Wenjun Xie, Massachusetts Institute of Technology, United StatesVinicius Contessoto, Rice University, United States
Copyright © 2021 Zhang and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jing Huang, huangjing@westlake.edu.cn