Nucleosome plasticity is a critical element of chromatin liquid–liquid phase separation and multivalent nucleosome interactions

1The Eukaryotic nucleus is a highly compartmentalized system that achieves its internal organization entirely without the use of membranes. Inside the nucleus, hundreds of millions of DNA base pairs are densely packed into a highly dynamic and heterogeneous structure known as chromatin. The basic building blocks of chromatin are nucleosomes: 10-nm wide nanoparticles composed of approximately 147 base pairs (bp) of DNA wrapped around a histone protein octamer (two copies each of H2A, H2B, H3, and H4),. To assemble chromatin, nucleosomes are first joined together by free DNA linker segments of varying lengths—measured in units of nucleosome repeat lengths (NRL = 147 bp + linker DNA length)— forming a “beads-on-a-string” structure termed the 10-nm fiber. Subsequently, using their charged and contoured surfaces, and charged and flexible protruding “arms” (histone tails), nucleosomes establish interactions with one another and with the DNA. These interactions trigger folding of the 10-nm fiber and the dense packing of DNA inside cells,.

2The structure of chromatin, beyond the 10-nm fiber, remains an intense topic of research and debate,–. The traditional textbook view suggests that 10-nm chromatin folds into a regular and rigid 30-nm zigzag fiber, where nucleosomes interact preferentially with their second-nearest neighbors (i.e., nucleosome i with i ± 2). While this zigzag fiber model is supported by in vitro studies of reconstituted chromatin arrays–, experiments interrogating chromatin inside cells have consistently failed to detect 30-nm fibers–. Instead, accumulating evidence is shifting the structural paradigm of chromatin in vivo in favor of the “liquid-like” or “fluid-like” model,, which proposes that 10-nm fibers condense into an irregular and dynamic polymorphic ensemble,,. The term “liquid” here is used to emphasize a structure that is absent of long-range translational order, where nucleosomes can flow and relax easily, and engage in interactions with a wide range of neighbors,,, like in a “sea of nucleosomes”. Nucleosomes within disordered chromatin are proposed to form heterogeneous groups of variable sizes and densities—interspersed with nucleosome-free regions. Furthermore, the density of nucleosomes, rather than the structure of the 30-nm fiber, is what seems to distinguish different chromatin regions (e.g., nucleosome density is higher in heterochromatin than in euchromatin),, and change upon differentiation. Sequencing-based methods that can resolve chromatin interactions in situ at single-nucleosome resolution—i.e., micrococcal nuclease chromosome conformation assay (Micro-C),, ionizing radiation-induced spatially correlated cleavage of DNA with sequencing (RICC-seq), and high-throughput chromosome conformation capture with nucleosome orientation (Hi-CO)—suggest that the irregular organization of chromatin is underpinned by dominant interactions among i and i ± 2 nucleosomes. Consistent models where chromatin exhibits a large-scale disordered organization but contains strong short-range zigzag contacts include the hierarchical looping model, and the multiplex higher-order folding model.

3The disordered behavior of nucleosomes is not surprising if one considers the notable intrinsic heterogeneity of chromatin physicochemical parameters in vivo (e.g., varying DNA sequences and epigenetic marks, heterogeneous distributions of post-translational histone modifications, non-uniform NRLs, presence of nucleosome-free regions, and dynamic nucleosome breathing and sliding motions). Many of these parameters can independently enable chromatin polymorphism, triggering the folding of 10-nm fibers into irregular loops, hairpins, and bends. Indeed, irregular nucleosome spacing,, nucleosome-free regions,,–, heterogeneous on/off dyad binding of linker histone proteins to the nucleosome,, low linker histone concentrations or subtype variations, inhomogeneous distributions of post-translational modifications, and the disordered nature of the linker histone protein can independently give rise to a plethora of nucleosome orientations and interactions. In concert, these factors can further amplify or control chromatin polymorphism,,.

4In the past three years, the paradigm of dynamic liquid-like behavior of nucleosomes within cells has gained significant traction due to the realization that chromatin and its associated multivalent biomolecules can undergo liquid–liquid phase separation (LLPS) in vitro and in cells–. LLPS is now postulated as a mechanism, alongside others, to explain genome compartmentalization without the use of physical membranes. Notable examples include the formation of constitutive heterochromatin via the transcriptional repressor HP1,,, as well as the emergence of important liquid-like structures such as the nucleoli, Cajal bodies, and nuclear speckles. In addition, the formation of super-enhancer regions has been recently associated with the phase separation of a combination of transcription factors and co-factors, chromatin regulators, non-coding RNAs, and RNA Polymerase II–,–.

5The origin of intranuclear phase separation is intricately linked to the complex and crowded biomolecular environment of the cell nucleus,; i.e., the nucleoplasm is a highly multi-component mixture of proteins and nucleic acids with varying compositions across different regions. Indeed, orientation-independent differential interference contrast microscopy reveals that the density of the total material in the nucleus is as high as 208 mg/ml within heterochromatin and 136 mg/ml in the surrounding euchromatin regions. The formation of diverse phase-separated chromatin compartments becomes thermodynamically stable in specific genomic regions when the concentration of key biomolecules, termed scaffolds, surpasses a threshold. These conditions allow scaffolds—normally dissolved in the nucleoplasm—to drive LLPS by minimizing their free energy through the formation of numerous attractive interactions with one another. Therefore, the features that affect binding among nucleosomes, and between nucleosomes and their chromatin-binding proteins (e.g., their chemical makeup and mechanical properties, along with the prevailing microenvironment) are expected to be crucial regulators of intranuclear LLPS. In particular, the capacity of biomolecules (e.g., proteins, RNA, DNA, and nucleosomes) to interact with at least three binding partners (i.e., multivalency) is essential for forming sufficient transient interconnections to compensate for the entropic loss due to the reduced number of microstates upon demixing,.

6Our work focuses on assessing an important feature of nucleosomes that likely pertains to chromatin LLPS: their inherent plasticity. That is, rather than static building blocks—as routinely considered in large-scale chromatin structural models—at short length-scales nucleosomes are highly dynamic and structurally irregular entities,–, and are better described as a “dynamic family of particles”. Nucleosomes in vivo can pack a broad range of DNA base pairs around the histone core (~100–170 bp), and can also have varied histone compositions and stoichiometries. In addition, thermal fluctuations cause some of the DNA–histone core interactions to spontaneously break and reform from one end of the nucleosome, while the majority of nucleosomal DNA remains wrapped around the histone core,—these phenomena, known as “nucleosome breathing” (also referred to as DNA breathing), is favored at physiological salt concentrations,. In vivo, the probability of nucleosome breathing across chromatin is likely very diverse too, as it can be sensibly altered by post-translational modifications of nucleosomes– and based on their constituent DNA sequences,,. Furthermore, the presence of H3K56Ac combined with DNA sequence changes can increase the rate of nucleosome unwrapping by one order of magnitude or more. The recruitment of chromatin-binding proteins can also affect the plasticity of nucleosomes. For instance, binding of multiple Swi6 molecules to a nucleosome has been shown to disrupt the DNA–histone bound state, increasing the exposure of the buried histone core to solvent, and crucially, promoting chromatin LLPS.

7Nucleosome structural fluctuations provide a transient opportunity for the binding of transcription factors to DNA, and hence, for transcription,. However, prior to the nucleosome barrier, pioneer and other transcription factors need to surpass the steric hindrance imposed by nucleosome–nucleosome interactions. This seems particularly challenging if we consider the apparent rigidity of nucleosomes within the 30-nm fiber structural models. Therefore, this begs a fundamental question that we target in this work: What are the physical and molecular factors that modulate the accessibility of nucleosomes within compact assemblies? To resolve this conundrum, assessing how the mesoscale properties of chromatin (e.g., density, flexibility, shape, and size) are impacted by the dynamic behavior of nucleosomes is crucial. In this regard, computational modeling of chromatin structure—both top-down polymer models trained on experimental datasets, and bottom-up mechanistic descriptions—offers an ideal complement to experiments. Polymer chromatin models trained on experimental datasets, are currently the best tools for investigating chromatin organization at the whole-nucleus scale; these approaches are essential to propose plausible structural ensembles of chromatin and test whether hypotheses on the structural behavior of chromatin and/or the physical mechanisms that dictate such structure are consistent with the experimental data. In partnership, mechanistic computational models can generate hypotheses and link them to the physicochemical properties of chromatin. A wide range of mechanistic coarse-grained models with nucleosome and sub-nucleosome resolution have been developed in the past few years,,,,,,– to bridge molecular and physicochemical information of nucleosomes to the mesoscale properties of chromatin. Future integration of both data-driven polymer models with mechanistic chromatin descriptions holds great potential for providing a complete view of chromatin organization.

8In the present study, we develop an advanced mechanistic multiscale chromatin model designed to investigate the connection between the fine molecular details of nucleosomes (including amino acid and DNA sequence, specific distributions of post-translational modifications and epigenetic marks, protein flexibility, DNA mechanical properties, and nucleosome plasticity) and the mesoscale (up to sub-Mb scale) organization of chromatin. By bridging atomistic and sub-Mb chromatin scales, our multiscale approach can dissect biophysical properties of chromatin that emerge from the dynamic formation and breakage of interactions among a few thousand nucleosomes with diverse molecular features; thereby, it can uncover the molecular and biophysical mechanisms that explain the self-organization and intrinsic LLPS of chromatin with diverse chemical makeups. When we model chromatin at residue/base-pair resolution, the dominant role of electrostatics, the high dimensionality of the system, and the resemblance of chromatin to a poly-branched polymer make attaining sufficient sampling highly non-trivial. To overcome these challenges, we develop, in tandem, a powerful Debye-length replica-exchange molecular dynamics (MD) simulation approach that is transferable and can be used to explore the energy landscapes of other intractable charged systems.

9Our multiscale model and Debye-length enhanced-sampling technique reveal that nucleosome breathing is promoted at physiological salt conditions, in agreement with experiments,. Whereas enhancement of nucleosome breathing at physiological salt concentrations drives chromatin to populate a highly dynamical, but compact, liquid-like structural ensemble,, inhibition of breathing at low salt gives rise to 30-nm zigzag fibers. This modulation of nucleosome breathing with salt, and its impact on chromatin self-assembly might help reconcile longstanding differences between fiber-based and in vivo chromatin models. Our simulations further explain that liquid-like chromatin organization is characterized by short-lived and orientationally diverse internucleosome interactions, which are mediated by transient non-specific DNA–histone tail contacts. In contrast, 30-nm fibers are sustained by long-lived regular face-to-face nucleosome interactions. Importantly, nucleosome plasticity promotes both liquid-like folding of individual chromatin systems and LLPS of chromatin arrays via the same physical mechanism: namely, it enhances the multivalency of nucleosomes and, therefore, the connectivity and stability of both compact chromatin and the phase-separated chromatin. The stochastic organization of nucleosomes within compact chromatin that we observe, both within single arrays and condensates, paints a much more permissive picture of nucleosome targeting than that offered by the fiber-based models. The realization that nucleosomes can be simultaneously stochastically organized and tightly packed might have important implications in reshaping the molecular mechanisms used to link chromatin structure to modulation of DNA accessibility.

10We have developed a multiscale model (Fig. 1 and Supplementary Fig. 1) that describes the physicochemical heterogeneity of chromatin and the plasticity of in vivo nucleosomes, as well as their impact on functionally relevant length scales (i.e., up to sub-Mb scales). This model integrates three complementary levels of resolution: atomistic representations of nucleosomes (Level 1), a chemically-specific chromatin model (Level 2), and a minimal chromatin model (Level 3). Integration of these spatiotemporal scales is necessary to probe the biophysical and molecular forces underpinning the modulation of large-scale chromatin organization by subtle chemical changes (e.g., in charge, hydrophobicity and flexibility) originating, for instance, from histone/DNA mutations, post-translational modifications of histones, DNA epigenetic marks, and binding of regulatory proteins.

11Our chemically-specific coarse-grained model features representations of breathing nucleosomes with all histone proteins resolved at the residue level (preserving the sequence-dependent charge, hydrophobicity, size, and flexibility of the atomistic histones). Double-stranded DNA is described at the base-pair step level (with charge and sequence-dependent mechanical properties described with a modified version of the Rigid Base Pair (RBP)– model with added phosphate charges; see “Methods” section). Parameters for our chemically-specific model are obtained from experimental amino-acid pairwise contact propensities–, large datasets of atomistic MD simulations of DNA strands, and bias-exchange metadynamics atomistic simulations of 211-bp nucleosomes.

12Notably, the physicochemical and molecular fidelity of histones and DNA within our chemically-specific coarse-grained model gives rise to a DNA polymer that spontaneously wraps ~1.7 times around the histone core, and that adopts the correct topology and left-handed chirality under weak negative DNA supercoiling, and the chiral inverted right-handed counterpart under weak positive supercoiling, consistent with experiments– (Supplementary Notes and Supplementary Fig. 8). The molecular resolution of our model also results in nucleosomes that naturally exhibit spontaneous breathing motions (i.e., without being primed to do so) and display force-induced unwrapping in quantitative agreement with experiments at single base-pair resolution (see “Model validation” section, Supplementary Notes, and Fig. 2).

13Further coarse-graining is essential to reduce the system dimensionality and investigate chromatin LLPS. Accordingly, from our chemically-specific coarse-grained simulations, we derive a consistent minimal chromatin model that describes each nucleosome with just a few particles. Such a dimensionality reduction enables the simulation of chemically heterogeneous sub-Mb scale chromatin regions and LLPS, while also considering nucleosome thermal fluctuations. Specifically, we use one ellipsoid for each histone core, and develop a “minimal RBP-like model” for DNA at a resolution of 5 bp per bead, with minimal helical parameters extracted from chemically-specific coarse-grained simulations of 200 bp DNA strands (see “Methods” section and the Supplementary Methods). Nucleosome–nucleosome and nucleosome–DNA interactions are modeled with orientationally dependent potentials fitted to reproduce internucleosome potentials of mean force calculated with our chemically-specific coarse-grained chromatin model (Supplementary Fig. 4 and Supplementary Methods). Significantly, a comparison of the rate of exponential decay of the autocorrelation functions of the chromatin radius of gyration in both models (Supplementary Fig. 7 and Supplementary Notes) indicates that the timescales in our minimal chromatin model are 10 times faster than those in our chemically-specific model.

14The detailed description of our models and the rationale behind their designs, resolutions, and parameters are given in the “Methods” section and the Supplementary Methods. A comparison of our model predictions with experiments is discussed in the “Model validation” section below and in the Supplementary Notes.

15We begin by comparing various quantities in our chemically-specific coarse-grained chromatin model with corresponding experimental observables. First, we find that the persistence length of DNA estimated in our simulations (i.e., performed using our modified RBP model with phosphate charges; see “Methods” section) agrees well with salt-dependent values from high-throughput tethered particle motion single-molecule and light scattering experiments, and with sequence-dependent measurements from cyclization assays (Supplementary Fig. 9 and Supplementary Notes). The residue-resolution protein model that we use to describe histones has been shown to reproduce well the experimental radii of gyration of several intrinsically disordered proteins.

16Next, we assess the accuracy of the critical dynamical unwrapping behavior of nucleosomes in our model by performing single-molecule force-extension simulations and comparing our results with those from force spectroscopy experiments–. Pulling a single nucleosome at the equilibrium speeds typically used in force spectroscopy experiments (e.g., ~0.1 mm/s for magnetic tweezers) is computationally unfeasible; that would require millisecond-long trajectories, which are presently not easily achievable. With the current computing power, force-extension steered MD would need to be performed at pulling speeds significantly above those required to maintain equilibrium conditions (e.g., 10–100 times faster than in experiments). To overcome this computational limitation, following the procedure of Lequieu et al., we perform equilibrium umbrella sampling to estimate the potential of mean force (PMF) of single nucleosome unwrapping, using the DNA end-to-end distance as the order parameter. Subsequently, we derive force-extension plots by taking the numerical derivative of the PMF curves (see force-extension curves in Supplementary Fig. 10).

17Consistent with experiments, the force-induced nucleosome unwrapping behavior predicted by our model can be separated into three equilibrium regimes (Fig. 2a, b), each spanning a force-extension region (Supplementary Fig. 10) that matches the experimental values in refs. – (see further discussion in Supplementary Notes). Most notably, by modeling nucleosomes with chemical and mechanical accuracy, we obtain quantitative agreement between our model predictions and the free-energy landscape for nucleosome unwrapping derived from mechanical unzipping experiments at single DNA base-pair resolution (Fig. 2c).

18Further, our model shows that unwrapping of the outer DNA turn is associated with a free-energy barrier (ΔG1 in Fig. 2) of ~11.5 kBT at 0.15 mol/L NaCl, which closely matches estimates from force spectroscopy experiments at similar salt conditions (~9–11.1 kBT),,. In agreement with magnetic tweezer experiments, we estimate an increase of ~10.5 kBT in the free-energy barrier for outer DNA turn unwrapping as the monovalent salt concentration decreases from 0.15 to 0.05 mol/L of NaCl (Fig. 2a, c). This increase in free energy can be explained from the enhanced electrostatic attraction between DNA and histones at low salt. Hence, as previously shown, nucleosome breathing is feasible at physiological salt concentrations but becomes increasingly challenging in low salt conditions.

19The free-energy landscape for nucleosome unwrapping computed with our model agrees well with several other experimental findings. For instance, in line with the dependency of nucleosome unwrapping on DNA sequence, nucleosomes with an unfavorable polyA DNA sequence have a free-energy barrier for unwrapping the outer DNA turn that is 1.5 kBT smaller than in 1KX5 nucleosomes. Histone tail clipping leads to the unwrapping of the outer turn with a negligible energetic penalty, in agreement with mechanical disruption experiments where unwrapping of the outer DNA turn occurs at near-zero forces after histone tail removal.

20Regarding chromatin behavior, our force-extension curves computed using an extension clamp (Supplementary Fig. 10) exhibit the typical saw-toothed pattern of optical tweezer experimental curves, where the force displays an abrupt drop accompanied by an increase in the extension due to the partial unwrapping of individual nucleosomes (see further discussion in the Supplementary Notes). Such behavior is also consistent with the step-like patterns emerging from magnetic tweezer experiments, where a force clamp is used instead. As we describe below, our sedimentation coefficients for 12-nucleosome 165-bp chromatin arrays are in quantitative agreement with the experimental values of Grigoryev and colleagues. In addition, the liquid-like chromatin ensembles that we report are in qualitative agreement with the disordered organization of nucleosomes within clutches derived from super-resolution microscopy experiments of chromatin in nuclei, and the chromEMT polymorphic nucleosome organization within chromatin in situ.

21Our chemically-specific coarse-grained chromatin model contains sufficient physical details to examine the effects of nucleosome breathing on the structure of small (<10 kb) chromatin systems. However, simulating chromatin arrays at high resolution (i.e., one bead per protein residue and DNA base pair) in physiological conditions is computationally challenging for various reasons. When represented at high resolution, chromatin is a high-dimensional system made up of a large number of oppositely charged particles (e.g., mainly lysine and DNA phosphate beads) that establish strong long-range interactions with each other. Such features give rise to a rugged energy landscape of chromatin that is populated by many competing low-lying minima separated by high energy barriers. A rugged energy landscape is difficult to sample with standard MD simulations, as transitions across the high energy barriers are rare within the accessible simulation timescales. Although Monte Carlo (MC) simulations are effective at overcoming high energy barriers, chromatin at high-resolution and under physiological conditions resembles a highly dense poly-branched polymer; consequently, most MC moves have a low acceptance probability due to steric clashes even after small displacements, just like in a dense liquid. To overcome these challenges and achieve sufficient sampling of chromatin at high resolution, we developed a Hamiltonian replica-exchange MD scheme that varies the Debye length across replicas. An advantage of this “Debye-length replica-exchange” approach, over standard temperature replica-exchange MD (T-REMD), is that it allows us to use a much smaller number of replicas (i.e., 16 instead of 80 for a 12-nucleosome system). This is achieved by implementing larger differences within the Hamiltonians of neighboring replicas but focused on selected degrees of freedom; such degrees of freedom are chosen to precisely modulate the electrostatic interactions, which contribute most strongly to the high energy barriers (Supplementary Methods). In addition, unlike the various temperature-replicas in T-REMD, each of the Debye-length-tuned replicas in our approach directly explores the behavior of chromatin at conditions easily accessible to experiments (e.g., each replica at a different salt concentration within the 0.01–0.15 mol/L NaCl range).

22Using this Debye-length replica-exchange approach to sample our chemically-specific model, we can compare the behavior of chromatin with nucleosomes that exhibit spontaneous DNA unwrapping (i.e., nucleosome breathing) versus cases where the nucleosomes are constrained to remain permanently wrapped (i.e., no nucleosome breathing is allowed). As a benchmark for our model, we focus on 12-nucleosome chromatin arrays with a regular NRL of 165 bp (or a linker length of 18 bp); since in vitro sedimentation coefficients are available for validation, and because short regular NRLs favor folding into ideal zigzag structures well-characterized by near-atomic resolution in vitro experiments,. Furthermore, using a regular NRL across the array allows us to exclude the structural heterogeneity stemming from linker DNA variability, and to focus on the effects of nucleosome thermal breathing.

23Our simulations reveal that constraining the nucleosomal DNA to remain fully and permanently attached to the histone core directs chromatin with short DNA linkers to fold into 30-nm rigid ladder-like zigzag fibers at 0.15 mol/L NaCl (see “Non-breathing” in Fig. 3a); i.e., where nucleosomes stack perfectly face-to-face with their second-nearest neighbors and rarely interact with other nucleosomes or in alternate orientations (top panel in Fig. 3b). The regular zigzag fiber structure we observe is analogous to that derived from cryo-EM experiments for 12-nucleosome chromatin systems with short linker lengths, and the 167-bp tetranucleosome crystals.

24More strikingly, our simulations show that the thermal breathing motion of nucleosomes destabilizes the formation of regular 30-nm zigzag fibers, and favors instead the organization of chromatin into a liquid-like ensemble (“Breathing” in Fig. 3a). This liquid-like ensemble encompasses a wide range of compact structures where nucleosomes interact with a multiplicity of neighbors in diverse orientations (i.e., face-to-face, side-to-side, and face-to-side) (bottom panel in Fig. 3b), as had been postulated by Maeshima and collaborators. These ensembles are consistent with the disordered organization of nucleosomes observed with super-resolution nanoscopy and chromEMT experiments of chromatin inside cells. Furthermore, despite the high degree of variation in internucleosome interactions, the torsional and bending rigidity of the DNA intrinsically directs nucleosomes to engage in frequent, but orientationally diverse, interactions with their second-nearest neighbors. Therefore, chromatin forms a structure that lacks long-range order but that is underpinned by dominant interactions among i and i ± 2 nucleosomes, as observed by RICC-seq, Micro-C,, and Hi-CO. As a result of nucleosome disorder, liquid-like chromatin manifests a higher degree of compaction and flexibility than the 30-nm fiber. These features are evident from the wider distributions of the sedimentation coefficients (shifted towards the right) of chromatin with breathing nucleosomes, compared to those of chromatin with non-breathing nucleosomes (Fig. 3c). Importantly, the sedimentation coefficient predicted by our model for 165-bp chromatin at physiological salt (0.15 mol/L NaCl) is in quantitative agreement with the experimental value. Our simulations also capture qualitatively the progressive decondensation of chromatin with decreasing monovalent salt observed experimentally.

25To rationalize the impact of nucleosome breathing in chromatin self-assembly, we quantify the average number of DNA base pairs that unwrap from each nucleosome due to thermal fluctuations as a function of the NaCl concentration (Fig. 3d). In vitro, single-molecule (sm) fluorescence resonance energy transfer (FRET) experiments show that nucleosome breathing is hindered at 0.02 mol/L NaCl (i.e., ~10% probability; 1.4 ms unwrapped versus 14 ms wrapped), but is promoted as the salt concentration is increased to 0.1 mol/L NaCl (i.e., ~ 30% probability; 1.5 ms unwrapped versus 3–4 ms wrapped). Consistently, we observe that the average number of unwrapped DNA base pairs increases significantly with salt (i.e., from 7 ± 2 bp at 0.01 mol/L NaCl to 22 ± 5 bp at 0.15 mol/L NaCl). Such enhanced unwrapping at physiological salt concentrations originates from the weakening of the DNA–histone core attraction as the electrostatic screening increases. Notably, enhanced unwrapping of nucleosomes at physiological salt implies that the DNA linker regions continuously lengthen and shorten (i.e., the average length of the linker DNA doubles from 18 to ~38 bp). Beyond simply increasing the fluctuations in the internucleosome distances, dynamic variations in linker DNA length give rise to a very important concomitant effect: diversification of the rotational angles among immediately linked nucleosomes. In other words, because the DNA is a helix that twists 360° every ~10.5–10.7 bp, when the length of the DNA between two nucleosomes changes by a few bases, the nucleosomes are not only spaced out differently but also variably rotated with respect to one another. It is precisely the emergence of this marked heterogeneity in internucleosome distances and rotational angles that explains, from a molecular point of view, the loss of long-range translational order in the organization of nucleosomes within liquid-like chromatin. Consistently, mesoscale simulations showed that large DNA linker variations increase the structural heterogeneity of chromatin, and experiments have shown that internucleosome rotational angle variability favors chromatin structural heterogeneity,. Moreover, mesoscale modeling revealed that binding of non-histone proteins, which can locally bend linker DNA, also destabilizes the regular 30-nm fiber folding. The structural effects of nucleosome breathing that we observe are more pronounced in chromatin arrays with short NRLs, as the internucleosome orientations within such systems are otherwise highly restricted by their short linkers. As the linker DNA lengthens, fluctuations become more energetically favorable and begin to intrinsically promote the heterogeneous nucleosome–nucleosome organization that sustains the liquid-like behavior of chromatin (Supplementary Fig. 11 and Supplementary Notes).

26The modulation of nucleosome breathing with salt and its impact on chromatin structure may help explain why ordered and disordered structural chromatin models have been derived from in vitro and in vivo data, respectively. Such differences have already been attributed to the low salt concentrations used in many of the in vitro experiments,, and to the regularity of reconstituted chromatin arrays—i.e., with strong nucleosome positioning sequences and homogeneous linker DNA sequences, uniform NRLs, homogeneous histone protein compositions, and a relatively small number of nucleosomes (~4–100). Our work further suggests that considering the disparities in the dynamic behavior of nucleosomes at physiological versus low salt is important in such debate.

27The ability of our chemically-specific model to resolve the motions of individual amino acids and DNA base pairs within compact chromatin enables us to examine the precise contributions of each of these species in directing chromatin organization. Specifically, we compute the fraction of time each amino acid or DNA base pair in a given nucleosome mediates internucleosome interactions. We categorize these interactions into three main groups: DNA-, globular histone-, and histone tail-mediated interactions. This analysis reveals that the molecular driving forces that stabilize the liquid-like organization of chromatin versus the regular 30-nm fiber folding are strikingly different (Fig. 4).

28Electron microscopy and single-molecule force spectroscopy experiments on H4-tail cross-linked chromatin, show that folding of 30-nm zigzag fibers is driven by face-to-face interactions between the H4-tail of one nucleosome and the H2A histone on the surface of another. In agreement, we observe that 30-nm fibers are sustained by interactions between globular histones on the surfaces of the two stacked nucleosomes (including those within the acidic patch), between the H4-tail and DNA, and more modestly between the H4-tail and the acidic patch (Fig. 4a, b). Besides the H4-tail—which has an ideal location on the nucleosome face—the other histone tails play a less prominent role in the folding of the 30-nm fiber. Our simulations further illuminate that the face-to-face stacking of nucleosomes with their second-nearest neighbors within zigzag fibers is long-lived (i.e., once stacked, nucleosomes rarely unstack; Fig. 5a).

29A fascinating insight stemming from our simulations is that the molecular interactions sustaining the compact state of liquid-like chromatin are instead highly heterogeneous (Fig. 4c, d) and transient (i.e., nucleosomes bind and unbind dynamically; Fig. 5b); being most strongly contributed by non-specific electrostatic interactions between the various disordered histone tails (via their lysines and arginines) and the DNA (both nucleosomal and linker DNA). Hence, unlike in the 30-nm fiber, within liquid-like chromatin the important acid patch region is free to recruit a wide range of chromatin-binding factors present in cells. Such diverse short-lived interactions are consistent with chromatin dynamically transitioning between the heterogeneous compact structures that make up the liquid-like ensemble. The key importance of histone tail–DNA electrostatic interactions is supported by experiments demonstrating that chromatin with tail-less nucleosomes fails to condense. In addition, mutation of all H4 arginine and lysine residues to alanine, or histone tail acetylation inhibits nucleosome–nucleosome interactions and the intrinsic LLPS of chromatin.

30Interestingly, despite the huge differences in the molecular interactions that stabilize the liquid-like and the fiber-like structures of chromatin, the dominant role of H4K16 is uncontested in both cases (Fig. 4). Among the contacts established by the H4-tail, the electrostatic interactions mediated by K16 are the strongest (most significantly with DNA in both types of chromatin structures, and more modestly with the acidic patch in the zigzag fibers). This dominance of the H4K16 residue is in agreement with the well-known decompaction triggered by H4K16 acetylation,, and the observation that reversible acetylation of H4K16—one of the most frequent post-translational modifications across organisms—has diverse functional implications.

31Besides characterizing the molecular features of liquid-like chromatin, we were interested in understanding the underlying physical principles that drive chromatin to adopt such a disordered organization, in terms of its thermodynamics. The structural heterogeneity of the liquid-like ensemble (Fig. 3c) indicates that such organization decreases the free energy of chromatin by expanding the number of accessible microstates (entropy gain), compared to those available in a regular 30-nm fiber. Next, to assess the variation in enthalpy, we define the “valency” of a nucleosome as the average number of nucleosome–nucleosome contacts it establishes (Supplementary Methods). We observe that nucleosomes within liquid-like chromatin have a significantly higher valency at physiological conditions than nucleosomes within 30-nm fibers (Fig. 5b). Hence, beyond the entropic driving force, nucleosomes within a liquid-like ensemble decrease their free energy by establishing numerous, but still transient, attractive interactions that maximize the enthalpy gain upon chromatin compaction. Numerous weak and transient internucleosome contacts are preferred over fewer strong longer-lived face-to-face interactions, as the latter represents a greater entropic cost.

32Recent groundbreaking experiments discovered that 12-nucleosome reconstituted chromatin undergoes intrinsic LLPS—i.e., without the aid of additional proteins—under physiological salt concentrations in vitro and when microinjected into cells. Extensive studies characterizing the LLPS of proteins and nucleic acids have demonstrated that multivalency is the dominant driving force for their LLPS–,–. That is, proteins with high valencies can stabilize LLPS by forming numerous weak attractive protein–protein,,, protein–RNA, and/or protein–DNA, interactions that compensate for the entropy loss upon demixing. Furthermore, binding of multiple Swi6 (the Schizosaccharomyces pombe HP1 protein) molecules to nucleosomes was recently shown to reshape the nucleosome in a manner consistent with nucleosome breathing—i.e., increasing the solvent exposure of core histones—and subsequently promote HP1-chromatin LLPS. These ideas, together with our observations of nucleosome valency enhancement from spontaneous breathing, led us to hypothesize that nucleosome plasticity could be crucial in facilitating the intrinsic LLPS of chromatin.

33To investigate this phenomenon and gain molecular and thermodynamic insight, we use our minimal coarse-grained chromatin model, as it simultaneously predicts chromatin ensembles in quantitative agreement with those of our chemically-specific model (Supplementary Fig. 4) and can simulate a solution of hundreds of interacting chromatin arrays. Specifically, we perform direct coexistence simulations of systems containing 125 independent 12-nucleosome chromatin arrays with a uniform NRL of 165 bp at different conditions. The direct coexistence method involves simulating two different phases—the condensed (chromatin-enriched) liquid in contact with the diluted (chromatin-depleted) liquid—in the same simulation box separated by an interface– (Supplementary Methods). From these simulations, we compute full liquid–liquid phase diagrams of chromatin at constant room temperature in the “NaCl concentration” versus “chromatin density” space (Supplementary Methods and Supplementary Fig. 5), and compare the results for chromatin with breathing nucleosomes and, as a control, for nucleosomes that are artificially constrained to remain fully wrapped (i.e., non-breathing).

34Our phase diagrams allow us to compare the conditions under which chromatin LLPS takes place spontaneously for the two representations. This analysis reveals that, besides promoting the liquid-like behavior of individual chromatin arrays, nucleosome breathing increases the range of stability of chromatin LLPS by enhancing the valency of nucleosomes (Fig. 6a). That is, chromatin spontaneously forms condensates above a critical monovalent salt concentration; i.e., where screening by counterions is sufficiently strong to eliminate the DNA–DNA repulsion and encourage the formation of numerous weak and transient attractive internucleosome interactions. Compared with breathing nucleosomes, when nucleosomes are constrained to remain fully wrapped, chromatin phase separation requires higher NaCl concentrations to become thermodynamically stable (Fig. 6a); the limited valency of non-breathing nucleosomes necessitates the formation of stronger nucleosome–nucleosome interactions to obtain sufficient enthalpic gain for LLPS. Beyond the critical solution salt concentration, the size of the liquid–liquid coexistence region is significantly larger for chromatin with nucleosomes that breathe spontaneously (Fig. 6a); i.e., under the same solution conditions, breathing nucleosomes yield a more dense and, hence, more stable condensed liquid.

35The physical forces governing LLPS of chromatin can be analyzed by computing the average density of molecular connections (i.e., bonds) that nucleosomes form per unit of volume within the condensed phase as a function of the salt concentration. Similar to its effect in enhancing the valency of nucleosomes within single chromatin arrays, we observe that nucleosome breathing promotes LLPS because it fosters a much higher nucleosome–nucleosome connectivity within the condensed liquid (Fig. 6b). A densely connected chromatin liquid phase fulfills two crucial requirements: making LLPS thermodynamically favorable and ensuring that the condensate mantains liquid-like properties. For instance, a condensed chromatin liquid and a gel would exhibit a similar percolating network structure (i.e., where each chromatin array is bound to at least one other array). However, while chromatin arrays diffuse freely within liquids, they are dynamically arrested within gels, as a result of increasingly strong and long-lived nucleosome–nucleosome interactions. Therefore, to preserve the liquid-like properties of the chromatin condensate, nucleosome–nucleosome interactions must remain sufficiently weak, and thus short-lived. Importantly, these weak nucleosome-nucleosome interactions must establish a densely connected liquid network to compensate for the entropy loss upon demixing, and make LLPS thermodynamically stable. Hence, liquid chromatin condensates are favored by the dynamic formation and rupture of a large number of weak attractive nucleosome–nucleosome interactions, which in turn are facilitated by the dynamic breathing of nucleosomes.

36We introduce a mechanistic multiscale model of chromatin designed to investigate the connection between the fine-atomistic details of nucleosomes and the emergence of chromatin self-organization and LLPS in systems with over a thousand nucleosomes. By combining atomistic simulations, a residue/base-pair resolution model, and a minimal representation of chromatin, our multiscale approach enables the study of collective effects of amino acid sequence and mutations, post-translational modifications, histone secondary structural changes, DNA sequence, and nucleosome dynamics in the modulation of the mesoscale structural properties of chromatin.

37Our simulations put forward nucleosome plasticity as a key driving force of the intrinsic liquid-like behavior of chromatin at physiological salt concentrations. In vivo, nucleosome plasticity can originate not only from nucleosome breathing, but also from nucleosome sliding, binding of proteins like HP1, nucleosome remodeling, post-translational histone modifications, histone replacement, and other mechanisms. We observe that nucleosome plasticity transforms nucleosomes from the uniform and static disc-like repeating units needed to sustain rigid zigzag fibers, to highly heterogeneous and dynamical particles that engage in promiscuous nucleosome–nucleosome interactions, sample a wide range of internucleosome rotational angles, and spontaneously self-assemble into disordered structures. Such a disordered organization of nucleosomes is in agreement with the liquid-like model of Maeshima and colleagues,,, and the ultrastructure of in vivo chromatin visualized by chromEMT and super-resolution nanoscopy.

38Regardless of the marked heterogeneity and promiscuity of internucleosome interactions, we find that nucleosomes within liquid-like chromatin engage in frequent i and i ± 2 contacts, as observed in sequencing-based experiments of chromatin in situ–. However, unlike those sustaining zigzag fibers, the dominant second-nearest neighbor interactions found within liquid-like chromatin are highly orientationally diverse and transient, and do not resemble a stack of tetranucleosomes. Furthermore, in contrast to the 30-nm zigzag fibers, which are mostly assembled via long-lived face-to-face H4-tail to acidic patch interactions, liquid-like chromatin is instead stabilized by a diverse range of short-lived non-specific DNA–histone tail electrostatic interactions. Therefore, within liquid-like chromatin, the nucleosome acidic patch region can be more easily accessed by the myriad of chromatin-binding factors that have been proposed to target it in vivo,, which may include HP1. Furthermore, controlled access to the acidic patch region has been hypothesized to play a crucial role in the modulation of chromatin remodeling motors that regulate nucleosome sliding.

39We also find that the attractive interactions that maintain the DNA wrapped around the histone core are predominantly electrostatic in nature. Hence, nucleosome breathing is boosted by electrostatic screening at physiological salt concentrations and hindered at lower salt concentrations, as also shown experimentally. Our work further reveals that significant nucleosome breathing at physiological salt favors the liquid-like behavior of chromatin even in artificially homogeneous oligonucleosomes (e.g., with uniform DNA linker lengths and DNA sequences and regular histone compositions). In contrast, lowering the salt concentration, as done in some in vitro experiments, inhibits nucleosome fluctuations and drives chromatin to form ordered 30-nm zigzag fibers. Accordingly, the modulation of nucleosome breathing with salt may aid in reconciling discrepancies between the ordered and disordered chromatin structural models derived from in vitro and in vivo experiments, respectively.

40Importantly, our work demonstrates that liquid-like chromatin is simultaneously stochastically organized and tightly packed. This suggests that chromatin compaction does not immediately imply steric barring of enzymes and that indeed chromatin as a compact liquid-like system is optimum for DNA searchability; i.e., liquid-like chromatin may provide easier and more homogeneous DNA access for processes like “scanning and targeting genomic DNA” without the need for chromatin to undergo decompaction, as postulated by Maeshima and collaborators.

41The strong link between the stochastic organization of chromatin—which we show might facilitate access to nucleosomal DNA—and nucleosome plasticity might have important functional implications. Nucleosomes represent a fluctuating barrier for the binding of transcription factors to DNA, and hence, for transcription,; the maximally repressive state is that of the fully wrapped nucleosome, while the maximally non-repressive state is that of a nucleosome-free region. However, although the majority of transcription factors seem to bind to nucleosome-free DNA regions,, spontaneous nucleosome breathing may provide transient access to some of these proteins to their nucleosome binding sites and, hence, facilitate transcription initiation–, especially at physiological ionic strengths. The increased access of chromatin due to nucleosome breathing might be most relevant for rationalizing the mechanisms of binding of pioneering factors, which are a special class of proteins that target nucleosomal DNA and in compact chromatin,. Additionally, in vitro, nucleosome breathing has been suggested to play a role in transcription elongation, by facilitating the movement of RNA polymerase II ternary elongation complex across the nucleosomal DNA,. In vivo, spontaneous nucleosome breathing, combined with the action of chromatin remodellers, allows the Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)-associated protein 9 (Cas9) to access the nucleosomal DNA. Our observations also support one potential mechanism influencing the control of gene transcription kinetics: spontaneous thermal breathing. This breathing behavior has been proposed to allow nucleosomes to adopt a wide range of promoter configurations, some of which transiently facilitate transcription and others that momentarily inhibit it; hence, modulating the rate of discontinuous transcription of genes that includes bursts of activity. Our work strongly suggests that enhancement of nucleosome breathing at promoters would indeed foster their structural heterogeneity and favour LLPS.

42Significantly, we demonstrate that the same molecular and biophysical driving forces that sustain the liquid-like behavior of nucleosomes within single compact chromatin arrays also promote the intrinsic phase separation of a solution of chromatin arrays. Specifically, nucleosome breathing fosters not only the formation of disordered and flexible chromatin ensembles but also LLPS and the emergence of phase-separated chromatin compartments. Interestingly, when nucleosome breathing is artificially suppressed (analogous to the behavior of strongly positioning artificial sequences), chromatin shows a reduced propensity for LLPS. Thus, intranuclear conditions that can spontaneously tune intrinsic nucleosome breathing via modulation of electrostatic interactions (e.g., changes in ionic salt concentration, pH, DNA/histone mutations) are expected to alter the spatial organization, degree of compaction, and compartmentalization of chromatin. Modulation of nucleosome breathing may, therefore, represent a key mechanism used by cells, not only to organize chromatin but also, to modulate its function.

43Our work also reveals that nucleosome breathing and multivalency are intricately linked and are positively correlated. Multivalency has been previously identified as an important property governing protein LLPS, both in the cytoplasm and nucleoplasm. Furthermore, protein multivalency and, therefore, LLPS can be modulated by several external (e.g., salt concentration, temperature, pH, multi-component composition) and intrinsic factors (e.g., protein mutations, post-translational modifications, and disorder-to-order transitions). For chromatin, we show that an increase in the monovalent salt concentration facilitates nucleosome breathing, and subsequently, enhances nucleosome multivalency. Moreover, our work strongly suggests that, in the absence of significant environmental changes, within a fiber-like chromatin model, the valency of nucleosomes can be considered as roughly static; i.e., for short linker DNAs, the torsional rigidity of the DNA locks second-nearest nucleosome neighbors in an arrangement where they bind exclusively to one another and rarely to other nucleosomes. In contrast, within a plastic nucleosomal framework, nucleosomes, like proteins, possess an inherent capacity (i.e., nucleosomal breathing) to dynamically modify their valency, and ultimately regulate their functionality.

44Recent landmark experiments by Sanulli et al. demonstrated that the binding of many HP1 molecules unexpectedly reshapes nucleosomes and makes the histone core more accessible to the solvent; consistent with nucleosomal DNA unwrapping, which in turn can increase the availability of the core for interactions that facilitate LLPS. Our work further suggests that the increase in the plasticity of nucleosomes induced by HP1 is consistent with the amplification of: (a) the local flexibility of chromatin, (b) the range of accessible nucleosome–nucleosome pair orientations, (c) the effective valency of nucleosomes, and (d) the transient nature of the inter-nucleosomal attractive interactions that lead to chromatin LLPS.

45Together, our work postulates that nucleosome plasticity is an intrinsic property of nucleosomes that facilitates chromatin stochastic self-assembly and LLPS, and hence, contributes to regulating the organization and membrane-less compartmentalization of the genome. These findings advance the molecular mechanisms and biophysical understanding of how the liquid-like organization of the genome is formed and sustained, and how it can be regulated. Modulation of nucleosome plasticity might have important implications in the functional organization of the genome and in the control of gene transcription parameters.

46In this work, we develop a bottom-up multiscale modeling technique—combining atomistic representations (Level 1) with two levels of coarse-graining (Levels 2 and 3)—to link the chemical heterogeneity of chromatin and the spontaneous breathing motions of nucleosomes to chromatin self-assembly and its LLPS. We have chosen a multiscale strategy because it allows us to take advantage of (1) the ability of atomistic models to elucidate how chemical changes transform the local behavior of proteins and DNA, describe the binding of proteins to chromatin, and establish how DNA sequence transforms its mechanical properties, and (2) the capacity of coarse-grained descriptions to reduce the system dimensionality markedly (e.g., represent a 300 kb chromatin region, or about 30 million atoms plus solvent, with as little as ~50,000 beads). The need for combining two levels of coarse-graining, instead of just one, stems from our interest in describing systems with thousands of nucleosomes (Level 3) while retaining essential physicochemical information mapped from the all-atom level (our intermediate resolution, or level 2). The two interconnected levels of coarse-grained models we have developed are described below along with the rationale behind their designs. An animated illustration of our multiscale strategy can be found in https://sef43.gitlab.io/.

47When designing our chemically-specific coarse-grained model (Level 2), our goal was to fulfill two opposing requirements. The first requirement was to describe proteins and DNA at high-enough resolution to capture the effects of DNA and amino-acid sequence variations in nucleosome–nucleosome interactions, and also the intrinsic nucleosome breathing motions that emerge naturally from the histone–DNA interactions and the mechanical properties of the DNA. The second requirement was to reduce the number of degrees of freedom within chromatin as much as possible to efficiently simulate oligonucleosome systems. Describing proteins at a resolution of one bead per amino acid, and the DNA at a resolution of one bead per base pair, although expensive computationally, is the ideal choice to realize these two goals; i.e., such a resolution is needed to retain the full chemical composition of a heterogeneous chromatin array and map, from the bottom up, the diverse physicochemical properties of the distinct amino acid and nucleobases that make up each of the different nucleosomes within chromatin. This resolution is also needed to consider amino acid point mutations, if desired, interrogate the role of specific histone residues in mediating chromatin organization, probe binding of additional proteins to specific chromatin residues, and describe the contributions of distinct amino acids to nucleosome unwrapping. Furthermore, while histone core proteins are largely α-helical and exhibit relatively small structural fluctuations in crystallographic studies, and in atomistic MD simulations, histone tails are largely disordered and highly flexible. A resolution of one bead per residue is the coarsest resolution that can describe the detailed topology of the histone globular domains and the flexibility of the histone tails, and consider disorder-to-order transitions that might be triggered by post-translational modifications. In terms of DNA, a resolution of one bead per base pair is the coarsest resolution that can capture in full the influence of DNA sequence on the mechanical properties of DNA (i.e., twist, roll, and twist). Collectively, our chemically-specific model reduces the number of particles in a 5 kb chromatin region, or in a ~25 nucleosome system, from ~0.5 million atoms (plus solvent) to only ~25,000 beads. Below, we provide more details of this model.

48Following the work of Dignon and colleagues, we represent each amino acid explicitly—both within the histone globular regions and histone tail regions— by using a single bead that carries the charge, hydrophobicity, and size of its atomistic counterpart. Each bead is defined as a point particle with an excluded volume centered on the Cα atom of the amino acid it represents. For each amino acid, the bead diameter is calculated from the experimentally measured van der Waals volumes and assuming that amino acids have a spherical shape, as done previously. We also assign a charge to each bead corresponding to the total charge of the related amino acid. The sequence-dependent hydrophobic attraction between specific amino acid pairs is accounted for by the Kim–Hummer model,, which consists of a shifted and truncated Lennard–Jones potential with parameters derived from experimental amino-acid pairwise contact propensities. Using our previous microsecond-long Bias-exchange metadynamics molecular dynamics simulations of 211-bp nucleosomes, we differentiate between the globular histone core, residues that retain their secondary structure throughout the simulation, and the disordered histone tails. We then treat amino acids within intrinsically disordered regions as fully flexible polymers (i.e., with no energetic penalty for bending) using a harmonic potential with a stiffness bond constant kb of 20 kcal/mol/Å2 and a resting length r0 of 3.5 Å (Supplementary Methods), as proposed by Dignon and colleagues. We describe the relatively small structural fluctuations within the histone globular domains by building an elastic network model, which avoids the need of including internal non-bonded terms in these regions. In practice, we take a representative structure from the highest populated cluster in our atomistic simulations as the reference structure and connect all the globular histone core beads that are within 7.5 Å of each other with harmonic springs and form a Gaussian elastic network model (GNM). For the harmonic bond interaction of the GNM, we use a spring constant kGNM of 20 kcal/mol/Å2, and the equilibrium distance among amino-acid pairs is set equal to its value in the reference atomistic structure. Linker histones, and other chromatin-binding proteins of interest, can be described with our model in the exact same way as histones. That is, first, an initial coarse-grained model is built starting from a high-resolution structure of the protein (e.g., from ref. and adding the intrinsically disordered regions as described in ref. for H1). Next, using either experimental structural data or atomistic simulations, the degree of order/disorder in the protein regions is defined. Finally, each amino acid is represented with one bead, and residues within globular domains are connected by an elastic network model, while those within disordered regions are described as fully flexible chains.

49The second crucial feature of our model is that it considers sequence-dependent DNA mechanical properties by using the rigid base pair (RBP)– model with added phosphate charges. The RBP model represents one DNA base-pair step with one single ellipsoid and depicts the DNA conformational changes in terms of harmonic deformations of six helical parameters (three angles: twist, roll, and tilt, and three distances: slide, shift, and rise) that account for the relative orientations and positions of neighboring base-pair planes. The DNA mechanical potential energy is computed from the sum of harmonic distortions of equilibrium base-pair step geometries. We used the Orozco group parameters, computed from MD atomistic simulations—i.e., the equilibrium values by fitting Gaussian functions to the distributions of helical parameters, and the elastic force constants by inversion of the covariance matrix in helical space. In practice, our model represents single base pairs, both within the nucleosomal and linker DNA, by one coarse-grained bead defined by a position vector, r, and an orientation quaternion,. We add two virtual charge sites to each DNA ellipsoid (i.e., one per phosphate approximating the shape of the DNA phosphate backbone) to consider the crucial electrostatic interactions that drive chromatin self-organization. We implement this RBP plus charged virtual sites in LAMMPS (http://lammps.sandia.gov) with an ellipsoid defined by two-point particles with a negligible but non zero mass. While the combined three-particle base-pair bead is treated as a single rigid body for the dynamics, the individual components each contribute to the calculation of the potential and forces.

50Besides the excluded volume and hydrophobic non-bonded interactions, we consider electrostatic interactions between the charged beads by means of a Debye–Hückel potential, but omit all non-bonded interactions between directly bonded beads. The binding of the nucleosomal DNA to the histone protein core is achieved through these protein–DNA electrostatic and hydrophobic interactions, resulting in a nucleosomal DNA that wraps ~1.7 times around the histone core and exhibits spontaneous unwrapping/re-wrapping. The force-induced unwrapping behavior of these nucleosomes is in quantitative agreement with experiments at single base-pair resolution (see “Model validation” section and the Supplementary Notes). To model nucleosomes that are artificially constrained to be “non-breathing”, we describe the histone core together with the bound nucleosomal DNA as a single GNM, using the same 7.5 Å threshold and bond parameters. This results in the nucleosomal DNA being constrained to remain permanently bound to the histone core, and inhibits both nucleosome breathing and sliding.

51The Debye–Hückel approximation is a mean-field theory where the ion density is given by the linearized Boltzmann distribution and effects like ion condensation, ion correlations, ion heterogeneity, and specific ion binding are ignored. By invoking such approximation, we assume that the effects of monovalent counterions in the solution can be reduced to simply screening the mean electrostatic potential from the chromatin charges, and as such, we describe the decay of charge–charge interactions with distance by a Yukawa function. Such an approximation is crucial to reduce the high dimensionality of the chromatin phase space and to enable us to sample it efficiently. Importantly, this approximation is exact in the low salt limit, and previous chromatin coarse-grained models have shown that it captures well the salt-dependent compaction of chromatin at the low to moderate monovalent salt concentrations present inside cells (≤0.15 mol/L of NaCl),. However, its numerical implementation becomes progressively challenging as the salt concentration reaches very low values because estimating the interactions requires larger and larger Debye lengths and, hence, larger and larger cutoffs (Supplementary Methods).

52An additional limitation of invoking the Debye–Hückel approximation is that it neglects the effects of Mg2+ ions in solution, which are highly abundant inside cells (i.e., ~10–20 mmol/L but with the fraction of free Mg2+ ions estimated at <5% of this). Mg2+ ions are thought to outcompete monovalent ions and interact preferentially with the multiple closely-spaced negative charges on the nucleosome surface and on the DNA. Reassuringly, at such low concentrations, Mg2+ ions have been shown to only modestly enhance chromatin compaction in vitro, e.g., by ~2% and ~13% when comparing the sedimentation coefficients of a 165-bp 12-nucleosome array in 0.15 mol/L NaCl versus in 1 and 2 mmol/L MgCl2, respectively. Importantly, landmark work of the late Jonathan Widom showed that in the presence of 0.1 mol/L Na+, the addition of 0–10 mmol/L Mg2+ did not substantially change the equilibrium constant for nucleosome unwrapping. More recent sm-FRET studies observed that increasing the ionic strength of the buffer, by raising instead the NaCl concentration to up to 0.4 mol/L, had a minimal impact on nucleosome unwrapping. Hence, the effects that increased ion screening by Mg2+ ions at physiological concentrations could exert in the modulation of nucleosome unwrapping are likely subtle and expected to take place at much finer scales than what we set to capture with our chemically-specific model. Furthermore, we note that the experimental data that we have used to quantitatively validate the force-induced unwrapping response of our model (i.e., mechanical unzipping experiments of Wang and colleagues) is for nucleosomes already immersed in a buffer containing both 0.1 mol/L NaCl and 0.5 mmol/L MgCl2, suggesting that the balance of parameters in our model captures well this regime.

53In summary, our chemically-specific coarse-grained model of chromatin preserves the atomistic shape and size of the nucleosome core, the length and flexibility of the histone tails, the explicit charges and hydrophobic nature of all the different amino acids within the histone protein (including those in the important acid patch region), the sequence-dependent mechanical properties of the DNA, and the thermal breathing motions of nucleosomes. A list of all the parameters we use in our chemically-specific model, along with the choice of values and justification is given in Supplementary Table 9.

54In a similar fashion to the design of our chemically-specific model, our goal when developing the minimal model (Level 3) was to find the coarsest possible representation of nucleosomes that simultaneously allows us to (1) account for the breathing motions of nucleosomes and quantitatively reproduce the chromatin structural ensembles that we observe with our chemically-specific model, and (2) model chromatin systems with thousands of nucleosomes to reach functionally relevant scales or phase-separating systems. To map the breathing motions of nucleosomes from our chemically-specific simulations, and also the torsional and bending properties of the linker DNA (which we find crucially determine the configurational ensembles of chromatin), we require an explicit representation of the DNA—this is indeed one of the limiting factors defining the resolution of our minimal model. The exact resolution of our minimal DNA model is dictated by the physical diameter of a canonical B-DNA strand and the type of potentials we use to describe interactions among minimal DNA beads. Specifically, because the excluded volume of a DNA bead is determined by the physical diameter of B-DNA (which is ~20 Å), and we use the repulsive part of a Lennard–Jones potential to model DNA–DNA repulsion (see below), the coarsest resolution that a DNA bead can have in our model is 5 bp (as this corresponds to ~17 Å in length). Such a resolution guarantees that the distance among sequentially bonded DNA beads is smaller than their excluded volume size, and thus, makes the DNA a self-avoiding polymer. A higher resolution would work too, but would introduce unnecessary computational expense, and would reduce the size of the chromatin systems we could study. The resolution of histone proteins, in contrast, can be reduced significantly more. We choose a resolution of one bead per histone octamer, as this is sufficient to represent the geometry of the nucleosome, the topological distribution of DNA beads around the histone core, and the chemically-specific strength of the nucleosome–nucleosome interactions. Together, the resolution of our minimal model allows us to reduce the number of particles in a 300 kb chromatin region, or 1500 nucleosomes, from ~30 million atoms (plus solvent) to only ~50,000 beads.

55When devising our minimal model, we took care in capturing the manner in which nucleosomes interact with one another and the mechanical properties of the DNA (torsional and bending flexibility). Both features are crucial in defining the conformational space of nucleosomes. To account for these two features, while enabling efficient sampling of a few thousand nucleosomes, our minimal chromatin model approximates the shape and size of the histone core with a single bead (using an ellipsoid of 28 × 28 × 20 Å radii), and represents linker and nucleosomal DNA with one finite-size orientable sphere (using an ellipsoid of 12 × 12 × 12 Å radii) for every 5 base-pair steps. Inspired by the success of the RBP model in adequately approximating the atomistic mechanical properties of DNA at a single base-pair-step resolution, we propose a “minimal RBP-like model” for DNA at a resolution of 5 base pairs per bead. Accordingly, we optimize the minimal helical parameters from chemically-specific coarse-grained simulations of 200 bp DNA strands (Supplementary Methods). Then we describe nucleosome–nucleosome and nucleosome–DNA interactions by a series of orientationally dependent potentials fitted to reproduce internucleosome potentials of mean force calculated with our chemically-specific coarse-grained chromatin model (Supplementary Fig. 4 and Supplementary Methods). To represent breathing nucleosomes, we analyze our chemically-specific chromatin simulations to determine the fraction of DNA that remains predominantly bound to the histone core and define only that fraction as nucleosomal DNA. For each snapshot in the chemically-specific model trajectories, we will get a different definition which builds up a set of structures representative of the thermodynamic equilibrium distribution of nucleosome breathing states. By using diverse structures from this set, we can incorporate the effects of nucleosome breathing in the model without having to directly simulate completely free DNA, significantly reducing the degrees of freedom, enabling superior sampling. In the case of the non-breathing nucleosomes, we simply define the nucleosomal DNA following the same definition used in the chemically-specific model. Then, in both cases, breathing and non-breathing nucleosomes, we attach permanently the nucleosomal DNA to the histone core bead. Hence breathing and non-breathing nucleosomes are constructed using equilibrium configurations from different chemically-specific simulations (Supplementary Fig. 3). In addition, to account for the slightly exposed histone core in breathing nucleosomes, we add an additional anisotropic potential to the total energy. We parameterize this anisotropic term to be consistent with experimental force spectroscopy experiments on single nucleosomes– (Supplementary Methods and Supplementary Figs. 2, 4). A list of all the parameters we use in this minimal chromatin model, along with the choice of values and justification is given in Supplementary Table 10.

56To drive chromatin systems over the free-energy barriers and achieve complete sampling, we developed a Hamiltonian exchange method that attempts exchanges between replicas with different values of the Debye-length of the screened-coulomb interaction (λD) while keeping the temperatures the same. As in standard Hamiltonian exchange, the exchange probability is given by:where xi are the chromatin coordinates of the ith replica, and the potential energy function at Debye length —the original Debye length of the ith replica. The exchange is accepted or rejected based on the Metropolis criteria, and upon exchange, the potential energy functions (or coordinates) are switched. For our 12-nucleosome chromatin systems, we find that at the value of λD = 8.0 Å (corresponding to 0.15 mol/L salt), the chromatin structures are compact and suffer from sampling issues, increasing the Debye length to 15 Å gives rise to open structures that can sample effectively. We find that a range of 8.0–15.0 Å requires 16 replicas to get an acceptance probability of 30 %. This is significantly less than the ~80 replicas that we found would be required for standard T-REMD to sample the range of 300–600 K with a similar exchange probability. Additionally, the different Debye lengths are all at physically relevant salt concentrations and 300 K, therefore while increasing the sampling we can investigate the salt-dependent behavior of the system. Therefore, without the need for reweighing, the information extracted directly from each of the different replicas can be compared with experimental observables at different salt concentrations, whereas for T-REMD only the replica at 300 K (or the experimental temperature) is directly (ie., without reweighing) relevant for experimental comparison.

57Simulations were performed using LAMMPS (version 3rd March 2020) with our custom code (see Code availability). We used the program 3DNA (version 2.3) within our model building methods. All data analysis was done using Python (version 3.8.5) with NumPy (version 1.19.2) and SciPy (version 1.5.2). All data was plotted using Matplotlib (version 3.3.2). Images were rendered using the Open Visualization Tool (OVITO) software (version 3.0.0). We used the Weighted Histogram Analysis Method (WHAM) program, (version 2.0.9) to calculate PMFs.

58Further information on research design is available in the Nature Research Reporting Summary linked to this article.