Order and interactions in DNA arrays: Multiscale molecular dynamics simulation

While densely packed DNA arrays are known to exhibit hexagonal and orthorhombic local packings, the detailed mechanism governing the associated phase transition remains rather elusive. Furthermore, at high densities the atomistic resolution is paramount to properly account for fine details, encompassing the DNA molecular order, the contingent ordering of counterions and the induced molecular ordering of the bathing solvent, bringing together electrostatic, steric, thermal and direct hydrogen-bonding interactions, resulting in the observed osmotic equation of state. We perform a multiscale simulation of dense DNA arrays by enclosing a set of 16 atomistically resolved DNA molecules within a semi-permeable membrane, allowing the passage of water and salt ions, and thus mimicking the behavior of DNA arrays subjected to external osmotic stress in a bathing solution of monovalent salt and multivalent counterions. By varying the DNA density, local packing symmetry, and counterion type, we obtain osmotic equation of state together with the hexagonal-orthorhombic phase transition, and full structural characterization of the DNA subphase in terms of its positional and angular orientational fluctuations, counterion distributions, and the solvent local dielectric response profile with its order parameters that allow us to identify the hydration force as the primary interaction mechanism at high DNA densities.


I. Interaction phenomenology
Monovalent salts, apart from their specific interaction with the DNA surface 1,2 , provide usually a screening electrolyte medium for electrostatic repulsion between highly charged dsDNAs, while trivalent cations such as CoHex 3+ and/or Spd 3+ at sufficiently large concentrations can lead to dsDNA condensation by turning electrostatic repulsion into a counterion mediated attraction, despite DNA's high negative charge [3][4][5][6][7] , enforcing a phase separation from the bathing electrolyte.
While high valency is a necessary condition for attractive interaction, it is not exhaustive as e.g.
CoHex is more efficient as a condensing agent then spermidine, though they are both equally charged 1 . On the other hand, monovalent cations such as K + can in their turn act as condensing agents for quadruple stranded DNA 8 ; furthermore, whereas some divalent ions such as Cd 2+ and Mn 2+ can under certain circumstances condense the double stranded DNA, other divalent salts, containing e.g. Mg 2+ cation, again act plainly as a screening electrolyte medium for dsDNA, but at the same time invariably condensing the triple stranded DNA 9 . These features make the unique identification of the ultimate mechanism of counterion-mediated DNA condensation difficult to pin down.
The phenomenology of DNA-DNA interactions, as revealed by the EoS experiments, is complicated and defies a single classification scheme, depending on a multiplicity of conditions. While it is clear that purely mean-field electrostatic theories based on the Poisson-Boltzmann phenomenology cannot predict anything but repulsion 5 , venturing beyond this framework can be pursued on different levels. By concentrating on the non-mean-field correlation effects, as first suggested by Rouzina and Bloomfield 10 , the Wigner crystal approach 11 and the strong coupling approach 12 as well as their varieties, such as the charge density wave 13 and ion bridging model 14,15 , all single out the universal features of higher valency counterions, based on their valency alone, as responsible for the attractions that they mediate. On the other hand, the tight counterion binding approach 1 was the first to introduce non-electrostatic interactions as governing the distribution of the counterions on the surface of DNA, while preserving the electrostatics as the mechanism mediating the interactions between thus decorated DNAs. The furthest deviation from the predominant role of electrostatics was the proposition by Rau and Parsegian 16 that water-mediated interactions present the most important source of DNA attractions and that the way to understand them is to abolish the implicit dielectric continuum model for the solvent.
II. Additional simulation details DNA molecules are always described on the atomistic level and modelled with the AMBER 03 force field 17 . On the atomistic scale, the water molecules are modelled with a TIP3P model 18 .
For the Na + and Cl − ions we first test the standard AMBER 03 parameters, but find that the DNA molecules substantially deviate from their reference lattice positions, in accordance with Ref. 19 . There, DNA molecules under similar conditions were observed to form clusters, contrary to the experimental observations. Furthermore, discrepancies were also found for the osmotic pressure of the DNA arrays. Thus, we follow the guidance of Ref. 19 and instead use the Joung and Cheatham parameters 20 with additional corrections to Na + -phosphate interaction. The Spd 3+ molecules are modeled only on the atomistic level (not leaving the atomistic domain) with the parameters obtained with AmberTools program 21 , employing a standard procedure for obtaining custom molecule force field parameters. The water molecules and ions (Na + and Cl − ) can change their resolution adaptively on the fly between the atomistic and coarse-grained regions. The effective interactions on the coarse-grained scale are determined via Iterative Boltzmann Inversion method 22 that is incorporated into the STOCK web toolkit 23 .
For each lattice type and initial DNA-DNA separation the DNA molecule is translated and rotated to form an array. In the hexagonal configuration, the DNA molecules are rotated by 2π/3 with respect to each other. In the orthorhombic configuration, the four nearest neighbors are rotated by π/2, whereas the two further away are at the non-optimal orientation ( ϕ = 0). The array is is applied also to the DNA helix. In particular, we use patches, i.e., corresponding intramolecular DNA interactions defined by bond, angle, and dihedral interaction potentials, to connect each strand to its periodic image along the z-axis. While it is more customary to use the experimentally determined DNA structures for the initial structure, we opt here for a generic one. In particular, the initial atomic coordinates of a B-form DNA with 10 base-pair sequence (5'-CTCTCGCGCG-3') molecule are generated using the 3D-DART Web server 24 . With this choice we avoid the possible strains that could occur at the periodic links of both chains. The atomistic region is a rhombic prism and for the hexagonal latices coincides with the semipermeable membrane. For the orthorhombic lattices the size is set to encompass the whole membrane, which is tilted at a different angle (see Figure S1). The width of hybrid and coarse-grained regions are set to 0.9 nm for all cases. The resolution region boundaries are static during the simulation run. The thermodynamic forces applied in the adaptive resolution scheme (AdResS) are depicted in Figure S2 for water, Na + and Cl − . They are slightly modified form Ref. 25 due to different model parameters.

III. Osmotic pressure calculation
The osmotic stress methodology to determine the osmotic pressure in the DNA subphase relies on the equivalence between the mechanical and the osmotic forces acting on a semipermeable membrane, the latter controlled by a water soluble polymer (usually polyethyleneglycol -PEG) 26 , that exerts a uniform compression on the DNA subphase 16 .
In MD simulations, the osmotic pressure between pure water and salt solution region can be calculated by introducing a semi-permeable membrane 27,28 . We calculate the osmotic pressure of a DNA array in a similar way by placing a semi-permeable membrane around the set of DNA molecules. The membrane is modelled by a repulsive part of the 10-4 Lennard-Jones interaction of the following form where x is a perpendicular distance to the wall, σ and are the 12-6 Lenard-Jones parameters and ρ is the wall density. The interaction is applied to the DNA backbone atoms and thus prevents a given DNA molecule from crossing the wall, whereas the solvent molecules (water, Na + , and Cl − ions) can freely pass through. The system inside walls is thus open, i.e., semi grand canonical, even though the total system is closed. The osmotic pressure of the DNA array is then calculated as F wall /A, where F wall is the mean force on the wall and A is its area.
In DNA arrays of hexagonal symmetry, the osmotic pressure can be equivalently plotted as a function of the interaxial spacings between nearest neighbors or the DNA density, while in the orthorhombic phase the DNA-DNA spacings are not all equal and thus the osmotic pressuredensity plot offers a better comparison between the two lattices. Note that the DNA density is defined as an average inside the semi-permeable membrane and depends on its size. Additionally, the DNA-DNA separation varies during the simulation and deviates (especially for low densities) from the initially set value.

Size variation
The experimentally observed DNA liquid crystals contain a large number of DNA molecules, while in the simulations we can only simulate a few. Thus, it is important to first evaluate the size effects on the osmotic pressure. We tested four sizes of DNA assemblies that contain 9, 16, 25, and 36 DNA molecules. In Figure S3,

Determination of phase transition via Lindemann criterion
The translational and orientational orders of the DNA assemblies are also quantified by the normalized root-mean-square deviation σ r /a and the rotational root mean square deviations σ ϕ /π as a function of pressure, respectively (see the main article). The results are depicted in Figure S4.

Local pressure / Hydration force
The p yy component of the local pressure tensor within the membrane walls is calculated along the y-axis (see Figure S6) using the Method of Planes (MOP) definition 29,30 . The profiles of the pressure are calculated by splitting the axis into slabs. The slab component is composed of a kinetic contribution p kin yy (y) and a virial part p pot yy (y) defined as and The . . . denotes the configurational average and the m i , p i and f i are the mass, velocity and total force on particle i, respectively. Figure S5 shows the difference in pressure contributions of the local pressure tensor within the membrane walls is computed by averaging over y. Figure S7 showsp the total osmotic pressure and the pure DNA contribution, that can even become negative, in complete analogy with the case of interacting hydrated phospholipid membranes 31 . This difference indicates that water contributes fundamentally to the force equilibrium in the DNA array and the identification of the hydration force as its major component is thus vindicated.

Order parameters
The water order parameters are shown in Figures S8-S10. In the order parameter η (3) = Q ij = l q l r il r jl , the index l runs over the water hydrogens with charge q l and r is the distance between the hydrogen and oxygen atom of the considered water molecule. . . . denotes the average over trajectory and over all the DNA pairs.
Water in the DNA subphase tends to limit to its bulk properties for concentrations below ∼ 300 mg ml −1 , and deviates significantly from the full bulk properties at concentrations above ∼ 600 mg ml −1 .  Figure S8. FIG. S10. Water order parameters η (3) between pairs of DNA molecules. The representation is the same as in Figure S8.

Tetrahedral order
The tetrahedrality order parameter Q 4 is shown in Figure S11 and is defined as 32 where i is the oxygen atom of the reference water molecule and θ ijk the angle between vectors r ij and r ik , where j and k are the oxygen atoms of nearest neighbors of central water i. The sum runs over distinct pairs of the four closest neighbors, i.e., over six oxygen-oxygen-oxygen angles. The spatially varying Q 4 is computed by discretizing distances from the DNA molecule into bins and averaging over water molecules that fall into a corresponding bin.  Relative permittivity of water Figure S12 shows the relative permittivities of water as defined in Ref. 25 .

Occupancy and residential times of solvent entities
We compute the average occupancy and residence times of counterions Na + and Spd 3+ and oxygen atoms of water as detailed in Ref. 25  FIG. S12. Relative permittivities of water. The representation is the same as in Figure S8. FIG. S14. Same as Figure S13 for the orthorhombic lattices. FIG. S15. Same as Figure S13 for the hexagonal lattices with Na + and Spd 3+ counterions. FIG. S16. Same as Figure S13 for the orthorhombic lattices with Na + and Spd 3+ counterions.