help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wedemann, G.
Right arrow Articles by Langowski, J.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Wedemann, G.
Right arrow Articles by Langowski, J.

Biophys J, June 2002, p. 2847-2859, Vol. 82, No. 6

Computer Simulation of the 30-Nanometer Chromatin Fiber

Gero Wedemann and Jörg Langowski

German Cancer Research Center (DKFZ), Division Biophysics of Macromolecules (H0500), Im Neuenheimer Feld 280, 69120 Heidelberg, Germany


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

A new Monte Carlo model for the structure of chromatin is presented here. Based on our previous work on superhelical DNA and polynucleosomes, it reintegrates aspects of the "solenoid" and the "zig-zag" models. The DNA is modeled as a flexible elastic polymer chain, consisting of segments connected by elastic bending, torsional, and stretching springs. The electrostatic interaction between the DNA segments is described by the Debye-Hückel approximation. Nucleosome core particles are represented by oblate ellipsoids; their interaction potential has been parameterized by a comparison with data from liquid crystals of nucleosome solutions. DNA and chromatosomes are linked either at the surface of the chromatosome or through a rigid nucleosome stem. Equilibrium ensembles of 100-nucleosome chains at physiological ionic strength were generated by a Metropolis-Monte Carlo algorithm. For a DNA linked at the nucleosome stem and a nucleosome repeat of 200 bp, the simulated fiber diameter of 32 nm and the mass density of 6.1 nucleosomes per 11 nm fiber length are in excellent agreement with experimental values from the literature. The experimental value of the inclination of DNA and nucleosomes to the fiber axis could also be reproduced. Whereas the linker DNA connects chromatosomes on opposite sides of the fiber, the overall packing of the nucleosomes leads to a helical aspect of the structure. The persistence length of the simulated fibers is 265 nm. For more random fibers where the tilt angles between two nucleosomes are chosen according to a Gaussian distribution along the fiber, the persistence length decreases to 30 nm with increasing width of the distribution, whereas the other observable parameters such as the mass density remain unchanged. Polynucleosomes with repeat lengths of 212 bp also form fibers with the expected experimental properties. Systems with larger repeat length form fibers, but the mass density is significantly lower than the measured value. The theoretical characteristics of a fiber with a repeat length of 192 bp where DNA and nucleosomes are connected at the core particle are in agreement with the experimental values. Systems without a stem and a repeat length of 217 bp do not form fibers.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

DNA of most eukaryotic organisms is complexed with proteins into highly compact structures designated as chromatin (van Holde, 1989). DNA and histone protein assemble into nucleosomes; the nucleosome cores are connected by the linker DNA. At low salt concentration, chromatin adopts a beads-on-a-string-like appearance (Olins and Olins, 1974; Kornberg, 1974). At higher salt concentrations, chromatin compacts into a fiber of ~30-nm diameter (Finch and Klug, 1976). The consensus value for the linear mass density of this fiber is 6 to 7 nucleosomes/11 nm (Gerchman and Ramakrishnan, 1987; van Holde, 1989; and references therein). Assuming nucleosome repeat lengths of 180 to 220 bp, this corresponds to a DNA density of 1080 to 1320 bp/11 nm, or a compaction ratio of 35 to 43 compared with linear B-DNA. Models of the chromatin fiber have to take into account this mass density and the fiber diameter.

Many models for the structure of the 30-nm chromatin fiber have been developed (van Holde, 1989). Early analysis of x-ray diffraction, electron microscopy, and solution scattering data led to the development of the solenoid model (Finch and Klug, 1976). In this model, the nucleosomes are arranged in a regular helix with a period of six nucleosomes per turn and a pitch of 11 nm. The linker DNA follows a continuous spiral path between nucleosomes with a correspondingly large bend. Because of the tight bending of the linker DNA, the solenoid model requires a substantial input of elastic energy. Whereas this energy could be obtained, e.g., from an association of the linker DNA with the positively charged histone tails, models that do not require this energy input might be an interesting alternative.

In particular, such models have been developed recently based on new data from cryoelectron microscopy and scanning force microscopy-measurements and a reanalysis of the original data. That work resulted in the proposal of structural models that are compatible with the existing experimental data but lack the regularity of the solenoid model (van Holde and Zlatanova, 1995; Woodcock and Horowitz, 1995). Based on the so-called cross-linker models where the linker DNA forms a straight path between successive nucleosomes (Bordas et al., 1986; Kubista et al., 1990), Woodcock et al. (1993) proposed the zig-zag model. Here the nucleosomes form a more or less random zig-zag structure that can be mathematically described by the length of the linker DNA, the distance between the entry and exit point of the DNA emerging from chromatosomes, the entry-exit angle of the linker DNA (alpha ), and the tilt-angle between connected nucleosomes (beta ).

More support for the zig-zag model comes from a recent analysis of radiation-induced in vivo fragmentation of chromatin (Rydberg et al., 1998). Here a characteristic length distribution of the DNA fragments was found that could be correlated with model predictions of different structural arrangements of the chromatin fiber. The best agreement was found for the zig-zag model, whereas a solenoidal arrangement of nucleosomes gave fragment distributions that deviated considerably from their experiment.

There are, however, some questions that are not yet resolved by the zig-zag model. It describes correctly some qualitative and quantitative aspects of the fiber, such as its general aspect in microscopic images and its diameter; however, other experimentally known quantities like the linear mass density are not reproduced very well. Also, the inherent randomness of a large macromolecular structure like the chromatin fiber could only be taken into account by introducing random variations in geometrical parameters, such as nucleosome-nucleosome bending and twisting angles. The problem with all cited models is that they only describe the static, geometrical structure whereas effects of thermal fluctuations due to Brownian motions are neglected. Moreover, long-range forces such as electrostatic interactions are not contained in these "static" models. It would therefore be desirable to develop a description of the chromatin fiber that gives a correct picture of its accessible conformations at thermodynamic equilibrium and at the same time agrees with the known structural aspects. Such a description requires a numerical simulation.

Modeling the chromatin fiber on a computer is a challenging task because the molecule is very large. The use of atom-scale molecular dynamics is out of the question because the system is far too complex: a chromatin chain of 50 nucleosomes already contains approximately one million atoms, not counting the solvent molecules that must be accounted for in the simulation. Therefore, any model of the chromatin chain will necessarily be "coarse-grained," i.e., the basic units are much bigger than a single atom.

Earlier, a Brownian dynamics model that included DNA elasticity, electrostatic interactions between linker DNA segments, and thermal motions, but no internucleosome interactions, was developed in our group (Ehrlich et al., 1997). We could reproduce the hydrodynamic properties of polynucleosomes, but because of the missing internucleosome interactions and the approximation of nucleosomes as spheres, the model could not predict the correct mass density of the 30-nm fiber at physiological ionic strength.

Based on our earlier work we developed a new computer model that is presented here. The nucleosome cores, which in the crystal structure have an approximate flat disk shape, are modeled as ellipsoids. Their interaction is described using the Gay-Berne potential, which is a generalization of the Lennard-Jones potential for objects of ellipsoidal symmetry (Gay and Berne, 1981). This potential describes attractive as well as repulsive contributions of the internucleosomal interaction. The linker DNA is modeled as a chain of segments. DNA segments and the chromatosomes are coupled to each other through harmonic stretching, bending, and torsion potentials. The electrostatic repulsion between the DNA segments is modeled by line charges that interact via a Debye-Hückel potential. The parameters of all potentials are derived from experimental data. Typical configurations of the fiber at a given temperature are sampled using the classical Metropolis-Monte Carlo algorithm (Metropolis et al., 1953). Our model reproduces quantitatively the experimental data of the 30-nm fiber at physiological salt conditions such as diameter, mass density, tilt of nucleosomes and DNA to the fiber axis, and persistence length of bending.

Another important result from our simulations is that a variation of the twist angle between successive nucleosomes (beta ) has a stronger influence on the overall structure of the fiber than varying the strength of the internucleosomal interaction. The shape of the fiber can therefore be regulated much more effectively by changing geometrical parameters than through the internucleosomal interaction. Furthermore, we observed that for longer nucleosome repeats the stem motif is essential for the formation of the fiber.


    METHODS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Discretization of the chromatin chain

DNA that is not bound to the chromatosome is described by segments. They are chosen much smaller than the persistence length, so they are approximately straight: for linker DNAs 10 bp or shorter, one segment was used to connect the chromatosomes, and for longer linker DNAs, the length was divided into two segments. The position of each segment is determined by a position vector <A><AC>r</AC><AC>&cjs1164;</AC></A>i and a local coordinate system (fi, ui, vi) describing its orientation in space. The coordinate vector <A><AC>u</AC><AC>&cjs1164;</AC></A>i corresponds to the segment vector <A><AC>u</AC><AC>&cjs1164;</AC></A>i <A><AC>r</AC><AC>&cjs1164;</AC></A>i+1 - <A><AC>r</AC><AC>&cjs1164;</AC></A>i. The segment length is bi =  ui .

The position of a chromatosome is described by a vector to its center of mass and a local coordinate system where <A><AC>u</AC><AC>&cjs1164;</AC></A>i points from the center of mass to the geometrical center of the connecting points of the linker DNA, and <A><AC>v</AC><AC>&cjs1164;</AC></A>i is parallel to the axis of the nucleosome core.

Linker DNA is coupled to the chromatosome at two distinct points at a distance Dentry-exit. A chromatosome therefore forms a short segment by itself. The distance of the two DNA coupling points to the symmetry axis of the cylinder is Daxis-dna. If the linker DNA is coupled directly to the chromatosome, the value of Daxis-dna is equal to the radius of the nucleosome. In the case of a nucleosome stem this radius is larger. Fig. 1 shows how in both cases DNA segments and chromatosomes are connected.



View larger version (21K):
[in this window]
[in a new window]
 
FIGURE 1   Sketch of the geometrical structure of the model. Both types of coupling of DNA to the chromatosomes are shown. Coupling at distinct points (1) without and (2) with stem. DNA segments are denoted by straight lines, the gray rectangles symbolize the chromatosomes.

Elastic energies

The potentials of the elastic interactions are assumed to be harmonic. The strength constants of the interactions are named alpha <UP><SUB>(Y)</SUB><SUP>(X)</SUP></UP> in which X denotes the type of interaction, s denotes the stretching, b denotes the bending, t denotes the torsion, and Y denotes the type of interacting partners, DNA and Nuc.

The stretching potential of one segment is computed as
U<SUB><UP>stretch</UP></SUB>=&agr;<SUP>(<UP>s</UP>)</SUP><SUB><UP>Y</UP></SUB><UP>/</UP>b<SUP>0</SUP><SUB><UP>i</UP></SUB>(b<SUB><UP>i</UP></SUB><UP>−b</UP><SUP><UP>0</UP></SUP><SUB><UP>i</UP></SUB>)<SUP><UP>2</UP></SUP> (1)
The bending angle is computed between the segment vector <A><AC>u</AC><AC>&cjs1164;</AC></A>i+1 and its equilibrium position Bi relative to the connected segment. Bi is defined by the angles theta i and phi i:
<A><AC>B</AC><AC>&cjs1164;</AC></A><SUB>i</SUB>=<A><AC>f</AC><AC>&cjs1164;</AC></A><SUB>i</SUB><UP> sin</UP> &thgr;<SUB>i</SUB> <UP>cos</UP> &phgr;<SUB>i</SUB>+<A><AC>v</AC><AC>&cjs1164;</AC></A><SUB>i</SUB> <UP>sin</UP> &thgr;<SUB>i</SUB> <UP>sin</UP> &phgr;<SUB>i</SUB>+<A><AC>u</AC><AC>&cjs1164;</AC></A><SUB>i</SUB> <UP>cos</UP> &thgr;<SUB>i</SUB> (2)
In the case of no intrinsic bending theta i = phi i = 0 thus Bi = <A><AC>u</AC><AC>&cjs1164;</AC></A>i. The bending potential is then
U<SUB><UP>bend</UP></SUB>=&agr;<SUP>(b)</SUP><SUB>Y</SUB>/b<SUP>0</SUP><SUB>i</SUB>&thgr;<SUP>2</SUP><SUB>i</SUB>  <UP>with cos</UP>(&thgr;<SUB>i</SUB>)=<A><AC>B</AC><AC>&cjs1164;</AC></A><SUB>i</SUB><A><AC>u</AC><AC>&cjs1164;</AC></A><SUB>i+1</SUB> (3)
The coordinate systems of two connected segments i and i + 1 can be mapped on each other by an Euler-transformation with the angles (alpha i, beta i, gamma i) in which the rotation of the first angle is around <A><AC>u</AC><AC>&cjs1164;</AC></A>i. Then alpha i + gamma i - tau i is the total torsion. tau i is the intrinsic torsion. The torsion potential is then
U<SUB><UP>tors</UP></SUB>=&agr;<SUP>(t)</SUP><SUB>Y</SUB>/b<SUP>0</SUP><SUB>i</SUB>(&agr;<SUB>i</SUB>+&ggr;<SUB>i</SUB>−&tgr;<SUB>i</SUB>)<SUP>2</SUP> (4)

Electrostatic energy of DNA

The electrostatic interaction between free DNA double helix segments is described by integrating the solution of the Debye-Hückel equation for a point charge over two charged line segments:
E<SUP>(<UP>e</UP>)</SUP><SUB><UP>ij</UP></SUB>=<FR><NU>&ngr;<SUP>2</SUP></NU><DE>D</DE></FR> <LIM><OP>∫</OP></LIM> d&lgr;<SUB>i</SUB><LIM><OP>∫</OP></LIM> d&lgr;<SUB><UP>j</UP></SUB> <FR><NU><UP>exp</UP>(<UP>−&kgr;r<SUB>ij</SUB></UP>)</NU><DE>r<SUB><UP>ij</UP></SUB></DE></FR> (5)
D is the dielectric constant of water, and kappa  the inverse of the Debye length. rij is the distance between the current positions at the segments to which the integration parameters lambda i, lambda j correspond. The charge per unit length nu  is chosen such that the potential at the radius of the DNA coincides with the solution of the Poisson-Boltzmann equation for a cylinder with charge per length nu *0. For DNA in the presence of the Gouy layer of immobile counterions, this can be computed as nu *0 = qnu 0 in which nu 0 = -2e/Delta is the charge per length of the naked DNA (Schellman and Stigter, 1977), e is the proton charge, and Delta  = 0.34 nm is the distance between base pairs. Following Stigter, the value of q is 0.73 (Stigter, 1977). To save computation time, a tabulation of the double integral (5) is used. The table is parameterized by the distance of the segments and three values describing its relative orientation. During the simulation a linear interpolation of the tabulated values was used.

Internucleosomal potential

The internucleosomal interaction is modeled by a Gay-Berne potential with modifications by Kabadi (Gay and Berne, 1981; Kabadi, 1986a,b). The interaction energy of two particles with a center-to-center distance r is
V(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)=4&egr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)<FENCE><FENCE><FR><NU>&sfgr;<SUB>0</SUB></NU><DE>r−&sfgr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)+&sfgr;<SUB>0</SUB></DE></FR></FENCE><SUP>12</SUP></FENCE> (6)

<FENCE>−<FENCE><FR><NU>&sfgr;<SUB>0</SUB></NU><DE>r−&sfgr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)+&sfgr;<SUB>0</SUB></DE></FR></FENCE><SUP>6</SUP></FENCE>
The vectors ûi and û2 point into the direction of the symmetry axis of the particles. sigma 0 scales the potential width, and epsilon 0 scales the potential depth. The parameter chi  defines the anisotropy of the potential width, chi ' defines the anisotropy of the potential depth:
&sfgr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)=&sfgr;<SUB>0</SUB> <FENCE>1−<FR><NU>1</NU><DE>2</DE></FR> &khgr; <FENCE><FR><NU>(<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>+<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)<SUP>2</SUP></NU><DE>1+&khgr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)</DE></FR></FENCE>+<FENCE><FR><NU>(<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>−<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)<SUP>2</SUP></NU><DE>1−&khgr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)</DE></FR></FENCE></FENCE><SUP>−1/2</SUP> (7)

&khgr;=(&sfgr;<SUP>2</SUP><SUB>∥</SUB>−&sfgr;<SUP>2</SUP><SUB>⊥</SUB>)/(&sfgr;<SUP>2</SUP><SUB>∥</SUB>+&sfgr;<SUP>2</SUP><SUB>⊥</SUB>) (8)

&egr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)=&egr;<SUP>&ngr;</SUP>(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)&egr;′<SUP>&mgr;</SUP>(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>) (9)

&egr;(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)=&egr;<SUB>0</SUB>[1−&khgr;<SUP>2</SUP>(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)<SUP>2</SUP>]<SUP>−1/2</SUP> (10)

&egr;′(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>, <A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>, <A><AC>r</AC><AC>ˆ</AC></A>)=1−<FR><NU>1</NU><DE>2</DE></FR> &khgr;′ <FENCE><FR><NU>(<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>+<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)<SUP>2</SUP></NU><DE>1+&khgr;′(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)</DE></FR>+<FR><NU>(<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB>−<A><AC>r</AC><AC>ˆ</AC></A><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)<SUP>2</SUP></NU><DE>1−&khgr;′(<A><AC>u</AC><AC>ˆ</AC></A><SUB>1</SUB><A><AC>u</AC><AC>ˆ</AC></A><SUB>2</SUB>)</DE></FR></FENCE> (11)

&khgr;′=(&egr;<SUP>1/&mgr;</SUP><SUB>s</SUB>−&egr;<SUP>1/&mgr;</SUP><SUB>e</SUB>)/(&egr;<SUP>1/&mgr;</SUP><SUB>s</SUB>+&egr;<SUP>1/&mgr;</SUP><SUB>e</SUB>) (12)
sigma || is the relative potential width for particles oriented parallel and sigma perp for particles oriented orthogonal. epsilon s defines the relative potential width for particles in lateral (s, side-to-side) and epsilon e in longitudinal (e, end-to-end) orientation. nu  and µ are dimensionless parameters; generally one uses nu  = 1 and µ = 2.

Monte Carlo simulation

The classical Metropolis-Monte Carlo procedure is used to generate randomly a statistically relevant set of representative configurations of the system at temperature T (Metropolis et al., 1953). Starting with an arbitrary configuration, a new configuration is created randomly from the previous configuration (see below). The difference Delta E between the energy of the old and the new configuration is computed. If Delta E < 0, the new configuration will be accepted; if Delta E > 0, the new configuration will be accepted if a random number 0 < zeta  < 1 is below the threshold e-Delta E/kBT. For details see the extensive literature (e.g., Binder and Heermann, 1988; Allen and Tildesley, 1987, 1993).

As a trial move we used the pivot and the rotation move (Freire and Horta, 1976; Baumgärtner and Binder, 1979). In the literature both moves were shown to be correct and effective. In the pivot move, one segment is selected randomly as the origin of a rotation of a random angle of the interval [-delta , delta ] around a random axis. For a rotation move, a segment is chosen randomly. The end point of this segment is rotated around an axis through the start point of this segment and the end point of the next segment by an angle chosen randomly from the interval [-phi , phi ].

Structure and testing of the simulation program

The simulation program was developed in an object oriented manner in the programming language C++ (Stroustrup, 1997). All functions of each object were tested separately for correctness. The cooperation of the objects was tested thoroughly. We used the high performance debugging tool insure++ (Parasoft, Monrovia, CA) for verifying the formal correctness of the code. We implemented two different mechanisms for checking the self-consistency of the objects. In the testing phase and the early production phase we used them to check all parts of the program during run time. Furthermore, we computed simplified models like chains with only stretching interaction or chains with stretching and bending interaction by setting the parameters of the other interactions to 0. The computed values of the end-to-end distance agreed with the theoretical exact values.

The program was made parallel using POSIX threads, yielding a typical speedup of 1.9 using 2 processors and 3.3 using 4 processors.

To check for a possible influence of the used R250 random number generator (r250) we computed also some fibers with the built-in C-random number generator and the R16807 random number generator. Comparing the end-to-end distances of the different calculations, we found no differences within the range of the statistical error.

Correlation lengths and number of simulation steps

Configurations computed with the given Monte Carlo algorithm are correlated. To check whether enough simulation had been done to ensure good statistics, we computed the autocorrelation time of the energy, end-to-end distance, and mass density of the fibers. The autocorrelation function CAA function of a quantity A is defined as
C<SUB><UP>AA</UP></SUB>(j)=<FR><NU>⟨A(k)A(k+j)⟩−⟨A⟩<SUP>2</SUP></NU><DE>⟨A<SUP>2</SUP>⟩−⟨A⟩<SUP>2</SUP></DE></FR> (13)
j and k are indices that number the computed configurations. The integrated autocorrelation time tau int, A is defined as
&tgr;<SUB><UP>int,A</UP></SUB>=<FR><NU>1</NU><DE>2</DE></FR>+<LIM><OP>∑</OP><LL><UP>j</UP>=1</LL><UL>∞</UL></LIM> C<SUB><UP>AA</UP></SUB>(j). (14)
Two configurations at a distance of more than 2tau int, A can be considered statistically independent (Sokal, 1996). We computed tau int, A for a typical trajectory where each configuration was stored and for trajectories where each one-thousandth was stored to detected long range correlations. The maximal value of the computed tau int, A in Table 1 is 2500, which means that after 5000 simulation steps the configuration can be considered statistically independent from the previous one. Similar results were found in the other calculations.


                              
View this table:
[in this window]
[in a new window]
 
TABLE 1   Integrated correlation times tau int in number of Monte Carlo steps from a simulation with alpha  = 26°, beta  = 110°

For each simulation we performed at least 2.5 × 106 simulation steps. The first 0.5 × 106 steps corresponding to 100 independent configurations were not included in the analysis of the data because they are needed for relaxation. A typical computation of a polynucleosome with 100 nucleosomes took ~2 weeks of computer time on a 500-MHz Intel Pentium III processor.

Data analysis

First we define an approximate backbone of the fiber as the connecting line of the geometric centers of all chromatosomes and DNA segments in a window of 40 segments, which is moved in steps of eight segments over the chain. To avoid artifacts at the ends, the eight outermost segments are ignored. The contour length Lc is defined as the sum of the distances between adjacent backbone points. On both ends of the fiber, 20 chromatosomes are checked whether they are located inside a half space with the boundary surface orthogonal to the last segment of the backbone. These chromatosomes plus the number of the remaining chromatosomes are the effective number of chromatosomes Neff. Thus, the mass density is Neff/Lc. The mean radius of the fiber is defined as the average distance of all chain segments to the nearest segment of the backbone. The mean tilt of the DNA segments or the chromatosomes is computed from the angle between the DNA, respectively, chromatosome axis, to the nearest segment of the backbone. The persistence length lB of a fiber can be computed as
⟨u(s+s′)u(s′)⟩=e<SUP><UP>−s/l<SUB>B</SUB></UP></SUP> (15)
in which u(s) is the tangent vector to the fiber at the contour length s (Doi and Edwards, 1986). We compute this by interpolating natural cubic splines through the points of the backbone (Press et al., 1992). From these splines we determine the tangent vectors at the ends of the backbone. From the mean value of the product of these two vectors, we compute lB according to Eq. 15.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Parameters used in the simulation

Before describing the results, we shall discuss the parameters used for the simulation, in particular the internucleosome interaction potential. An overview of these parameters is given in Table 2. The bending and torsion elasticities correspond to a persistence length of 50 and 75 nm, respectively. For the DNA stretching modulus, we used the value determined in stretching experiments with laser tweezers (Smith et al., 1996). The elasticity of the chromatosome has not yet been determined experimentally; therefore we assumed that the chromatosome is stiff. Because a stiff segment is numerically difficult to handle, we chose a torsion potential 5 times larger than for DNA and the same stretching potential as for DNA. The DNA is assumed to be tightly bound to the chromatosome (Morse and Cantor, 1985).


                              
View this table:
[in this window]
[in a new window]
 
TABLE 2   Elastic and interaction parameters used in the simulations

The anisotropy of the potential depth of the internucleosomal interaction potential chi ' can be estimated as described for prolate ellipsoids (Gay and Berne, 1981). The disk is approximated by six Lennard-Jones spheres with radius 5.5 nm as in Fig. 2, and the minimum of the potential is computed for longitudinal and lateral orientation of the disks. Nearly independent of the radial orientation of the spheres one gets epsilon s/epsilon e = 1/5 and chi ' = 0.383. 



View larger version (32K):
[in this window]
[in a new window]
 
FIGURE 2   Approximation of a disk by spheres for calculation of the anisotropy of the internucleosome interaction potential. The total potential is approximated by summing over the Lennard-Jones contributions of the individual spheres. Two different possibilities of radial orientation are shown.

The actual depth and the position of the potential minimum can be obtained from the properties of liquid crystals of nucleosome core particles (Leforestier and Livolant, 1997). Under suitable conditions, core particles form a hexagonal-columnar phase with a distance of 11.55 ± 1 nm between the columns and a mean distance of 7.16 ± 0.65 nm between the particles in one column. We assume that these distances correspond to the minima of the internucleosomal potential. The minimum of Eq. 6 in lateral-lateral orientation is at r<UP><SUB>min</SUB><SUP>lat</SUP></UP> = 21/6 sigma 0 and in longitudinal orientation at r<UP><SUB>min</SUB><SUP>long</SUP></UP> sigma 0 (21/6 - 1 + (1 - 2chi /(1 + chi ))). This leads to sigma 0 = 10.3 nm and chi  = -0.506. The potential depth is estimated from a comparison of the experimental data with a computer simulation of a Gay-Berne discogen (kappa  = 0.345, kappa ' = 0.25, nu  = 1, µ = 2) at a scaled density of rho * = sigma <UP><SUB>0</SUB><SUP>3</SUP></UP>/V = 2.5. In the simulations the discogen undergoes a phase transition from the nematic phase to the hexagonal-columnar phase at T* = kBT/epsilon 0 = 4. Because the parameters of this Gay-Berne-discogen are similar to the parameters used above, we expect a phase transition at approximately the same T*. Thus we can chose epsilon 0 = (1/4) kBT. This value is equal to a potential depth of 5 × (1/4) kBT = 1.25 kBT for both core particles oriented axially, which is not too far from the value of 2 kBT determined recently using laser tweezers (Cui and Bustamante, 2000).

The geometry of the connection between the DNA and the chromatosome is sketched in Fig. 3. We assume that the linker DNA is either connected directly to the chromatosome (defining the chromatosome as a histone core with exactly two superhelix turns) or alternatively to a nucleosome stem, which is rigidly coupled to the rest of the chromatosome. In the following, we first analyze the situation where the linker DNA is connected via a stem ("stemmed" chromatosome). The entry and exit points of the DNA have a vertical distance (in the direction of the cylinder axis) of 3.1 nm and a radial distance of 8 nm from the symmetry axis of the core particle.



View larger version (15K):
[in this window]
[in a new window]
 
FIGURE 3   Geometry of the linker DNA-chromatosome connection. alpha  is the angle between linker DNA arms, beta  is the twist angle between adjacent nucleosomes, and d is the vertical offset between incoming and outgoing linker DNA.

The total length of the free linker DNA segments between two chromatosomes in our model depends on the nucleosome repeat length and on the coupling between flexible linker DNA and chromatosome (i.e., the presence or absence of the nucleosome stem). The precise path of the DNA outside the chromatosome is unknown; however, the length of the flexible linker may be estimated from the geometry of the core particle and the stem motif. The core particle contains 146-bp DNA in a superhelix with 1.65 turns, thus two turns correspond to 177 bp. Therefore, without a stem the length of the DNA connecting two chromatosomes equals the repeat length less 177 bp. The nucleosome stem ends at a distance of ~8 nm from the center of the core particle (Bednar et al., 1998). Assuming that the DNA follows a straight line from the exit point of the core particle to the end of the stem (Fig. 4), the total length of the DNA is 146 bp + 2 × 22 bp = 190 bp. The length of the DNA connecting two stemmed chromatosomes is then the repeat length less 190 bp.



View larger version (9K):
[in this window]
[in a new window]
 
FIGURE 4   Sketch for computation of the DNA content of a chromatosome with stem. R is the nucleosome core radius, d is the distance between chromatosome axis and the end of the stem, and l is the length of the DNA that is not bound to the core particle.

We defined the entry-exit angle of the linker DNA (alpha ) and the twist angle between adjacent nucleosomes (beta ) for elastic equilibrium in the absence of electrostatic forces and thermal fluctuations. Studies by cryoelectron microscopy have determined alpha  to 35° at 80 mM Na+ (Bednar et al., 1998). Because electrostatic and entropic forces change this angle in the simulation, we had to adjust the initial value of alpha  in the absence of external forces so that the effective angle measured on the simulated structure agrees with the experimental value. To this aim we performed simulations of polynucleosomes with 100 nucleosomes and a linker length of 10 bp, corresponding to a repeat length of 200 bp (e.g., in rat liver) at 100 mM NaCl, varying alpha  and keeping beta  = 80° constant (Fig. 5). For an elastic equilibrium angle alpha  = 26°, the effective entry-exit angle alpha sim in the simulations coincides with the experimental value of 35° to 40°. In subsequent simulations, this value was used. beta  is not known directly from experiments and had to be determined indirectly. Therefore, we performed simulations at different beta  and compared the mass density with the experimental value (Fig. 6). At beta  = 110°, the mass density of the fiber is 6.1 nucleosomes per 11 nm, which is in agreement with the experimental value of 6 nucleosomes/11 nm (Gerchman and Ramakrishnan, 1987). We therefore used this value in the subsequent simulations, if not stated otherwise.



View larger version (7K):
[in this window]
[in a new window]
 
FIGURE 5   Effective entry-exit angle of the linker DNA alpha sim in the simulated fiber as a function of the elastic equilibrium angle alpha . beta  is 80°.



View larger version (8K):
[in this window]
[in a new window]
 
FIGURE 6   Mass density of a simulated chromatin fiber (100 nucleosomes with stem, linker length 11 bp, alpha  = 26°) as a function of the twist angle between adjacent nucleosomes (beta ).

Structural properties of the simulated chromatin fibers

A visualization of the simulated structure (Fig. 7) shows that the computed fiber has an external aspect very similar to the solenoid model, although the nucleosomes are not ordered in a superhelix. The diameter of this fiber is 31.8 ± 0.1 nm, the mean tilt angle between the nucleosome axis and the fiber axis is 50.7° ± 0.7° and between linker DNA and fiber axis 92.5° ± 0.6°. The persistence length is 265 nm.



View larger version (28K):
[in this window]
[in a new window]
 
FIGURE 7   Visualization of two configurations of a simulated 100-nucleosome fiber, starting from a fully condensed straight configuration. Nucleosomes with stem, linker length 11 bp, alpha  = 26°, beta  = 110°.

If the angle beta  is not kept fixed along the fiber, but is set for each pair of chromatosomes randomly according to a Gaussian distribution, very similar results are found within a large range of the width of the distribution for most parameters (Table 3). We performed six simulations for systems with a width (1sigma ) of 10°, 20°, and 30°. The radius and the mass density are nearly unchanged. However, the value of the persistence length decreases dramatically while increasing the width of the beta  distribution (Fig. 8). At a width of 30°, the persistence length is only 28 nm, which is ~10% of the value at 0°.


                              
View this table:
[in this window]
[in a new window]
 
TABLE 3   Mean radius, mass density, and persistence length of simulated polynucleosomes with random beta  along the fiber



View larger version (9K):
[in this window]
[in a new window]
 
FIGURE 8   Persistence length of the simulated fiber as a function of the 1sigma -width of the beta  distribution. One-hundred nucleosomes with stem, linker length 11 bp, alpha  = 26°, beta  = 110°.

The previous calculations were started from a condensed structure, where the sum of the elastic potentials is zero. To check whether the final result depends on the initial structure or represents some local minimum, we performed simulations starting from a configuration where all segments are ordered in a straight line (Fig. 9) In this case, a structure very close to the final equilibrium configuration is reached after 160,000 Monte Carlo steps (Fig. 9 e); energy equilibration, however, takes ~5 × 105 steps (data not shown). For three simulations over 5 × 106 steps each, we obtained values for all tested properties (persistence length, diameter, mass density, total energy) that agreed with those starting from the condensed initial configuration within statistical error.



View larger version (9K):
[in this window]
[in a new window]
 
FIGURE 9   Simulation of the condensation of a 100-nucleosome fiber, starting from a stretched bead-string structure (a), and after (b) 15,000, (c) 30,000, (d) 45,000, (e) 60,000, and (f) 160,000 steps. Nucleosomes with stem, linker length 11 bp, alpha  = 26°, beta  = 110°.

Because the strength of the internucleosomal interaction was not measured directly, we had to estimate it from experiments on liquid crystals of nucleosome core particles (Leforestier and Livolant, 1997) as described above. Because this estimation is rather crude, we studied the sensitivity of the model to changes in the strength of the internucleosomal interaction. We performed simulations at values of epsilon 0 of the Gay-Berne potential from 0.1 kT up to 0.4 kT. In the axial direction, this corresponds to a potential depth of 0.5 to 2 kT. The data exhibit only small changes of the mass density (Fig. 10) and the persistence length (data not shown) depending on the value of epsilon 0.



View larger version (10K):
[in this window]
[in a new window]
 
FIGURE 10   Mass density of the simulated fiber as a function of the depth of the internucleosomal interaction potential epsilon 0. One-hundred nucleosomes with stem, linker length 11 bp, alpha  = 26°, beta  = 110°.

The calculations described so far were done for stemmed chromatosomes with a linker length of 10 bp, corresponding to a repeat length of 200 bp. The mass density of polynucleosomes with stem and a linker length of 22 bp, corresponding to a repeat length of 212 bp (e.g., in chicken erythrocytes) has its highest value of 4.6 ± 0.1 nucleosomes/11 nm at beta  = 110°. The mean diameter of this fiber is 36.1 ± 0.4 nm, the tilt of the nucleosome axis to the fiber is 55 ± 2°, and of the DNA 84 ± 2°, the persistence length is 49 nm.

At a repeat length smaller than 200 bp, it is unlikely that nucleosomes have a stem (see above). Therefore, we performed calculations of stem-less polynucleosomes with linker DNA of a length of 15 bp directly connected to the chromatosome, corresponding to a repeat length of 192 bp (e.g., in HeLa cells). Here we chose alpha  = 30°. At beta  = 120°, the mass density is 5.6 ± 0.1 nucleosomes/11 nm. If beta  is increased above this value, a stable fiber is no longer formed. The mean diameter of the fiber is 28.0 ± 0.1 nm, the persistence length 265 nm, the mean tilt of the nucleosome axis to the fiber axis is 55.2° ± 0.7° and between the linker DNA and the fiber axis 85.6° ± 0.8°. These values are in agreement with the experimental data. However, for stemless polynucleosomes with a linker length of 40 bp corresponding to a repeat length of 217 bp, it was not possible to obtain stable fibers at any chosen value of beta  and epsilon 0.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Several experimental structural properties of the chromatin fiber at physiological ionic strength could be reproduced in our simulation.

The general aspect of the fiber resembles that of the classical solenoid model: nucleosomes stack to form a helix with a diameter of ~30 nm. However, other than in the solenoid model, the linker DNA is straight and crosses the fiber, so that next-nearest neighbors stack onto each other. Experimentally, the radius of the fiber was determined to ~30 nm by electron microscopy (Finch and Klug, 1976). The computed values of all fibers are between 28 and 36 nm, which is in good agreement with this value.

Data on the linear mass density of the fiber are controversial in the literature. An overview of older data can be found in the textbook by van Holde (1989). The most reliable data of linear mass density of the fiber were gained with neutron scattering and scanning transmission electron microscopy experiments. At physiological ionic strength one finds a value of 6 nucleosomes per 11 nm (Gerchman and Ramakrishnan, 1987). In the cited experiments, spermine was used to prepare the samples. Spermine binds tightly to DNA; it might actually lead to an additional compaction of the chromatin fiber because older neutron scattering experiments without spermine show lower values of the linear mass density at corresponding ionic strength (Suau et al., 1979). On the other hand, light scattering data (Campbell et al., 1978) agree well with the data from Gerchman and Ramakrishnan (1987). Furthermore, in all cited experiments chromatin was prepared by micrococcus nuclease digestion; one could speculate that preferentially open chromatin structures are digested and only the more condensed structures are left. In summary, a value of 6 nucleosomes per 11 nm or somewhat smaller seems to be compatible with the experiments. The values of the computed fibers are between 4.6 and 6 nucleosomes per 11 nm, which is in good agreement with these data. A value of 6 nucleosomes per 11 nm was also used to fit the tilt angle between adjacent nucleosomes beta . However, it is not trivial that this value of compaction of the fiber can be reached, as it can be seen from the fact that for a fiber with a repeat length of 212 bp the maximal value found in our simulations was 4.6 nucleosomes per 11 nm.

The angle between nucleosomes and linker DNA to the fiber axis can be determined from experiments using linear dichroism. In the past, rather different results were reported by various groups, but most were later attributed to preparation artifacts (Dimitrov et al., 1990). Through a consistent reinterpretation of the valid data of flow linear dichroism and electric linear dichroism, the angle between the linker DNA and the fiber axis was estimated to be 50° to 90° and the angle between the axis of the nucleosomal DNA superhelix and the fiber axis to be 45° to 60°. These data are in very good agreement with the values found for the simulated fibers.

The experimental value of the persistence length of the chromatin fiber depends on the technique that was used to determine its value. Analysis of the distances between genetic markers in nuclei from human fibroblasts (repeat length approx 190 bp) obtained from experiments using fluorescence in situ hybridization (van den Engh et al., 1992; Yokota et al., 1995) using a Flory (statistical segment) model results in a value of 100 to 140 nm, an analysis using the Porod-Kratky wormlike chain model in 70 to 110 nm (Mehring, 1998). The persistence length was also determined by an analysis of the end-to-end distances of chromatin fibers of unknown origin bound to a mica surface and measured with scanning force microscopy. The analysis results in a value of 30 to 50 nm (Castro, 1994) as cited in Houchmandzadeh et al. (1997). It is unclear whether the binding of the fiber to mica causes artifacts. Very recent experiments stretching single chromatin fibers from chicken erythrocytes with laser tweezers resulted in a value of 30 nm (Cui and Bustamante, 2000) at low salt concentrations; no data for the persistence length was given at physiological salt. In summary, it seems that for short repeat lengths the persistence length is in the range of 30 to 140 nm. The values from our simulations of regular fibers with shorter repeat lengths (192 and 200 bp) are 200 and 265 nm, which is larger than the experimental value, whereas the value of 49 nm for larger repeat length (212 bp) agrees well. The persistence length of fibers with smaller repeat length decreases dramatically if we increase the width of the distribution of beta ; a 1-sigma width of 20° to 30° yields a persistence length that fits very well to the experimental values. This irregularity in beta , i.e., in the twist between successive nucleosomes, should in principle be implemented in the model together with a corresponding change in linker length (to keep the correct orientation of the DNA on the nucleosome). However, in our model DNA is approximated by a homogeneous elastic rod whose rotational orientation on the histone core does not enter the model. Therefore, for simplifying the simulation, we have chosen to neglect the linker length change, taking into account the variation in beta  only. Because the twist variation of 20° to 30° necessary to reproduce the experimental persistence length of the chromatin fiber would correspond to a variation in linker DNA length of only ±1 bp, neglecting this length change should not have dramatic consequences on the structure of the fiber as long as the twist variation is accounted for correctly.

The ionic environment of the solvent is known to have an important influence on the conformation of the chromatin fiber. To assess the dependence of our model on the ionic strength, we performed simulations of a stemmed polynucleosome with alpha  = 26°, beta  = 110°, and a linker-DNA of 11 bp. The mass density of the simulated fiber does decrease slightly with lowering the ionic strength (data not shown); however, the measured decrease below 70 mM NaCl, which is much more pronounced (Gerchman and Ramakrishnan, 1987), could not be reproduced quantitatively. This is not surprising because we neglect in our model all effects of the ionic strength on the nucleosomal interaction and on model parameters such as alpha , beta . Thus, the structural conclusions that we draw from our model are for the moment only valid at the near-physiological ionic conditions used in most of our simulations. Here, the total contribution of the electrostatic interaction (i.e., DNA-DNA repulsion) to the chain energy is ~18% (data not shown); the elastic, Gay-Berne contributions, and therefore the geometrical effects dominate.

Changing the ionic strength will most probably act on the chromatin conformation by changing the nucleosome-nucleosome interactions; however, this would require a parameterization of the ionic strength dependence of this interaction. This can in principle be implemented in our model by changing, e.g., the parameters of the Gay-Berne potential in such a way that the experimental dependence of the linear mass density on salt concentration is reproduced. However, a systematic investigation of the electrostatics will require another extensive set of simulations; together with a study of the dependence of the force-extension curve on salt concentration, these calculations are currently underway and will form the subject of a forthcoming publication.

In a recent Brownian dynamics simulation where the interaction potential between nucleosomes was described by a point-charge approximation based on the known crystal structure, Beard and Schlick (2001b) could demonstrate an opening of the chromatin fiber structure when the ionic strength was lowered from 50 to 10 mM. However, the complexity of the interaction potential (Beard and Schlick, 2001a) as well as the Brownian dynamics algorithm will make the simulation much more computationally intensive than a Monte Carlo approach. For systematic studies of effects of linker DNA geometry on the structure of the chromatin fiber, as well as for simulation of much larger structures (such as entire chromatin loops of a size of 100-200 kb), we therefore prefer the Monte Carlo procedure used here.

Another Monte Carlo model for the chromatin fiber was proposed recently (Katritch et al., 2000). The nucleosomes were approximated by spheres interacting through simple steplike potentials, and electrostatic interactions were not included explicitly. Whereas the authors were able to reproduce the stretching experiments on single chromatin fibers with laser tweezers at 5 and 40 mM NaCl, no data at physiological ionic strength was shown. Also, no results of the model for other well-known properties of the fiber, such as mass density or diameter, were given in that work.

A long-standing debate in chromatin research is whether the shape and properties of the fiber are regulated by the internucleosomal forces (Luger et al., 1997) or geometrical constraints (Krajewski and Ausió, 1994). Our simulations suggest that changes in the connection geometry between subsequent nucleosomes might have a much larger influence on the mass density than changes in internucleosomal interaction forces. Further studies of this point, as well as a more precise analysis of the persistence length and the salt dependence of structural parameters will require the development of more detailed potentials for the electrostatic interactions of the DNA at short distances and for the internucleosomal interactions.

    ACKNOWLEDGMENTS

The authors are grateful to Peter Grassberger (BUGH, Wuppertal) for helpful discussions and support. We thank Christian Münkel (SAP AG, Walldorf, Germany) for important help in the beginning of this project, Karsten Rippe (DKFZ, Heidelberg) for stimulating discussions, and Katalin Tóth (DKFZ, Heidelberg) for critical reading of the manuscript. We thank the DKFZ (Heidelberg), the IWR (Heidelberg), and the University of Mannheim for supplying us generously with computer power. This work was supported by Grant 01 KW 9620 of the German Ministry for Education and Research (German Human Genome Project DHGP).

    FOOTNOTES

.

Address reprint requests to Jörg Langowski, German Cancer Research Center (DKFZ), Division Biophysics of Macromolecules (H0500), Im Neuenheimer Feld 280, 69120 Heidelberg, Germany. Tel.: 49-6221-423390; Fax: 49-6221-423391; E-mail: joerg.langowski{at}dkfz.de

Submitted January 25, 2001, and accepted for publication February 22, 2002.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES