Biophysical Journal 84:3607-3623 (2003)
© 2003 The Biophysical Society
Azimuthal Frustration and Bundling in Columnar DNA Aggregates
H. M. Harreis,
C. N. Likos and
H. Löwen
Institut für Theoretische Physik II, Heinrich-Heine-Universität Düsseldorf, Düsseldorf, Germany
Correspondence: Address reprint requests to Holger M. Harreis, Institut für Theoretische Physik II, University of Düsseldorf, Universitätsstr. 1, Düsseldorf, Germany 40225. Tel.: 49-211-811-3701; Fax: 49-211-811-2262; E-mail: harreis{at}thphy.uni-duesseldorf.de.
 |
ABSTRACT
|
|---|
The interaction between two stiff parallel DNA molecules is discussed using linear Debye-Hückel screening theory with and without inclusion of the dielectric discontinuity at the DNA surface, taking into account the helical symmetry of DNA. The pair potential furthermore includes the amount and distribution of counterions adsorbed on the DNA surface. The interaction does not only depend on the interaxial separation of two DNA molecules, but also on their azimuthal orientation. The optimal mutual azimuthal angle is a function of the DNA-DNA interaxial separation, which leads to azimuthal frustrations in an aggregate. On the basis of the pair potential, the positional and orientational order in columnar B-DNA assemblies in solution is investigated. Phase diagrams are calculated using lattice sums supplemented with the entropic contributions of the counterions in solution. A variety of positionally and azimuthally ordered phases and bundling transitions is predicted, which strongly depend on the counterion adsorption patterns.
 |
INTRODUCTION
|
|---|
Many biological systems contain densely packed DNA assemblies, as, for example viral phage heads and sperm. For the proper functioning of these biological systems, including humans, it is of extreme importance that the mechanisms carrying out the packaging of DNA in the cell work in a robust manner, since, for example, it is believed that DNA packing in chromatin plays an important role in gene regulation (Wolffe, 1992
). In light of the rapidly growing field of gene therapy it is of great interest to understand the mechanisms actually responsible in living organisms for condensing DNA into densely packaged assemblies. The first step to this end is a model of DNA which is able to capture its most significant characteristics, with the second step consisting of devising a theory for DNA assemblies. In the last few years, many efforts have been made on the theoretical side to understand the interaction of two DNA molecules and DNA condensation with a variety of methods, including molecular dynamics and Brownian dynamics simulations (Grønbech Jensen et al., 1997
; Kornyshev and Leikin, 1997
; Ha and Liu, 1997
; Podgornik and Parsegian, 1998
; Shklovskii, 1999
; Kornyshev and Leikin, 1999
; Sottas et al., 1999
; Nguyen et al., 2000
; Allahyarov and Löwen, 2000
; Ha and Liu, 2001
). The matter is complicated by the fact that due to its chemical structure DNA is a helical molecule, rendering solutions for the DNA-DNA interaction considerably complicated. Moreover, the overall electroneutrality condition dictates that counterions be present in the solution, and the latter screen the electrostatic repulsion between the DNA rods. Only far from their axes can DNA molecules be apprehended as uniformly charged cylinders: this is the simplest approximation possible in an investigation of the DNA-DNA interaction and one that neglects the helical symmetry completely (see, for example, Grønbech Jensen et al., 1997
; Levin et al., 1999
; Hansen and Löwen, 2000
; Levin, 2002
; Strey et al., 1999
, 1997
; Oosawa, 1971
; Stigter, 1977
; Manning, 1978
; Frank-Kamenetskii et al., 1987
; and references therein). It has to be expected that such an approximation works well for distances much larger than the scale of the helical symmetry of the DNA molecule, R >> H, where H
3.4 nm is the DNA pitch length. This approach amounts to calculating the interaction of two homogeneously charged cylinders, whereby the continuously smeared charges along the cylinders create an electrostatic repulsion of two DNA molecules (exponentially screened by the electrolyte). Indeed, predictions for force-distance curves on the grounds of a traditional Derjaguin-Landau-Verwey-Overbeek theory for homogeneously charged cylinders turned out to be accurate for separations larger than several nanometers, whereas significant deviations in the biologically more relevant range of smaller separations (Kornyshev and Leikin, 1997
) emerged. It can be concluded that apart from investigations where only the far-field behavior is of importance, it is crucial to consider the helical symmetry of DNA molecules, since the interaction potential in the relevant regime of intermediate distances is dramatically changed by the presence of a highly inhomogeneous charge distribution. An additional effect is provided by the fact that DNA is a polyelectrolyte molecule; in an aqueous solution, its cations dissolve into the solution, leaving behind a negatively charged DNA phosphate backbone. A major fraction of the cations condenses in the Bjerrum layer (Manning, 1978
) around the molecular surface. With cations specifically adsorbing onto the DNA surface present in the solution, however, the scenario changes; the DNA molecules can be fully neutralized (Wilson and Bloomfield, 1979
; Widom and Baldwin, 1980
; Heath and Schurr, 1992
) or even overcharged (Pelta et al., 1996
). The interaction potential is thus additionally influenced by the amount and type of counterions present in the solution.
To condense DNA in an aggregate, either osmotic stress (Rau et al., 1984
) or counterions specifically adsorbing on DNA have to be applied as condensing agents (Bloomfield, 1996
). The latter can be, e.g., salts with Mn2+, Cd2+, spermidin, protamine, or cobalt hexammine (Bloomfield, 1996
) cations, which are known to preferentially adsorb in the DNA grooves (Tajmir-Riahi et al., 1993
; Fita et al., 1983
; Hud et al., 1994
; Shui et al., 1998
). The sensitivity to the type of counterion for DNA aggregation (Bloomfield, 1996
) is manifest in the fact that other counterions, such as, e.g., Ca2+ or Mg2+, which are known to exhibit a high affinity to phosphates and thus predominantly adsorb on the strands, do not induce DNA aggregation. A model should thus incorporate/reproduce these subtle effects and be able to explain the mesomorphism (Podgornik et al., 1998
) of DNA aggregates stemming from the presence of different types of counterions.
Once the interaction of DNA molecules is derived by means of some theory, one can turn to the next step and calculate the properties of DNA assemblies. The structural organization and properties of such condensates in vivo are largely unknown but have been, in the last several years, under investigation in in vitro experiments (Robinson, 1961
; Wilson and Bloomfield, 1979
; Livolant and Bouligand, 1986
; Livolant, 1991
; Rau and Parsegian, 1992a
,b
; Rill et al., 1991
; Ma and Bloomfield, 1994
). Simple model systems able to predict the spatial as well as the orientational structure of these condensates are highly desirable for a better elucidation of the mechanisms occurring in vivo. Previous work has shown (Kornyshev and Leikin, 1998a
) that it is a reasonable approximation/simplification to focus on columnar assemblies, neglecting possible tilting effects, as we will explain later. Most of the work relied on approximating DNA as homogeneously charged rods (Grønbech Jensen et al., 1997
; Ha and Liu, 1997
; Podgornik and Parsegian, 1998
; Shklovskii, 1999
; Nguyen et al., 2000
). Only when taking into account, however, the helicity of DNA molecules, a relevant feature for the properties of such a columnar DNA assembly emerges: a nontrivial interplay between the torsional and translational degrees of freedom.
A mean-field calculation of this problem was presented in Lorman et al. (2001)
, whereas the full statistical mechanical problem of columnar DNA assemblies was recently solved in Harreis et al. (2002)
using a pair potential for the DNA-DNA interaction devised in Kornyshev and Leikin (1997)
. The motivation for the present article is twofold: first, we give more details and background for the calculations already published in Harreis et al. (2002)
. In this work, it was found that the dependence of the optimal azimuthal orientation angle of two DNA molecules on their interaxial separation gives rise to azimuthal frustrations in an aggregate, thereby inducing phase transitions between different ordered orientational structures. Furthermore, depending on the type and amount of counterions condensed on the DNA surface, strong attractions were found, resulting in DNA bundling transitions. More importantly, the second motivation for the present work is to discuss the effect of discretized charges along the DNA strands and the effect of the dielectric jump at the DNA surface on the phase behavior. We find that although the phase boundaries shift quantitatively, especially at high densities, the global topology of the phase diagrams remains unaffected. This gives evidence for the fact that the topology of the phase diagram itself is generic, i.e., will be stable also with respect to further changes in the interaction, including, for example, hydration forces that are sometimes modeled through a distance-dependent dielectric constant
(Lee et al., 2002
).
 |
THE MODEL
|
|---|
DNA is a helical biomolecule with two charged phosphate strands helically winding around a core region consisting of nucleotide basepairs. The two strands are not symmetrically distributed around the molecule's core region, but rather are separated by an azimuthal angle of 2
s
0.8
, see Fig. 1 for an illustration. Under physiological conditions, DNA is present in the B-DNA conformation, a right-handed helical molecule (Saenger, 1984
). In B-DNA, there are N = 10 nucleotides per helical turn with a helical pitch length of H
34 Å. Each nucleotide contains a negatively charged phosphate group, giving rise to a total charge of q = -10 e per helical pitch, which translates into a surface charge density of
= 16.8 µC/cm2. To model the interaction, we envision the molecules as long, rigid cylinders with a hard-core radius of a = 9 Å. Strictly speaking, this approximation is only appropriate for DNA fragments of contour lengths up to the persistence length Lp, which is typically found to be 500 Å1000 Å (depending on the ionic strength; see Kornyshev and Leikin, 2000
). Samples of parallel packed arrays have, however, been prepared for contour lengths of up to 100 Lp (Rau et al., 1984
; Podgornik et al., 1996
). In our model the phosphate backbone is accounted for by continuous helical line charges located on the surface of the DNA hard-core cylinder. We also calculated pair interactions for discrete charge patterns on the DNA surface, as we will discuss in detail later. Each DNA duplex furthermore carries a compensating positive charge stemming from the adsorbed counterions, which are modeled in the same way as the phosphate backbone as continuous line charges. The degree of charge compensation will be referred to as 0 <
< 1, whereas the fractions of condensed counterions in the minor and major grooves, and on the two strands, are accounted for by f1, f2, and f3, respectively, where f1 + f2 + f3 = 1 holds. The nonadsorbed, mobile counterions in solution screen the Coulomb interactions between the helices, causing at large separations an exponential decay of the latter with the Debye screening length
-1.
We wish, at this point, to discuss advantages and drawbacks of the present model that is characterized by a Debye-Hückel approach combined with the ion condensation model. Quite generally speaking, the two great advantages of the resulting Yukawa-type interaction are its great simplicity and its remarkable flexibility. Although not all situations, especially those of small separations might be accurately described quantitatively, possible deviations can be compensated by introducing the concept of effective, renormalized charges, as has been shown, for example, in Löwen (1994a)
. In the specific case of DNA-DNA interactions, previous work (Allahyarov and Löwen, 2000
) has investigated the question of DNA-DNA interaction in the framework of the primitive model of electrolytes. The authors thereby relied on microscopically resolved molecular dynamics simulations. The DNA molecules were modeled, as in the present study, as rigid cylinders having all structural parameters of B-DNA. More specifically, in the reference cited, the double helical charge pattern was incorporated via discrete charges exhibiting an effective diameter, placed on the cylindrical surfaces of the DNA molecules. The solvent was accounted for by its dielectric constant (
) as in the present work, but counter- and salt ions were explicitly included. It has been found in this reference that Yukawa-like effective interactions, resulting from a canonical tracing out of the microions, are capable of describing the potentials for DNA interaxial separations R > 25 Å. The authors furthermore showed that the behavior for interaxial separations R < 25 Å could be equally well described by a Yukawa potential; nevertheless, a different, separation-dependent effective charge has to be introduced for this range. To quantitatively reproduce the DNA-DNA pair potential, one is thus dealing with a Yukawa potential with a distant-dependent effective charge which saturates for R > 25 Å. It will, in general, be different from the input value of the Yukawa segment model as employed in the present work. It has thus to be expected that the R-dependent part of the pair interaction as shown in the present work will be affected by this charge renormalization. The angular part of the interaction will not be influenced, however, except for a scaling by an overall trivial factor, since the interaction is short-ranged and, as will be shown in what follows, the phase behavior is dominated by the nearest-neighbor interactions. The most important predictions of this study will not be affected: it is the location of the minimal azimuthal orientational angle between two DNA molecules that governs the frustrations and thus the equilibrium structure in the DNA aggregate. Furthermore, the phase diagrams in the case of repulsive interactions exhibit very small density jumps at the phase transitions, implying that the same effective charge, thus the same pair potential, can be employed in both phases. Although the absolute values of the free energies of the various phases will be affected by charge renormalization, the comparisons between those, and hence the location of the phase boundaries, will not. In the case of attractive interactions between two DNA molecules, the attractions occur only at a specific mutual azimuthal orientation. Since, as we have argued above, the latter remains unmodified by charge renormalization, it can be concluded that the same statement holds for the phase diagrams caused by these attractions. We could, as we will discuss in the following sections, furthermore illustrate that the predicted phase diagrams are, in their essential features, qualitatively robust against variations of the underlying pair potential, such as inclusion/exclusion of the dielectric jump at the DNA surface. We thus have good grounds to believe that although corrections to the pair potential are necessary, the predictions regarding macroscopic behavior are robust. Thus, the current approach captures the essential physics governing the biological phenomena at hand.
In our model, we study formally the two extreme cases of dielectric constants
1and
in the DNA core and in the solvent, respectively. The first case is that we assume no dielectric jump at all,
/
1 = 1, whereas the other limit is
/
1 =
. In the first case, it is more convenient to formulate the interaction in terms of a Yukawa-segment model, whereas the second case has been elaborated in a practical form by Kornyshev and Leikin (1999)
. The motivation to study different
/
1 is to check effects of the discontinuity formally. In reality one would expect
/
1
since the dielectric constant of bulk water is very high. Close to the DNA surfaces, however, it is not at all clear whether the effect of a dielectric discontinuity as described by macroscopic electrostatics is justified. More realistic dielectric effects were taken into account by a space-dependent dielectric constant
(Lado et al., 1998
). One could surmise that if the resulting interaction and phase behavior is similar for the two limiting cases
/
1 = 1 and
/
1=
, dielectric effects on this molecular scale are not actually very important at all. This in turn gives evidence for at least qualitative stability of our results under application of more realistic interactions stemming from more refined molecular calculations.
The main characteristics of the model DNA molecules are illustrated in Fig. 1. For clarity, possible condensed counterion strands have been omitted in the illustration. The azimuthal orientation of molecule i is referred to by its azimuthal angle
i, which is defined in the following way. A plane (shaded gray in Fig. 1) perpendicular to the parallel axes of the two DNA molecules hits the dark colored 5'3' strand (Sinden, 1994
) of each molecule at the point indicated by the vector originating from molecule's i axis, which we may formally call spin. The angle
i formed by this vector and some arbitrary reference direction on the plane, taken, for clarity, to be the vector connecting the two molecules' axes, is the azimuthal orientation angle of molecule i. We assume that the DNA molecules are parallel, as depicted in Fig. 1, which is justified by reasons given in A Theory for DNA Assemblies. If we furthermore assume the molecules to be infinitely long and their charge distributions to be described by helical line charges as illustrated in Fig. 1, their mutual state can be described by two parameters: their interaxial separation R as well as their mutual azimuthal orientation,
=
1 -
2. The problem thus reduces to an effective two-dimensional problem of x-y spins interacting via a potential U(R,
). We further illustrate this point in Fig. 2, which depicts the gray-shaded plane included in Fig. 1 in more detail. It has to be noted that the problem may only be viewed as effectively spatially two-dimensional under the assumption of continuous line charges. For discrete charge patterns, the orientations
1 and
2 both enter the pair potential. Let us assume discrete charges to illustrate the validity of this statement. The two molecules shall be separated by a vector R, as shown in Fig. 2, with molecule 1 at an angle
1 and molecule 2 at an angle
2 relative to R in a given plane P that perpendicularly cuts the molecular axes. The points where the 5'3' strands of molecules 1 and 2 hit the plane P shall be denoted by p1 and p2, respectively. Both in p1 and p2 discrete charges are located as is the case in Fig. 2. Were the interaction only to depend on the mutual azimuthal orientation,
=
1 -
2, a configuration with both molecules turned by 
should yield the same interaction. Obviously, after turning both molecules the 5'3' strands will hit P in new locations
and
. This means altered charge distances in the contribution of plane P to the total DNA-DNA interaction. Since the total interaction is the sum of the contributions of all charges (planes), the interaction might still be conserved if another plane along the molecule contributed the same value after the rotation as P did before the rotation and vice versa. The only plane capable to switch configuration with P through a rotation by 
is a plane P' shifted by
z = H 
/2
relative to P along the molecular axes. The 5'3' strands will then cut through P' in p1' = p1 and p2' = p2. With a discrete charge pattern, however, charges will only be located in p1' and p2' if
z is commensurate with the rise of the charge pattern along the molecular axes, or in other words, if 
= n 2
/N holds, with n
and N the number of DNA charges per helical pitch and strand. If on the other hand, continuous line charges are used, the original plane P and P' are equivalent without any further condition and contribute the same amount to the interaction. The only requirement to be met is that the molecules be at least one helical pitch long, so that the existence of P' is guaranteed. The mutual azimuthal angle
can, for continuous line charges under the additional condition of infinitely long molecules, equivalently be thought of as a relative vertical shift z = H
/2
of the two molecules. We will come back to a more detailed discussion on discrete charge patterns versus continuous line charges at a later point in the next section.
 |
THE PAIR POTENTIAL
|
|---|
As already sketched in the Introduction, the pair potential will be considered under different assumptions concerning dielectric jump and charge distributions. The approach is, on a general level, based on the linear screening theory picture, yielding a Yukawa-like, screened Coulomb interaction for any pair of charges on the two molecules (Schneider et al., 1985
, 1986
, 1987
). We will first resort to considering the case of no dielectric jump and refer to this situation as the Yukawa-segment-model potential. The Yukawa-segment idea has been tested against microion resolved simulations in Löwen (1994a
,b
) and has been used for calculating dynamical correlations in Tobacco-Mosaic Virus suspensions and phase diagram calculations of the latter in Kirchhoff et al. (1996)
and Graf and Löwen (1999)
, respectively. Here, the Yukawa-segment approach furthermore allows for testing the influence of a discrete charge pattern as opposed to continuous line charges. The second case includes the dielectric jump at the DNA surface, yet necessitates continuous line charges. We will refer to it in the following as Kornyshev-Leikin potential.
Yukawa-segment-model potential
The canonical starting point for the Yukawa-segment-model is to exactly mimic the discrete number of charges present in real DNA molecules. The second generic case, opposed to the former, is to assume the charge distributions to be continuous line charges. Although the first approach might, at first sight, seem superior to the latter, it has to be kept in mind that the real charge distribution will definitely be not pointlike, but rather smeared on the whole phosphate group, two charges of which are closely neighbored, so that a modulated continuous line charge distribution should be the most realistic way of modeling the DNA charge distribution. Such an approach, however, requires an input from quantum chemical calculations and is therefore beyond the scope of the present study. We will now first illustrate the general approach to the calculation of the pair potential and then come back to a discussion of the differences between the discrete and the continuous charge distribution version.
We assume linear screening to act between any two charge elements qi and qj on the continuous helical line charges of the DNA molecules, yielding a Yukawa interaction (Schneider et al., 1985
, 1986
, 1987
),
 | (1) |
Here,
=
is the inverse Debye screening length and
= 81 is the dielectric constant of the solvent (water). To access the total pair interaction of two DNA molecules, we have to integrate along each pair of interacting helical line charges (strands) (or sum in the case of discrete charge patterns).
Let molecule 1 be at the origin of the coordinate system and molecule 2 at R = Rx, see Fig. 2. In its most general form, a helix, parameterized by its helical angle
, furthermore depends on a set
of additional parameters. This set of parameters
consists of the helix radius a, the helical rise
= H/2
, the position (rx,ry) of the helix axis in the x-y plane and the angular offset 
of the helix, indicating where the helix starts to rise from the x-y plane. Making use of the special conditions present in our case, namely that we only consider molecules residing on the x-axis and that all helices exhibit the same radius as well as the same helical rise,
can be reduced to only consist of rx and 
,
= (rx,
). The corresponding helix parameterization for one single helix reads as
 | (2) |
The angular offset 
is set to
1 for the first strand on the first molecule. Thereby the angular offsets of all other strands involved are uniquely determined by the DNA geometry, for example, the second DNA phosphate strand on molecule 1 has 
= 2
s+
1, the counterion strand in the minor groove is characterized by 
=
s+
1 and the counterion strand in the major groove has 
=
- 2
s+
1 as angular offset. The charge strands on molecule 2 follow the same logic, except that their respective offsets have a term of
2 instead of
1, since the rotation of molecule 2 has to be accounted for; see again Fig. 2 for an illustration. The interaction between one strand on molecule 1 and another strand on molecule 2 is given by
 | (3) |
which is a diverging quantity, since the integral in Eq. 3 includes the two infinitely long strands. What we are interested in for our purpose is the interaction that segments of a given length L experience. As the Yukawa-type interaction between all charge segments decays exponentially and since, due to the periodicity, all helical pitches are the same, we may to this end, proceed in the following way. On molecule 1 one pitch length H serves as integration interval, whereas on molecule 2 we integrate from -
to
. Practically, due to the exponential decay in the potential, convergence of the integral is obtained after a maximum of 10 pitch lengths H has been integrated. The result is the interaction energy of one pitch on strand 1 with the total length of molecule 2. Multiplication of this quantity with the number of pitches L/H to be taken into account for a length L yields the interaction of a segment of length L on strand 1 with a segment of length L on strand 2, whereby endpoint effects are ignored via the integration from -
to
.
The total interaction of a segment of length L on molecule 1 with one on molecule 2, then, is the sum over the interactions of all strands on molecule 1 with all strands on molecule 2, including the DNA phosphate strands as well as the condensed counterion strands:
 | (4) |
where the symbolic notation above implicitly assumes i to be taken from the set of all strands on molecule 1 and j correspondingly from molecule 2. Inserting the Yukawa-segment interaction, Eq. 1 in Eq. 4 and carrying out the r and r' integrations in Eq. 3, together with the above consideration on the integration intervals, yields the expression
 | (5) |
Here and in Eq. 4, the index i is taken from the set i
{s1(1), s2(1), c1(1), c2(1)} and j covers j
{s1(2), s2(2), c1(2), c2(2)}, whereas
i shows the dependence of the given strand on the specific geometrical parameters determining its parameterization. By sk(l) the kth phosphate strand on the lth molecule is denoted, whereas ck(l) describes the corresponding counterion strand. In sk(l) the counterion strands which are condensed on the phosphate strands are included, since they only trivially renormalize the charge carried by the phosphate strands. This enters into the charge fraction parameters fi and fj in the following way:
 | (6) |
 | (7) |
 | (8) |
where f1, f2, and f3 are the fractions of counterions condensed in the minor and major grooves, and on the two strands, respectively, satisfying f1 + f2 + f3 = 1.
The differences of a discrete charge potential to a continuous line charge potential can be estimated by tuning the number of charges per pitch length, N. As we discussed in The Model, for discrete charges the interaction does depend on both molecules' orientations
1 and
2 and not only on the difference
=
1 -
2 as it is the case for continuous line charge distributions. For discrete charge patterns, this opens up two different routes: The first and simpler is to set
1 = 0 and look at U(R,
1 = 0,
2), whereas the second and more refined one is to vary
1 and
2 to then consider U(R,
) at
=
'1 -
'2, where
'1 and
'2 have been obtained as energetically optimal combination for a given mutual azimuthal orientation
of the two DNA molecules.
The first approach is taken in Fig. 3, where the pair interaction per persistence length Lp, U(R,
1 = 0,
2), is displayed as a function of the azimuthal orientation angle
2 with
1 = 0 fixed, at two fixed interaxial separations, R = 2.1 nm and R = 2.5 nm, for N = 10 and N = 20 charges, as well as for a continuous line charge. The counterion condensation parameters are f1=0.3, f2=0.7, and f3 = 0. It can be seen that already for N = 20 the obtained potential curve is indistinguishable from that for the continuous line charge potential at both interaxial separations. For N = 10 charges, on the other hand, deviations exist predominantly for R = 2.1 nm, but have decreased to a minuscule level for R = 2.5 nm. A more detailed structure of the pair potential as a function of the azimuthal orientation is apparent for the smaller separation. The differences mainly pertain to the region around the maximum and the two minima. The position of the global minimum, however, the most important parameter for the behavior in an assembly, is practically unchanged. This assertion is only based on the observation of the potential at two fixed interaxial separations. Its main point, however, is sustained by the data shown in the inset of Fig. 4, where we show the optimal azimuthal orientation angle
2,opt, again at f1 = 0 fixed, as a function of the DNA-DNA interaxial separation R. The corresponding potential is shown in the main graph of Fig. 4. The detailed behavior of the optimal azimuthal alignment angle is different for N = 10 charges from that for N = 20 charges and continuous line charges, whereas the latter are indistinguishable from one another. The key points for the behavior in an aggregate, however, remain unchanged for all cases: first, the optimal angle is nonzero for all interaxial separations smaller than R*
30 Å and second, for very small intersurface separations the optimal angle is
0.42
.

View larger version (19K):
[in this window]
[in a new window]
|
FIGURE 4 Yukawa-segment pair potential per length Lp as a function of the interaxial separation R of two DNA molecules, at the optimal angle 2,opt(R), depicted for N = 10, N = 20 charges as well as for continuous line charges. The dependence of the optimal angle on the interaxial separation R is shown in the inset.
|
|
The second approach to discrete charge patterns is to calculate the interaction for all combinations of
1 and
2 and then to minimize the obtained potential energy on curves of constant
=
1 -
2. This is the more realistic version of the approach shown above, yet it is still an approximation for an aggregate since the optimized combinations of
1 and
2 for a given
will not be possible with respect to all neighbors of a given DNA molecule. In Fig. 5 we compare this approach for N = 10 discrete charges at an interaxial separation R = 2.1 nm with the one presented above and with the continuous line charge version. We again have
= 0.9, f1 = 0.3, f2 = 0.7, and f3 = 0 for the counterion condensation parameters. The resulting potential curve is the lowest in energy, as one should expect from the procedure applied. The structure is close to the one induced by continuous line charge distributions and the minima are found at exactly the same loci as when keeping
1 fixed at
1 = 0; they are thus practically at the same positions as for the continuous line charge version. Again, from the analysis of one single interaxial separation we thus conjecture that the overall behavior of the pair potential will not present significant deviations from the reference continuous charge case. This statement is confirmed by analyzing the inset of Fig. 6, where we show the optimal azimuthal angles as a function of the interaxial separation R for the three different approaches to the charge distributions. Again, the dependence of the optimal azimuthal angle on the interaxial distance is very similar for the three cases studied, which will induce similar angular frustration behavior in an assembly. In detail, the optimal angle curve is closer to the one for a continuous charge distribution in the case where both angles
1 and
2 are free to rotate and the energetically optimal combination yielding the desired mutual azimuthal orientation
=
1 -
2 is chosen, as compared to the case where
1 is set to zero. As far as the behavior of the pair interaction at optimal azimuthal angle, shown in the main graph of Fig. 6, is concerned, both discrete charge versions fall on the same line, which shows a deviating course from the continuous version's behavior in the close-interaxial separation regime, whereas it approaches the continuous case's curve fast for larger R and both lines agree for R > 25 Å. We repeated the analysis of Figs. 5 and 6 for N = 20 charges. Here, no difference to the continuous charge distribution results could be discerned.
We can thus conclude that, first, the behavior of a DNA assembly will most probably not qualitatively differ for a discrete charge model with the real DNA charge number N = 10 and a continuous line charge model. The results will, however, in a quantitative manner depend on the underlying pair potential, especially for high concentrations, since, in Figs. 4 and 6 we found that for very close intersurface separations the pair interactions differed for a discrete and a continuous charge pattern on the DNA surface. Second, since already for N = 20 the results are indistinguishable from the ones for continuous line charges, we can furthermore surmise that a modulated continuous line charge distribution, as briefly discussed above to be the most realistic model, would not significantly differ even on the level of the pair potential. According to this reasoning, we will henceforth focus on continuous line charges, thereby avoiding the problem that for discrete charge patterns the potential depends on both molecules' azimuthal orientations
1 and
2, which significantly complicates matter for the strict analysis of an assembly.
Let us now investigate the effect of different amounts and types of counterions adsorbed on the DNA molecular surface. The type of counterion is herein modeled by the ratio of adsorbed charges in the minor and major grooves, as well as on the strands to the DNA phosphate backbone charge. We restrict our analysis to the most relevant cases: we will investigate
= 0.9 (meaning that 90% of the DNA charge is compensated by adsorbed counterions) with counterions adsorbing predominantly in the major groove, represented by charge fractions f1 = 0.3, f2 = 0.7, and f3 = 0, as well as with counterions exhibiting a high affinity to phosphates and thus condensing on the strands: f1 = 0, f2 = 0, and f3 = 1. A charge compensation value of
= 0.9 is known to be typical for DNA condensation (Kornyshev and Leikin, 1999
; Bloomfield, 1996
). Furthermore we calculate the potential for
= 0.7, which is a lower bound still occurring in DNA aggregation phenomena. Here, we also assume f1 = 0.3, f2 = 0.7, and f3 = 0. In Fig. 7 the potential is displayed as a function of the azimuthal angle
for two fixed interaxial separations, R = 2.5 nm and R = 3.0 nm, for f1 = 0.3, f2 = 0.7, f3 = 0, and
= 0.9 and
= 0.7. For both amounts of adsorbed counterions, the potential curves qualitatively agree. Due to the higher degree of charge compensation, however, the
= 0.9 potential values are smaller. In a subsequent step we minimize the potential with respect to the azimuthal alignment angle
, obtaining U(R,
opt). The result is displayed in Fig. 8. Both potentials being induced by situations where the majority of counterions condenses in the major groove are strongly attractive, with the one for
= 0.9 exceeding the one for
= 0.7. The potential stemming from a situation with all counterions condensed on strands, on the other hand, is purely repulsive.

View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 8 Yukawa-segment pair potential for two segments of length Lp as a function of the interaxial separation R of two DNA molecules, at the optimal angle opt(R), depicted for different values of the counterion condensation parameter and for different counterion adsorption patterns. The dependence of the optimal angle on the interaxial separation R is shown in the inset.
|
|
What is the origin of this qualitative difference? The mechanism can be thought of as a zipper (Kornyshev and Leikin, 1999
). Having a high charge compensation in the major groove creates a big charge separation: a negative helical line charge is located at the phosphate backbone position; a positive helical line charge rests in the adjacent major groove. With two opposing DNA molecules appropriately oriented, this allows for positive and negative charges to directly face each other as complementary parts in a zipper, creating a strong attraction between the two molecules. If counterion condensation solely occurs on strands, this mechanism is absent, creating a purely repulsive potential, as seen in Fig. 8 in the case of
= 0.9, f1 = 0, f2 = 0, and f3 = 1. In any case, the potential quickly decays toward zero for increasing interaxial separations so that in an assembly the dominant contributions to the total potential energy will stem from the nearest neighbors. The optimal angle, as a function of the interaxial separation, plotted in the inset of Fig. 8, is practically unaffected by this mechanism: in all three cases displayed, the optimal angle is nonzero for interaxial separations smaller than R*
28.25 Å, and zero else. Furthermore, a very similar increase from zero at R* to
opt(R = 20 Å)
0.47
is observed in all cases.
Let us finally remark that the Yukawa-segment model has the advantage of being very general and flexible. Any linearized field theory necessarily ends up with an effective Yukawa-type interaction. If hydration effects are included within a field theoretical description, the leading term for the effective interaction has again a Yukawa form. The electrostatic effects are well-described even at strong coupling provided the charges and screening lengths are suitably renormalized, as recently demonstrated in microion-resolved computer simulations of two parallel DNA strands (Allahyarov and Löwen, 2000
).
Kornyshev-Leikin potential
The Kornyshev-Leikin approach rewrites the result of linear screening theory in terms of a helical Fourier expansion (
1 <<
) (Kornyshev and Leikin, 1997
, 1998a
,b
). The pair interaction potential per unit length features a hard-core repulsion for interaxial separations R
2a, and for R > 2a reads as:
 | (9) |
The total interaction U(R,
) per segment of length L is simply U(R,
) = Lu(R,
). In Eq. 8
z denotes a vertical displacement, equivalent to the azimuthal alignment angle
= (2
/H)
z. Furthermore, u0 = 8
2/
2 (
2.9 kBT/Å at physiological ionic strength), and
. The function
n,m(x,y) is given by
 | (10) |
with the modified Bessel functions Kn(x) and Ij(y). The primes denote derivatives. As can be seen, the dependence of the pair potential on the mutual orientation angle
is affected by the distributions fi, i = 1,2,3 of the condensed counterions (Kornyshev and Leikin, 1999
). The dependence on the interaxial separation R is exponential. Keeping only the n = 0 term in the sum of Eq. 8 yields a pair potential of homogeneously charged cylinders, depending on R only. Summing up to |n| = 2 results in the approximation u(R,
)
C(R) - A(R)cos
+ B(R)cos 2
. Already at this level does the interaction potential u(R,
) show a peculiar dependence on the mutual azimuthal orientation angle, being a remarkable effect of DNA double-strandedness, as discussed above in the previous subsection. Here, A(R), B(R), and C(R) > 0 depend on the parameters of DNA structure as well as on the distribution of adsorbed ions, and A(R) > B(R) at large interaxial separations R. This potential has two symmetric azimuthal minima at
±
0 for distances smaller than a critical one at which A(R) = 2B(R), and one minimum at
= 0 for larger R. It thus already captures, to quite a good degree, the essential features of the full interaction potential as observed in the previous section in the framework of the Yukawa-segment model.
Let us now investigate the full potential. Due to rapid convergence of the sum in Eq. 8, truncation after the |n| = 5 terms suffices for the evaluation of the fully converged pair interaction potential. In Fig. 9 we show the KL potential U(R,
opt) at optimized azimuthal alignment angle,
opt, as plotted for the YS case in Fig. 8. It can be seen that the results are very similar to the ones discussed above for the YS potential. Counterion condensation on strands (f1 = 0, f2 = 0, and f3 = 1) gives rise to an exclusively repulsive potential, whereas condensation of a majority of the counterions in the major groove (f1 = 0.3, f2 = 0.7, and f3 = 0) results, at both charge compensations
= 0.7 and
= 0.9, in an attractive pair interaction. Differences in the KL approach to the YS approach can, however, also be inferred from a comparison of Fig. 8 and Fig. 9. These refer to the behavior at small intersurface separations. In the case of the YS model, the interaction decays monotonically to the contact value at R = 20 Å. Here, for the KL potential, however, the potential drops to its minimum value close to surface contact, but then rises again upon further approach. Furthermore a quantitative difference can be seen for
= 0.7 and f1 = 0.3, f2 = 0.7, f3 = 0, since in the KL case the attraction for this combination of parameters is much weaker than it was found to be for the YS model. Both observations can be attributed to the fact that in the KL case the dielectric jump is taken into account with
/
1 =
, where
is the dielectric constant of the solvent (water) and
1 is the dielectric constant of the DNA core. This allows for image charges at the DNA surface, bringing about a short-ranged repulsive part in the interaction, as evidenced in the potential curves in Fig. 9. Nonetheless, this short-range repulsion does not affect the behavior of the optimal angle as a function of the interaxial separation as compared to the one found in the YS case. We show the corresponding data in the inset of Fig. 9. The same functional form as for the YS potential is obtained, except for the fact that R*KL
29.5 Å is found to be slightly larger than R*YS
28.25 Å in the YS case.

View larger version (22K):
[in this window]
[in a new window]
|
FIGURE 9 Kornyshev-Leikin pair potential as a function of the interaxial separation R of two DNA molecules, at the optimal angle opt(R), depicted for different values of the counterion condensation parameter and for different counterion adsorption patterns. The dependence of the optimal angle on the interaxial separation R is shown in the inset.
|
|
We now have two realizations of the linear Debye-Hückel potential for the DNA interaction at hand stemming from different levels of modeling realized in the Debye-Hückel framework, which show differences with respect to the short-range behavior. In the following section we will present a theory to investigate the statistical properties of a columnar DNA assembly. We will thereby rely on the two YS and KL potentials as discussed above. The interesting question to be pursued apart from the main objective, being the general properties of such assemblies, is if and how the differences in the pair potentials affect the behavior of the assembly.
 |
A THEORY FOR DNA ASSEMBLIES
|
|---|
In the previous sections we showed that under the assumption of continuous line charges and infinitely long, rigid, parallel DNA molecules, the pair interaction potential U(R,
) of two DNA molecules only depends on the interaxial separation R and the mutual azimuthal orientation angle,
. The problem of statistical properties of columnar aggregates of long rigid DNA molecules may thus be mapped on a two-dimensional problem of particles that we may formally refer to as x-y spins, interacting via this unusual potential U(R,
) (see Fig. 2; see also Kornyshev and Leikin, 1997
). We repeat from previous subsections that the dominant contributions to the potential U(R,
) arise from the nearest neighbor interactions, as the R-dependent parts of the potential exponentially decrease with R. Before going into more detail on the theory for DNA assemblies we can, on the basis of the knowledge of the pair potential, already surmise a general trend in the behavior: we know that the potential has two symmetric azimuthal minima at
±
0 for distances smaller than a critical one and one minimum at
= 0 for larger R. Although the
= 0 case is compatible with any lattice,
0 results in frustrations of positional and orientational order (Strey et al., 2000
). Due to the R-
coupling in the interaction potential, one may expect peculiar positional and orientational structures in the aggregate, a feature known as the mesomorphism of DNA assemblies (Podgornik et al., 1998
). Carrying the formal analogy to spin systems further, we may refer to the orientational structure in the assembly also as spin or magnetic structures.
Lattice sums
For all cases studied in this article, the pair interaction U(R,
) is greater than kBT, so that the energy needed to destroy the translational or orientational order in an assembly must be more than several kBT at room temperature. Hence focusing on the ground state analysis of the basic structures of the assembly provides the representative thermodynamic states. This reasoning is further sustained by evidence from polymer crystallization, stating that upon compression the effective persistence length (this persistence length has to be distinguished from a smaller correlation scale which is decreasing due to deflections of the polymer within the tube; see Vroege and Lekkerkerker, 1992
) of polymers increases, bringing them into columnar alignment at high packing fractions. Since, as we already argued above, the problem is effectively two-dimensional, we consider the five two-dimensional Bravais lattices, i.e., the hexagonal (HEX), square (SQ), rectangular (REC), rhombic (RHO), and oblique (OBL) lattices to assess the representative thermodynamic states. As for the exploration of the ordered spin structures, we are, in principle, facing infinitely many degrees of freedom: every DNA molecule in the lattice has a continuous spectrum of possible orientations. We can, however, make use of a pair potential property that we noted in The Pair Potential, namely, that the pair interaction drops exponentially as a function of the interaxial separation R. Assuming that its range would solely encompass interactions contained in a fundamental unit cell (elementary plaquette), the approach could be much simplified in the following way. We restrict our analysis to finding the minimal energy state of this fundamental unit cell alone. This is achieved by minimizing the energy of the plaquette with respect to all spin angles residing on the elementary plaquette. Since no interactions beyond unit cells are assumed to be present, periodical repetition of this minimized unit cell along the lattice directions guarantees to give the ground state of the whole lattice. Due to the exponential decay of the R-dependent factors in the pair interaction potential this already presents an amazingly good approximation for our purposes. Since strictly speaking the range of the potential may extend beyond nearest neighbor interactions in some cases, we adopt a perturbation approach in the following way: the whole lattice is generated by periodical repetition of the elementary plaquette structure, involving two or three degrees of freedom depending on the lattice type under exploration, but interactions of higher order neighbors are nonetheless included in the calculation of the lattice sums.
The algorithms employed for generating the ordered spin structures on the whole lattice building on the fundamental unit cell differ depending on the lattice type. They are schematically illustrated in Fig. 10. One of the spins in the elementary plaquette is chosen as reference (
= 0). This leaves two degrees of freedom (
1,
2) in the case of the HEX lattice, and three degrees of freedom (
1,
2, and
3) for the REC and SQ lattices. The HEX lattice can be build up by periodically reflecting the unit cell across its edges, as is shown in Fig. 10 a. The same holds for the REC lattice with three free orientations per plaquette; see Fig. 10 b. In the case of the RHO and OBL lattices, however, employing the same procedure as for the REC lattice with three free spin angles per plaquette does not produce identical plaquettes: due to the fact that the geometrical symmetry of the unit cell is broken (a short and a long diagonal exist), mirror reflections across the edges generate different plaquettes on the whole lattice. The lattice may nonetheless be filled with identical plaquettes by employing two algorithms which are depicted in Fig. 10, c and d. In the first, spins of orientation
1 and
2 are placed along the edges, whereas the third free orientation angle is chosen to be
3 =
2
1. The whole lattice is then populated by successive mirror reflections ensuring that pairs of spins across all diagonals have the same relative angle of
2
1.

View larger version (17K):
[in this window]
[in a new window]
|
FIGURE 10 A schematic view of generating candidate ordered spin phases of the system. (a) for the HEX lattice; (b) for the REC and SQ lattices; and (c) and (d) for the RHO and OBL lattices.
|
|
The second algorithm is illustrated in Fig. 10 d. Again, spin angles of values
1 and
2 are chosen along the edges, whereas
3 is assigned a value of
1 +
2. The lattice is then built up by increasing the angular value by
2 along the oblique direction and by
1 along the horizontal lattice direction. The resulting lattice exhibits unit cells in which all pairs of spins across short diagonals have an angle difference of
1
2, whereas pairs of spins across long diagonals are separated by an angular difference of
1 +
2.
Lattice sums are then calculated and a minimization of the lattice energy with respect to the orientational degrees of freedom {
i}, the geometrical degrees of freedom (being the size ratios b/c for the REC lattice and/or the geometrical angle
for RHO and OBL lattices; see Fig. 10), is carried out. The result of the minimization procedure is the optimized lattice-sum energy, UX(
,
), where X stands for the lattice type, and
= (
1,
2,...,
N), denotes the configuration of the N spins in the system.
Three examples of lattice sums for fixed DNA density 
a2 and fixed salt concentration ns at a charge compensation
= 0.9 and f1 = 0, f2 = 0, and f3 = 1 are displayed in Fig. 11, ac. They depict the total energy stemming from the lattice sum as contour plots as a function of
1 and
2. They are representative of three different phases emerging for these parameters. The meaning of the three phases will be explained in detail in the next section. It can be clearly discerned from the contour plots that a certain symmetry prevails in the aggregate with respect to
1 and
2 whereby the symmetry axis is the line
1 =
2. The location of the minima evolves from 