Tension-Dependent Free Energies of Nucleosome Unwrapping.

Nucleosomes form the basic unit of compaction within eukaryotic genomes, and their locations represent an important, yet poorly understood, mechanism of genetic regulation. Quantifying the strength of interactions within the nucleosome is a central problem in biophysics and is critical to understanding how nucleosome positions influence gene expression. By comparing to single-molecule experiments, we demonstrate that a coarse-grained molecular model of the nucleosome can reproduce key aspects of nucleosome unwrapping. Using detailed simulations of DNA and histone proteins, we calculate the tension-dependent free energy surface corresponding to the unwrapping process. The model reproduces quantitatively the forces required to unwrap the nucleosome and reveals the role played by electrostatic interactions during this process. We then demonstrate that histone modifications and DNA sequence can have significant effects on the energies of nucleosome formation. Most notably, we show that histone tails contribute asymmetrically to the stability of the outer and inner turn of nucleosomal DNA and that depending on which histone tails are modified, the tension-dependent response is modulated differently.

Introduction

Eukaryotic genomes are packaged into a compact, yet dynamic, structure known as chromatin. The basic building block of chromatin is the nucleosome, a disk-like structure consisting of 147 base pairs of DNA wrapped into 1.7 superhelical turns around proteins known as histones.[1,2] These histone proteins form what is known as the histone octamer, a protein complex consisting of two copies of histone proteins H2A, H2B, H3, and H4. The surface of the histone octamer is highly positive, which interacts favorably with the negative backbone of DNA. As a result, at sufficiently low ionic conditions, nucleosomes are stable and spontaneously form.

The locations of nucleosomes along the genome play a central role in eukaryotic regulation. DNA segments incorporated into nucleosomes are inaccessible to other DNA binding proteins, including transcription factors and polymerases, and thus nucleosomes must be disrupted in order for the cellular machinery to access nucleosomal DNA. As such, the positions occupied by nucleosomes provide an additional, important mechanism by which eukaryotic genomes are regulated. Indeed, past work has demonstrated that deregulation of these processes is implicated in numerous diseases, including cancer.[3−5] Quantifying the strength of interactions within the nucleosome structure and the forces required to disrupt them is of fundamental importance to understanding the delicate dynamics of chromatin compaction.

Optical-trapping single-molecule techniques have been particularly effective at probing the many interactions within the nucleosome. In these experiments, chromatin fibers[6−9] or single nucleosomes[10−15] are subjected to pico-newton scale forces, thereby providing the ability to precisely perturb the native nucleosome structure. By analyzing the deformations that result from these forces, one can infer the underlying strength of binding energies within the nucleosome. Following the initial work by Mihardja et al.,[10] a consensus is emerging[11,13,15] in which a single nucleosome is disrupted in two stages. In the first, at 3 pN, the outer wrap of DNA is removed from the histone surface. This first wrap is removed gradually and is considered to be an equilibrium process, where spontaneous unwrapping and rewrapping events can be observed under a constant force. The second transition occurs at forces 8–9 pN and occurs rapidly via so-called “rips”, where the remaining wrap of DNA is suddenly released.[10] More recently, these transitions have been shown to depend on torque (i.e., DNA supercoiling via twist),[13] and to occur asymmetrically due to variability in the bound DNA sequence.[15]

Several theoretical and computational studies have sought to help interpret these experimental results. Following the initial work of Kulić and Schiessel,[16] most current treatments represent the nucleosome as an oriented spool, and the unbound DNA as a semiflexible worm-like chain.[17,18] While earlier studies were only able to detect a single distinct unwrapping transition,[16,19] consistent with the first experiments,[6] more recent work[17,18] has been able to reproduce the two transitions observed by Mihardja et al.[10] By relying on simple, primarily analytic models, these studies have provided significant insights into the fundamental physics that govern interactions within the nucleosome. Such approaches, however, have necessarily had to invoke assumptions and introduce adjustable parameters in order to describe experiments[17,18] (e.g., the DNA–histone binding energy). This limits their ability to predict nucleosomal behavior under different conditions, such as variations in DNA sequence or ionic environment, without resorting to additional experimental data. Additionally, these models cannot explicitly account for histone modifications, which are central to nucleosome positioning and higher-order chromatin structure.[9,20−23]

A complementary approach, which should in principle enable prediction of nucleosomal interactions under a wide array of situations, could rely on molecular models where the nucleosome can be assembled or disassembled explicitly. Though these approaches are particularly promising, their success has been frustrated by the inability to access the experimentally relevant time scales of stretching, typically rates of 100 nm/s. Clearly, these time scales are inaccessible to atomistic simulations, yet even a highly coarse-grained spool-like model of the nucleosome still employed stretching rates several orders of magnitude too fast.[19] There is therefore a need to develop models and methodologies to facilitate more direct comparisons between optical-trapping experiments and molecular-level calculations. If successful, such models could reveal the subtle, yet incredibly important effects of DNA-sequence and histone modifications on nucleosome stability.

In this work, we build on a recently proposed coarse-grained model of the nucleosome[24−26] to examine its response to external perturbations. A computational framework is proposed in which the tension-dependent response of the nucleosome is examined at equilibrium, thereby providing access to the free energy of nucleosome unwrapping under tension. Our results are found to be in agreement with experimental measurements by Mihardja et al.[10] and Brower-Toland et al.,[9] and serve to demonstrate that it is indeed possible to reproduce the absolute binding free energies of nucleosome formation in terms of purely molecular-level information, without resorting to additional parameters. Importantly, that model is then used to predict the impact of DNA sequence and histone modifications on the relative free energies of binding within the nucleosome.

Results

A schematic representation of our simulation setup is shown in Figure a. As with optical-trapping experiments, the “state” of the nucleosome is represented by two parameters: the tension (or force) exerted on the DNA molecule, τ, and the extension of the DNA ends, r. To facilitate comparison with experiments,[10] the ends of the DNA are not torsionally constrained. Figure b shows instantaneous configurations of the nucleosome model for five different values of extension, r. Consistent with previous observations, the outer wrap of the nucleosome is first removed (A → T1 → B), followed by the inner wrap (B → T2 → C).

Figure 1

Model of nucleosome unwrapping. (a) Coarse-grained topology of nucleosome. DNA is represented by 3SPN2.C,[24] and the histone proteins by AICG.[26] Both the end-to-end extension, r, and tension, τ, are constrained during a simulation. (b) Unwrapping process. During extension, the wraps of DNA around histone proteins are removed one by one. T1 and T2 denote the transition states separating the first (A ↔ B) and second (B ↔ C) unwrapping events. Figures were generated using VMD.[27]

In order to quantify these transitions, we examine the tension-dependent free energy of nucleosome unwrapping. By calculating the tension-dependent free energy instead of a simple force–extension curve, as done previously,[19,28] we can determine the true equilibrium behavior of the unwrapping process. Additionally, by performing simulations at equilibrium, we circumvent the issue of time scales that frustrate comparisons of traditional nonequilibrium molecular simulations to optical pulling experiments.

A representative two-dimensional tension-extension free energy surface for the 601 positioning sequence[29] is shown in Figure . Rather than increasing linearly with tension, the extension is quantized into three well-defined vertical bands, located at ∼120, 420, and 700, corresponding to states “A”, “B”, and “C” in Figure . At low tension (τ < 3 pN), a low extension (r < 200) is preferred. As tension is increased (τ ≈ 4–8 pN), the minimum free energy shifts to intermediate values of extension (r ≈ 420). At higher tension (τ > 8 pN), the minimum free energy shifts to larger values of extension (r ≈ 700). The free energy penalty of low tension and high extension (e.g., τ = 3, r = 700) or high tension with low extension (e.g., τ = 12, r = 200) results in large energy barriers >40 kT.

Figure 2

Tension-dependent free energy surface of nucleosome unwrapping for 601 positioning sequence. The free energy surface demonstrates minima at extensions of r ≈ 120, 420, 700, depending on tension. As tension increases, the minimum-energy extension shifts to larger values. Consistent with Mihardja et al.,[10] two transitions are observed.

The tension-dependent transition can also be visualized by plotting one-dimensional “slices” of the free energy surface at different values of tension (Figure a). Visualizing the data in this way clearly demonstrates that there are three basins of nucleosome extension: “fully wrapped”, “partially wrapped” and “unwrapped”. The basin that is favored depends on the tension applied to the DNA ends. As tension increases, the free energy minima shifts first from the “fully wrapped” to the “partially wrapped” basin, and then to the “unwrapped” basin. The boundaries of these basins are defined by the locations of the transition states (i.e., local maximum in the free energy) that separate neighboring basins. The transition states separating the A → B transition, T1, and the B → C transition, T2, are shown in Figure b.

Figure 3

(a) Free energy versus extension for different values of tension with the 601 positioning sequence. The locations of the transition states, T1 and T2, are used to define three basins: “fully wrapped”, “partially wrapped”, and “unwrapped”. L0 represents the contour length of the DNA molecule. (b) Probability of observing the nucleosome in each free energy basin for different tensions. The “fully” and “partially” wrapped states are at equilibrium (i.e., equal probability) when τ1* = 3.2 pN. The “partially” and “unwrapped” states are at equilibrium when τ2* = 8.9 pN. Error bars represent standard deviation across four independent simulations.

Once these three basins are defined, we can determine the precise tension at which the outer and inner DNA turns unwrap from the nucleosome. This is obtained by converting the tension-dependent free energy into probabilities and then integrating these probabilities to determine the total probability of finding the system in each basin (see Methods). The corresponding results are shown in Figure b; it can be appreciated that the probability of finding the system in the “fully wrapped” or “partially wrapped” basin is equivalent when τ ≈ 3.2 pN. Thus, when τ ≈ 3.2 pN the outer turn of nucleosomal DNA is in equilibrium (in a statistical mechanics sense) with its unbound state. We define this tension as τ1*. Similarly, the probability of the nucleosome in the “partially wrapped” and “unwrapped” basins is the same (i.e., the inner wrap is in equilibrium) when τ ≈ 8.9 pN, defined as τ2). These values are in quantitative agreement with those measured by Mihardja et al.,[10] who observed that the outer and inner DNA loops were removed at 3 pN and 8–9 pN, respectively.

A complementary approach to estimate τ1* and τ2* is to determine the tension at which the free energy barriers of the forward and reverse reactions are equal.[17]Figure shows the corresponding tension-dependent free energy barriers of the outer (A ↔ T1 ↔ B) and inner unwrapping (B ↔ T2 ↔ C) events. At low tension, the energy barriers for the forward reactions, A → B and B → C, dominate, and the forward (i.e., unwrapping) reaction rate is low. As tension increases, the energy barriers for the forward reactions decrease, while those for the reverse increase, thereby causing the unwrapping reaction to proceed at a higher rate. When the energy barriers of the forward and reverse reactions are equal, the two basins are at equilibrium (in a transition state theory sense), and τ1* and τ2* can be determined. These unwrapping tensions are estimated to be 3.3 pN and 8.5 pN, in excellent agreement with the probability-based analysis of Figure b.

Figure 4

Free energy barrier heights of nucleosome unwrapping for 601 positioning sequence. Solid lines represent the unwrapping (forward) reactions, dotted lines represent wrapping (reverse) reactions. When the unwrapping and wrapping barriers are equal, the two basins are at equilibrium with one another. This is found when τ1* = 3.3 pN for the outer wrap, and τ2* = 8.5 pN for the inner wrap. ΔA†(τ1*) = 4 kT and ΔA†(τ2*) = 16 kT. Error bars represent standard deviation across four independent simulations.

The magnitude of the free energy barrier also helps explain the observation by Mihardja et al.[10] that the outer turn of DNA can be removed reversibly, while the inner turn cannot. Since the energy barrier separating the “Fully” and “Partially” wrapped states is only ≈5 kT, the system can quickly transition between states when held at τ = τ1*. In contrast, the “partial wrap” and “unwrapped” states are separated by an energy barrier of ∼18 kT, indicating that even at equilibrium the P ↔ U transition occurs slowly. Thus, removal of the outer wrap may appear to be reversible on the time scales of a typical optical trapping experiment, while the inner wrap may not. Further, because force–extension curves are usually obtained via optical trapping by pulling a nucleosome at a fixed velocity, the experiments may not observe a P → U transition until τ > τ2*. This would cause the experiments to overestimate the value of τ2* and lead to a sudden, irreversible “ripping” event. We also note that the barrier estimates in this work (ΔA1† = 4 kT, ΔA2† = 16 kT) are in excellent agreement with those predicted by Sudhanshu et al.[17] (ΔA1† ≈ 6 kT, ΔA2† ≈ 15 kT).

Electrostatics, Sequence Dependence, Histone Modifications

Having validated the proposed model against experimental data,[10] we now examine the influence of ionic environment, DNA sequence, and histone modifications on the stability of the nucleosome. Such variations can have a significant impact on nucleosome formation, and the precise molecular origins of their impact is still poorly understood.

We first investigate the origins of the tension-dependent response by exploring the role of DNA–DNA electrostatic repulsion on the stability of the nucleosome structure. Past theoretical work[16,18] has suggested that DNA–DNA repulsion within the nucleosome is central to its tension-dependent response. Other studies, however, have observed that DNA–DNA electrostatic repulsion is unimportant and that the correct response can be achieved by accounting for the tension-dependent orientation of the free DNA ends.[17] Since our proposed model explicitly includes both contributions, we can directly evaluate the importance of DNA–DNA repulsion on nucleosome unwrapping. To examine this effect, we disable DNA–DNA electrostatic repulsion in our model between base pairs separated by more than 20 base pairs. Only disabling electrostatics between distant regions of DNA was necessary to avoid implicitly lowering the persistence length of DNA by neglecting Coulombic interactions between neighboring base pairs. All electrostatics responsible for DNA-histone affinity, however, remain intact.

Our results are summarized in Figure a,b. As anticipated,[16] removal of DNA–DNA repulsive interactions has a greater impact on the outer DNA loop (Δτ1* = +1.7pN) than on the inner DNA loop (Δτ2* = +1.1pN). However, in the absence of DNA–DNA repulsions, the qualitative features of the tension-dependent response remain unchanged. These results indicate that while DNA–DNA repulsions play a role in nucleosome disassembly, they are not primarily responsible for the two unwrapping steps observed in experiments. Our results are also consistent with prior experimental measurements, where the role of DNA–DNA repulsion on the stability of the outer turn was observed to be small.[30]

Figure 5

(a) Schematic representation of proposed model with DNA–DNA repulsion removed. (b) Tension-dependent free energy barriers for 601 positioning sequence with DNA-DNA repulsion removed. Δτ1* and Δτ2* represent change relative to complete model. Error bars represent standard deviation across three independent simulations.

We next examine the impact of DNA sequence on the relative binding free energies of nucleosome formation. Optical trapping experiments could in principle be used to probe the sequence-dependent energies within the nucleosome, but recent literature studies have been limited to the 601 positioning sequence[10,11,13] and slight variations.[15] Instead, competitive reconstitution assays are the dominant experimental technique for characterization of sequence-dependent relative binding free energies.[31,32] To compare model predictions to these experiments, we use the technique employed by Freeman et al.,[25] where the relative binding free energies of different DNA sequences are assessed computationally using alchemical transformations and thermodynamic integration (see Methods). A comparison of predicted and experimental free energies, shown in Figure , indicates that, as with previous work,[25] the model adopted here accurately reproduces the binding free energies of many different sequences. In general, the key predictor of binding free energy is the sequence-dependent shape of the DNA molecule (i.e., minor groove widths and intrinsic curvature). Sequences that bind strongly (low ΔΔA) possess periodic sequence motifs (e.g., TA base steps) that impart a shape that favorably “fits” underlying histone structure.[32] In contrast, weakly binding sequences (large ΔΔA) do not possess these periodic motifs.

Figure 6

Sequence-dependent binding free energies. Squares denote model proposed by Freeman et al.[25] (obtained at 300 K and 150 mM ionic strength). Circles denote model proposed in this work, obtained at 277 K and vanishing ionic strength (for consistency with ref (31)). Despite differing solution conditions and DNA–protein interactions, both models reproduce the relative binding free energies of nucleosome formation. The DNA sequences used here are given in ref (25). The dotted line corresponds to quantitative agreement.

In addition to DNA sequence, the modification of histone tails is widely considered to be the single most important determinant of chromatin structure.[20] Methylated and acetylated histones are enriched at promoters of highly expressed genes and are thought to play a role in the strong positioning of certain nucleosomes.[21,22] Histone tails are central to nucleosome–nucleosome interactions, and their modification has important implications on chromatin’s three-dimensional structure.[23] Experiments[9] have also established that removal of histone tails has a significant impact on the stability of the nucleosome.

To examine the role of histone tails on nucleosome stability at a molecular level, we return to our earlier analysis and calculate the tension-dependent free energy of nucleosome unwrapping. Our results can be compared to the optical trapping experiments of Brower-Toland et al.,[9] where arrays of 17 nucleosomes were disassembled for different histone tail modifications, including the complete removal via trypsin digest or post-translationally via acetylation. In the model, we perform this trypsin digest in silico to each histone (see Figure a) and calculate the resulting tension-dependent response. Figure b shows the change in the equilibrium tension of the outer, Δτ1*, and inner, Δτ2*, turn of DNA due to the removal of different histone tails. The experimental measurements are also included and correspond to the impact of histone modifications on the stability of the inner turn of DNA (i.e., Δτ2*). Our results are in excellent agreement with experimental measurements and predict the effect of each histone modification to within ±0.5 pN. Yet our results provide additional insight into these experiments, whose limited spatial resolution prevented the observation of the individual release of the outer DNA turn. Most importantly, we observe that tails of different histones contribute asymmetrically to the stability of each turn of nucleosomal DNA. The H3/4 tails dominate the stability of the outer DNA turn, Δτ1* = −1.2 pN, but contribute weakly to the stability of the inner turn, Δτ2* = −0.2 pN. In contrast, H2A/B tails have a small effect on the outer turn, Δτ1* = 0.4 pN, but a significant effect on the inner DNA turn, Δτ2* = 2.5 pN. Therefore, depending on whether histone modifications occur on H2A/B or H3/4, the stability of the nucleosome will be modulated differently. This observation suggests a potent mechanism by which each turn of nucleosomal DNA can be independently regulated and could help to explain the importance and role of H3/4 modifications relative to those of H2A/B.

Figure 7

(a) Molecular configuration highlighting histone tails removed by in silico trypsin digest. The exact residues removed are given in the original work by Brower-Toland et al.[9] (b) Change in equilibrium tension of outer, Δτ1*, and inner DNA turn, Δτ2*, resulting from removal of H3/4 tails (gH3/4), H2A/B tails (gH2A/B), and all histone tails (gAll). Tensions are reported relative to Δτ1,0* and Δτ2,0*, the values reported previously in Figure for the 601 positioning sequence. Experimental data correspond to removal of the inner turn of DNA (i.e., Δτ2*). Error bars represent standard deviation across three independent simulations or reported in ref (9).

Conclusion

In this work we have demonstrated that a molecular-model of the nucleosome, composed of two coarse-grained models of DNA and proteins,[24−26] can be combined parameter-free to accurately reproduce the tension-dependent response of nucleosome unwrapping. This model quantitatively reproduces the unwrapping forces observed in experiments[9,10] and the barrier heights predicted by prior theoretical studies.[17] We then demonstrated that this model can be used to examine, without adjustment, the role of subtle phenomena in nucleosome formation such as DNA–DNA Coulombic repulsion, DNA-sequence, and histone tail modifications.

Our proposed approach opens up a new avenue for theoretical examinations of nucleosome stability. As a first step, this model can aid the interpretation of recent optical pulling experiments where the nucleosome is subjected to torque[13] and is suitable for examining subtle features within the nucleosome such as sequence-dependent asymmetric unwrapping.[15] Further, analysis of experimental measurements can become increasingly sophisticated, because our model provides a tool for interrogating raw data, including the fluctuations, from FRET and optical pulling experiments. This combination of experiment and simulation could help to resolve nanometer-scale phenomena such as dynamic DNA–protein contacts and could effectively increase the spatial resolution of experimental measurements to the base-pair level. Yet the potential of our approach extends beyond single nucleosome experiments and can begin to elucidate many unsolved questions within chromatin biophysics. How does the methylation of specific histone tails (and not others) enhance the positioning of certain nucleosomes? What are the free energies of different folded chromatin structures, and how do histone modifications effect this energy landscape? What is the role of DNA sequence on these processes, and do certain DNA sequences dispose chromatin to different “folds”? The approach presented in this work represents an important step toward answering these questions.

Methods

The model adopted in this work relies on a coarse-grained model of DNA[24,25] and proteins,[26] which are combined to represent the nucleosome. Both models were developed independently, but they are implemented at the same level of description, thereby facilitating their concerted use. Specifically, for DNA we use the 3SPN coarse-grained representation, where each nucleotide is described by three force sites located at the phosphate, sugar and base.[24,33−35] For the histone proteins, we use the “atomistic-interaction based coarse-grained model” (AICG), where the protein is represented by one site per amino acid located at the center of mass of the side chain.[26]

Interactions between the 3SPN2.C and AICG models included electrostatic and excluded volume effects. Phosphate sites with 3SPN were assigned a charge of −0.6 as described previously.[35] Each protein site was given the net charge of that residue at physiological pH (i.e., +1 for Arg, His, and Lys; −1 for Asp and Glu, 0 for others). As with prior work,[25] the effective charge of interactions between DNA and protein sites was scaled by a factor of 1.67 to bring the local charge of the phosphates back to −1. We note that DNA–protein interactions in this work differ slightly from those employed by Freeman et al.,[25] where, in addition to electrostatics, a small Lennard–Jones attraction was added between all DNA and protein sites. The strength of this attraction was very weak (ϵPro–DNA = 0.25 kJ/mol) and was originally included to reduce fluctuations within the nucleosome structure. Here we demonstrate that this weak interaction is unnecessary; by omitting it, both the relative and absolute formation free energies of the nucleosome can be reproduced. The combined model is effectively parameter-free: both the model of DNA and protein are included as originally proposed without any additional terms. Electrostatic forces are approximated by Debye–Hückel theory. Debye–Hückel theory invokes many assumptions about the electrostatic environment, and is not strictly valid for the highly charged association of the histone proteins and DNA. Nonetheless, Debye–Hückel theory provides a useful first-order approximation of Coulombic forces and is employed here, without resorting to higher-order techniques. All simulations were performed in the canonical ensemble using a Langevin thermostat and 150 mM ionic strength.

As an initial condition, we combine the 1KX5 crystal structure[2] of the nucleosome core particle with a proposed configuration of exiting DNA[36,37] to form a 223 base-pair structure, with 147 base pairs bound to the histone proteins and 38 flanking bases on each side. When using the 601 positioning sequence,[29] the flanking bases were chosen as polyA. This configuration was only used as the initial configuration, and no information from either structure was directly encoded into the nucleosome model.

To extract the tension-dependent free energy surface, two constraints were applied to the nucleosomal model. First, a constant force (i.e., tension) was applied to each end of DNA in order to mimic the experimental setup of optical-trapping experiments. Then, harmonic constraints were applied to the end-to-end extension of the DNA molecule, and umbrella sampling was performed to determine the free-energy as a function of DNA extension. In umbrella sampling, many independent simulations are performed at specified locations in phase space, and molecular fluctuations are used to estimate the local free energy surface at that location. These many local estimates are then systematically combined to obtain the total free energy surface.[38,39] Because tension is held at a constant value during a simulation, the resulting free energy “surfaces” are not truly continuous functions of tension. They are instead a compilation of two-dimensional “curves” that are plotted cocurrently to construct the “surface” presented in Figure .

The relative free energy of binding for different DNA sequences (ΔΔA) was calculated as described in detail previously.[25] Briefly, a thermodynamic cycle was defined that represents the relative sequence-dependent free energy of nucleosome formation, ΔΔA, as the difference between the free energy difference of two DNA sequences in the bulk, ΔAbulk, and bound to the histone proteins, ΔAbound (i.e., ΔΔA = ΔAbulk – ΔAbound); ΔAbulk and ΔAbound are determined by thermodynamic integration. The DNA sequences analyzed are given explicitly in ref (25).