Molecular Dynamics Simulation of High Density DNA Arrays

Densely packed DNA arrays exhibit hexagonal and orthorhombic local packings, as well as a weakly first order transition between them. While we have some understanding of the interactions between DNA molecules in aqueous ionic solutions, the structural details of its ordered phases and the mechanism governing the respective phase transitions between them remains less well understood. Since at high DNA densities, i.e., small interaxial spacings, one can neither neglect the atomic details of the interacting macromolecular surfaces nor the atomic details of the intervening ionic solution, the atomistic resolution is a sine qua non to properly describe and analyze the interactions between DNA molecules. In fact, in order to properly understand the details of the observed osmotic equation of state, one needs to implement multiple levels of organization, spanning the range from the molecular order of DNA itself, the possible ordering of counterions, and then all the way to the induced molecular ordering of the aqueous solvent, all coupled together by electrostatic, steric, thermal and direct hydrogen-bonding interactions. Multiscale simulations therefore appear as singularly suited to connect the microscopic details of this system with its macroscopic thermodynamic behavior. We review the details of the simulation of dense atomistically resolved DNA arrays with different packing symmetries and the ensuing osmotic equation of state obtained by enclosing a DNA array in a monovalent salt and multivalent (spermidine) counterions within a solvent permeable membrane, mimicking the behavior of DNA arrays subjected to external osmotic stress. By varying the DNA density, the local packing symmetry, and the counterion type, we are able to analyze the osmotic equation of state together with the full structural characterization of the DNA subphase, the counterion distribution and the solvent structural order in terms of its different order parameters and consequently identify the most important contribution to the DNA-DNA interactions at high DNA densities.


Introduction
The properties of dilute (double-stranded) dsDNA solutions do not really suffice in order to understand their behavior in vivo [1,2].The internal environment of eukaryotic nuclei, bacterial nucleoids and/or viral capsids, is anything but dilute and the properties of DNA solutions in these environments can only be characterized by specifically invoking the fine details of molecular interactions between them at small intersurface spacings [3].The most direct method to probe the molecular interactions between DNA molecules at high densities is the osmotic stress method as introduced by Parsegian and coworkers [4,5], that yields the thermodynamic osmotic equation of state (EoS), connecting the DNA concentration, as measured e.g., by SAXS, and the DNA osmotic pressure.The detailed dependence of the EoS on the concentration of DNA then allows one to deconvolute the fundamental interactions between DNA molecular surfaces as a function of their separation and different bathing solution conditions, such as the ionic strength of the bathing salt solution, and the identity of the counterions [6].For monovalent salt solutions, the EoS exhibits a monotonic behavior as a function of the concentration of DNA, characterized by the hydration and fluctuation enhanced electrostatic regimes [3,7,8].It was only recently, through high resolution control of the osmotic pressure based on its known temperature variation [9,10], that a small discontinuous DNA density change became discernible, between the cholesteric and hexatic phases and between the hexatic and orthorhombic phases for long fragment DNA [10].The exact sequence and identity of the phases actually depends on the DNA length [11,12] and is apparent also from a sudden discontinuous change in the corresponding order parameters [13,14].The detailed sequence of mesophase ordering on increase of the DNA density has been identified as consisting of isotropic (I) → cholesteric (Ch) → line hexatic (LH) (or hexagonal columnar H) → orthorhombic (OR) phase [3,11,[15][16][17][18], with many details, including the complete DNA length dependence and the demarcation between the line hexatic and hexagonal order, still remaining to be systematically investigated [19].For multivalent counterions at concentrations above a critical value depending on their identity [20], the EoS exhibits pronounced van der Waals-like density discontinuities that are in general larger then in the case of monovalent salt solutions, signalling a buildup of attractive interactions [21,22], whose details can be inferred from complementary experiments [23], that eventually lead to DNA condensation [24].
The effect of monovalent and most of divalent salts on the EoS can be usually rationalized in terms of a screening electrolyte medium for electrostatic repulsion between highly charged dsDNAs, modified by an additional specific interaction with the DNA surface [25,26].The only exception to the screening effect are Mn 2+ ions at an elevated temperature that actually behave much more like multivalent ions and their effect on the interaction cannot be rationalized within the screening paradigm [21,22].Counter to the screening effects, the trivalent cations such as CoHex 3+ and/or Spd 3+ at sufficiently large concentrations can lead to dsDNA condensation by turning electrostatic repulsion into either a counterion mediated or solvent mediated attraction, despite DNA's high negative charge [21,22,24,[27][28][29][30], thus enforcing a phase separation from the bathing electrolyte.While multivalent counterions are a necessary condition for attractive interaction, it is not exhaustive as, e.g., CoHex is more efficient as a condensing agent then spermidine, despite both of them being equally charged [25].On the other hand, monovalent cations such as K + can act as condensing agents for quadruple stranded DNA [31], while some divalent ions such as Cd 2+ and Mn 2+ condense the dsDNA, other divalent salts, e.g., Mg 2+ cation, act plainly just as a screening electrolyte medium for dsDNA, but as a condensing agent for the triple stranded DNA [32].While it is generally agreed that purely mean-field electrostatic theories based on the Poisson-Boltzmann phenomenology [33] cannot predict anything but repulsion [28], venturing beyond this framework can be pursued on different levels.By concentrating on the non-mean-field correlation effects [34], the Wigner crystal approach [35] and the strong coupling approach [36] as well as their varieties, such as the charge density wave [37] and ion bridging model [38,39], all single out the universal features of higher valency counterions as responsible for the attractions that they mediate.On the other hand, the tight counterion binding approach [25] 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 [21,22] 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.

Simulating DNA Arrays
Since direct experiments probing the interactions and/or the associated molecular order can provide only a limited insight into the structural details that would allow an unequivocal deconvolution of the main physical mechanisms underlying these interactions, theoretical interpolations are therefore necessary [40].While the continuum approaches which mostly underpin the theoretical endeavors described in the paragraph above led to important insights, e.g., the significance of counterion correlations in reversing the sign of the electrostatic interactions between charged helices, backed up by detailed coarse grained simulations [41][42][43][44][45] the molecular details of the DNA solution exposed surface, the granularity of the molecular solvent, and the complicated interactions between both and the mobile charges in solution, preclude the understanding of all the relevant details.As a consequence, detailed all-atom molecular dynamics (MD) simulations appear to be the only vehicle that can bring forth a deeper understanding of the mechanisms and the relevant couplings between them [46], leading to the experimentally observed interaction and ordering phenomenology of DNA in high density mesophases [47][48][49][50][51]. On the downside, all-atom MD simulations of such complicated systems require huge computational resources and the majority of simulations aiming to describe details of the DNA solution phenomenology consist of 1-3 DNA molecules, where the focus is the counterion binding [48,52] and azimuthal dependence of DNA-DNA interaction [49].Only recently it became feasible to set up an all-atom MD simulation for a larger set of DNA molecules [51,53] characterized by a single packing geometry with only a partial characterization of the DNA countercharge and solvent ordering.This approach has been later extended by the applications of the multiscale MD technique AdResS (Adaptive Resolution Scheme) [54][55][56][57][58][59][60][61][62][63][64][65][66][67], which has been already successfully applied to various biological systems [68][69][70][71][72][73][74][75][76][77], enabling a concurrent and consistent coupling between the atomistic (AT) and the coarse-grained (CG) representations with a key feature of allowing molecules to freely move not only in real space but also in the resolution space across different regions and change their resolution on the fly according to their position in the computational domain.

Modeling the Solvent Grand Canonical (Osmotic Isobaric) Ensemble
In order to capture the salient properties of high density DNA arrays and their corresponding EoS, any simulation should address the problem of defining and maintaining the osmotic isobaric ensemble, i.e., to maintain the simulated system at a set osmotic pressure, or equivalently at a fixed chemical potential of water, equivalent to the solvent grand canonical ensemble with exchangeable water and solution ions with bulk reservoirs.Fixing the chemical potential of water can be accomplished generally in different ways.In the context of hydrated multilamellar lipid membrane assemblies, the constant water chemical potential is maintained either by including a large water reservoir [78], by direct grand-canonical simulation [79,80], or by the recently developed thermodynamic extrapolation method [81].The simulation setup corresponds either to two parallel planar lipid monolayers facing the water phase in the middle [82], or a single periodically replicated planar lipid bilayer in the multilamellar liquid-crystalline phase [83].The extrapolation method then allows for efficient simulations at prescribed chemical potential and yields the interaction forces between membranes with high precision.This enables also a quantitative comparison with measured hydration interactions between two overall charge-neutral lipid membranes in water at small separations and mimics the interactions in a densely packed phospholipid bilayer stack where the undulation forces are suppressed [84,85].While this type of modeling can be pursued for a one dimensional array of phospholipid membranes, it is nevertheless inapplicable for a two dimensional array of dsDNA molecules.In fact, contrary to the planar membrane case, the interactions between dsDNAs have a strong orientational component [16,25] that fundamentally affects the collective packing symmetry of the ordered array [13,17].Simulating just two dsDNA molecules in a bathing solution would entirely miss the ensuing orientational frustrations which lead to collective, non-pairwise additive effects.In addition to an implementation of the osmotic isobaric ensemble, simulation of a full array composed of as large a number of DNAs as possible is, therefore, a requirement.Applying the large water reservoir simulation methodology in order to implement the osmotic isobaric ensemble, the DNA subphase can then be segregated from the large water reservoir by way of either a cylindrical [51] or a combined orthorhombic-hexagonal semipermeable boundary [73].Changing the symmetry of the boundary confining the DNA subphase is a convinient approach to capture the phase transitions with a change of packing symmetry, such as the cholesteric (Ch) → line hexatic (LH) (or hexagonal columnar H) and/or the line hexatic → orthorhombic (OR) phase transitions [3,11,15,18].

Bathing Solution
Since the effects of monovalent and most of the divalent salts on the EoS, rationalized in terms of a generic screening electrolyte medium, is very different from the effects of multivalent cations, one needs to simulate at least two cases of the bathing solution composition: pure monovalent or divalent salt such as NaCl and MgCl 2 at a fixed concentration, and a mixed solution composed of a monvalent salt background in addition to CoHex 3+ , Spd 3+ or Sm 4+ counterions that are known to exhibit strong electrostatic correlation attractions between equally charged apposed DNA surfaces.In the two most recent simulations, the parametrizations were chosen such as to be close to experimental conditions in the osmotic stress experiments [4,23], with either 250 mM or 1 M pure NaCl [73], or with 200/50 mM NaCl-MgCl 2 ionic mixture [51], leading in both cases to the pronounced screening of the electrostatic interactions.On the other hand, the strongly coupled ions were simulated either as a combination of 1 M NaCl and spermidine Spd 3+ , where 6 Spd 3+ molecules are used per each DNA molecule [73] or sub-mM concentration of Sm 4+ and 200 mM NaCl [51].In both cases, the system was neutralized by pure Na counterions or a combination of pure Na counterions and multivalent cations.The monovalent salt concentrations correspond to Debye screening length of κ = 3.3 nm −1 for 1 M NaCl and to κ = 0.68 nm −1 for 200 mM NaCl.
The pertaining EoS phenomenology was investigated for a finite orientationally ordered array of 64 vs. 16 explicit atom-resolved DNAs within a solution-permeable membrane, modeled as a confining cylinder [51] or as a confining polygon allowing for hexagonal and/or orthorhombic packing symmetries [73].In the latter case, the modeling of the semipermeable membrane allows for comparison of the osmotic pressure for two high density phases with different packing symmetries.Characterization of the positional and angular orientational order of the DNA sub-phase together with the Lindemann criterion allows one to extract the osmotic pressure of the phase transitions between the phases of different order and symmetry, all concurrently with the distribution of counterions and positional, orientational and local tetrahedral coordination order of the water molecules.

Model Interaction Potential Parametrization in DNA Arrays
One of the challenges regarding the molecular simulations of DNA molecules is the reliability of the force fields.On the atomistic scale, simulations are most commonly performed employing the GROMOS [86], AMBER [87], CHARMM [88] or OPLS [89] force fields.The comparison of the two for the Dickerson-Drew dodecamer showed that the AMBER force field was more stable [90].On the CG level, the DNA molecule is typically modeled with three sites per nucleotide [50,[91][92][93][94].The interactions between CG beads can be obtained either with a top-down approach, i.e., fitted to reproduce key experimental data, especially thermodynamic properties, or with the bottom-up approach, employing methods such as inverse Monte Carlo [95], iterative Boltzmann inversion (IBI) [96], force matching [97][98][99], and relative entropy [100][101][102].Recently, also the popular MARTINI force field [103] was extended to DNA [104].On an even coarser level, the DNA can be modeled as a cylinder with charge located on the sites corresponding to the phosphate group in the B-form of DNA [41,42,105].These models were primarily used early on in Monte Carlo simulations of DNA arrays and were able to reproduce some essential features of the counterion-mediated DNA attraction.As already stated above, the DNA systems are highly sensitive to the bathing solution.Therefore, the calibration of the force field is not only relevant for the DNA molecule but also for the water and salt of the bathing solution.In a recent paper by Yoo et al. [53], the standard AMBER 03 parameters for the Na + and Cl − ions were actually found ill parametrized for simulations of DNA arrays.In particular, the DNA molecules deviated from their initial lattice positions and formed clusters, which is in disagreement with the experimental results.Better agreement with the experiments was found for the Joung and Cheatham salt parameters [106], chosen such as to reproduce the relevant experimental osmotic pressure of DNA arrays [53].In the context of AdResS, where any two AT and CG salt solution models can be coupled, the structure-based coarse-graining approaches, i.e., the Iterative Boltzmann Inversion method [96] that is incorporated into the STOCK web toolkit [68,107], are more convenient, as they allow a direct structural comparison with the underlying AT model.

Structural Characterization of the High Density DNA Subphase
Extensive MD simulations of DNA arrays with different local packing symmetries allow for a detailed structural characterization of the molecular order in, e.g., the H and OR phases and to identify the driving mechanism(s) of the phase transition between them.The EoS is obtained from multiple simulations with varied DNA-DNA densities corresponding to different average interaxial spacings for each bathing solution composition, either at lower densities [51] or high DNA densities [73] that are relevant for the H-OR phase transition.To conduct the simulations efficiently, the multiscale AdResS was adapted to the simulated system [73]; see the next sections.
The phase transition between the H and OR phases was identified and characterized by invoking the general Lindemann criterion as a proxy for the transition point [73].Resorting to the Lindemann criterion is an additional approximation and could be avoided if the phase transition were determined more precisely and reliably from the computed free-energy landscape bypassed by enhanced sampling techniques [108][109][110][111].However, this approach would require much more computational time and computer resources.The computed EoS in the high density DNA density regime in the presence of mono (Na) and multivalent (spermidine) counterions differ typically in the sense that the later exhibits a smaller osmotic pressure at the same background salt and DNA conditions, indicating the existence of an effective attractive interaction, possibly related to the bridging configurations of the multivalent ion between the vicinal DNA surfaces.In addition to the positional and orientational order of DNA molecules, the related local order parameters of water, its local dielectric response profile and the occupancy and the residence times of the ions and water can be obtained from the same simulations.The behavior of the local order parameters of water are then seen to be consistent with a predominant contribution of hydration forces to the DNA EoS at high densities as has been argued all along on purely experimental grounds by Parsegian, Rau and coworkers [26].

Adaptive Resolution Simulations of a DNA Molecule Solvated in Salt Solution
In previous sections, we have underlined the importance of performing DNA simulations in the solvent grand canonical ensemble that allows for the exchange of ions and water with the environment.Here, we start to describe our efforts to achieve this goal by employing AdResS, most straightforwardly implementable for a single DNA molecule embedded in a multiscale salt solution [70,71].We focus on (spatially varying) dielectric properties of solvent and the DNA molecule itself.In simulations, implementing the salt solution environment of DNA demands vast numbers of ions and associated large simulation box sizes.Hence, full-blown AT simulations of such systems are computationally demanding.A possible way to deal with this issue is to use CG models, while at the same time one has to keep in mind the fact, that the dielectric and solvent density profiles display bulk behavior only beyond the second solvent coordination shell of the DNA molecule [70].Therefore, it would be computationally advantageous to use a simplified CG representation for the salt solution in the "bulk" domain, whereas resorting to the AT level of detail for the proximal solvent.Multiscale simulations that concurrently couple the AT and CG resolutions provide the most efficient way to deal with the implementation of this varying simulation resolution, as they circumvent the limitations of both the AT and CG simulations by providing not only computational benefits but also enhanced physical insight [55].

AdResS-Single DNA
AdResS [54,55] seamlessly couples two or more levels of detail in the system.In AdResS, recapitulated below, the different resolutions are interlinked via a force interpolation with the total force F α on a given molecule α (or a bundle in the case where a supramolecular mapping is used) defined as The continuous weighting function w, which governs the transition between different resolution domains, depends on the position of a given molecule (bundle) in the simulation box and is defined such that w = 1 corresponds to the AT, w = 0 to the CG, and 0 < w < 1 to the hybrid (HY) domains, respectively.The HY region is sandwiched in between the AT and CG regions and enables an AT molecule to gradually find its energetically favorable orientation with respect to its neighboring molecules.A local thermostat, e.g., the Langevin or Dissipative Particle Dynamics (DPD) thermostat [112], is employed to supply or remove the free-energy associated with the switched on/off molecular degrees of freedom in the HY domain.AdResS allows for various geometric boundaries between the resolution regions.These boundaries can be, for example, one directional [113], cylindrical [70], spherical [69], or even self-adjusting [75].The relative vector R α − R is then a one, two or three dimensional vector, respectively, from the molecular position to the center of the AT domain R. The AT region can be centered at either a fixed or a mobile point, as for instance, in a simulation of a protein where it coincides with the protein's center of mass (CoM).Such setup ensures that the macromolecule and the proximal solvent are always simulated at the AT resolution.
AT and CG molecular models can differ substantially in their thermodynamic properties such as pressure and chemical potential.Pressure (chemical potential) differences at the resolution boundaries lead to density undulations across the HY domain.To mitigate this artefact external thermodynamic (TD) forces F TD α opposing the chemical potential gradient are applied to water molecules (bundles) and ions in the HY region [56,114,115].Employing the TD force, the AdResS methodology can be applied to couple, in principle, any AT and CG models.The structure-based CG approaches (see Section 2.3) are routinely used as they allow a direct structural comparison with the underlying AT model.On the other hand, top-down approaches such as MARTINI [103] focusing on thermodynamic properties can also be applied.However, due to the reduction of degrees of freedom associated with coarse graining, the CG interactions are softer leading to reduced viscosity [55,57,115].As a result, the CG models frequently display accelerated dynamics.In some cases, the speedup in dynamics is advantageous and can facilitate a more efficient equilibration of complex systems [116].If required, the transport properties in the CG and AT parts of the system can be matched.This can be accomplished using local thermostats [117,118] or via the Mori-Zwanzig projection formalism [119][120][121][122][123][124] incorporating memory effects in the generalized Langevin equation.CG potentials are also density and temperature dependent [55,[125][126][127].In the context of biomolecular simulations employing AdResS, this does not represent a major drawback as the simulations are manly performed at a fixed temperature of 300 K. Furthermore, the CG region is considered as an external buffer to the region of interest, i.e., the AT region.
We have demonstrated that using the above AdResS approach structural properties of a simulated DNA molecule and surrounding salt solution, presented in Figure 1, are faithfully captured in the region of interest [70].One can even increase the level of coarse-graining (and related computational efficiency) for the CG water and use a popular water model as the MARTINI [103], in which one CG bead corresponds to four water molecules [69,71,113], see Figure 2.

Dielectric Properties
It is known and also confirmed by our simulations that a DNA molecule substantionally influences the surrounding aqueous solvent, i.e., water is considerably ordered and has a slower dynamics owing to the motional restrictions caused by the DNA surface molecular groups and the solution ions.The dielectric constant is hence expected to be much lower than in the bulk [70].To validate this, we computed the spatially varying dielectric constant of water around a single DNA molecule, employing the Kirkwood's theory, in which the dielectric constant is related to the average vector sum of the dipole moments of a water molecule centered in a spherical free volume S inside the solvent continuum, with RF = 80, i.e., [70,128], µ is the external dipole moment of the water molecule, k B is the Boltzmann constant, T the temperature, and V ex,i the excluded volume due to DNA atoms and ions, respectively.The dipole density ρ is calculated as N i /(V − V ex,i ) i .Here, V is the volume of a Kirkwood sphere S, which is centered on a given water molecule i, with N i the number of water molecules with indices j = 1, N i in the Kirkwood sphere S around the reference water molecule i. ... i denotes an average over all water molecules in the set, and θ ij the angle between dipoles of water molecules i and j.
The spatially varying profile of is computed by discretizing distances from the DNA molecule into bins and calculating the average ... i , in Equation ( 2), over water molecules that belong to a particular bin.We found that the dielectric constant is increasing with increased distance from the DNA's CoM and reaches the bulk value at ∼1.7 nm.We also studied the dielectric behavior of water in specific domains around the DNA, namely, the backbone and the minor and major grooves.In agreement with experimental [129,130] and simulation data [128] we observe that (phosphate backbone) ≈ 66 > (major groove)≈ 46 > (minor groove) ≈ 37.

Simulating Osmotic Isobaric Ensemble: Densely Packed DNA Arrays
To study more complex systems, composed of several DNA molecules one could in principle replicate the one-DNA molecular systems, presented in the previous section.However, in the densely packed DNA arrays, the DNA molecules are positioned so close to each other that there is not enough space to place the different resolution domains in between them.Therefore, to mimic the solvent grand canonical simulations at dense DNA packing, we have embedded the whole DNA array within a reservoir of CG water molecules and ions, see Figure 3.The DNA array and reservoir are separated by a semi-permeable membrane that allows the water molecules and ions (but not the DNA molecules) to freely pass between the two domains [73].

Computation of Osmotic Pressure
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) [131], that exerts a uniform compression on the DNA subphase [22].In MD simulations, the fixed osmotic stress of the pure water and salt solution on the DNA subphase can be implemented via a semipermeable membrane [132,133] surrounding a set of DNA molecules.The osmotic pressure of the DNA array is then determined by measuring the force exerted by the DNA subphase on the confining membrane.The force is computed via a repulsive part of a 10-4 Lennard-Jones interaction between the semipermeable membrane and DNA backbone of the following form [73] U wall (x) = 2πρ 2σ 12 5x 10 − where x is a perpendicular distance to the wall, σ and are the 12-6 Lennard-Jones parameters and ρ is the wall density.The interaction acts only between the DNA backbone atoms and the membrane to prevent a given DNA molecule from crossing it.The aqueous solution molecules (water, Na + , and Cl − ions) can nevertheless freely pass through.Thus, the system inside the semipermeable membrane is open, i.e., semi grand canonical, even though the total system is clearly closed.The osmotic pressure of the DNA array can then be computed as F wall /A, where F wall is the mean force on the semipermeable membrane and A is its area.For low DNA concentrations, we found that the osmotic pressure is low and does not depend on a specific type or valency of the counterions present.The observed independence of the osmotic pressure on the ion type at low DNA densities is connected with the pronounced ionic screening of electrostatic interactions at 1 M salt conditions beyond the DNA intersurface separation of ∼0.3 nm and is consistent with experiments [134].In the opposite regime of high DNA densities, the osmotic pressure difference between two different ion types is substantial for both lattices, while the presence of multivalent Spd 3+ results in pronounced lowering of the osmotic pressure.At small DNA separations, our quantitative results compare quite well with the experimentally observed values at 2 mM Spd 3+ [23] for the hex lattice.At larger separations, however, the pressure obtained with simulations while small, is usually overestimated and attributable to the small size of the simulated array as well as to the possible inaccuracies in the force field.We also observed that the initial lattice in the simulated array is not preserved at the lowest examined DNA concentrations.Hence, it is questionable whether they can be meaningfully compared to the experimental data in this density regime.

Hexagonal Columnar (H) → Orthorhombic (OR) Phase Transition
Interactions between dsDNAs have a strong orientational component [16,25] that can affect the collective packing symmetry of the ordered array [13,17].The ensuing azymuthal rotational correlations between DNA molecules appear to be of crucial importance in OR lattices and the loss of azimuthal ordering is found at approximately the same DNA concentrations as the H → OR phase transition.Therefore, the rotational analog of the Lindemann empirical criterion [73] is a good proxy to establish the point of H → OR phase transition.We estimated [73] this transition to occur roughly at DNA concentrations between 700-800 and 500-700 mg/mL for the Na + and Na + -Spd 3+ counterions, respectively, which also compares favorably with the experimental estimations, i.e., ≈670 mg/mL [10,11].A more consistent determination of the phase transition based on no additional assumptions or approximations could proceed from the free-energy landscape computed by enhanced sampling techniques.While this is an enticing direction it is hampered by much heavier demands on the computational time and resources.

Orientational Order Parameters
The pronounced orientational dependence of the hydrogen bonding interactions between water molecules and the consequent ordering of the molecular solvent orientations can be quantified by introducing various orientational order parameters.Three different types of these order parametrers, η (1,2,3) , have been investigated in related contexts.They are defined as where α is the relative angle between a chosen direction and the local orientation of the water molecule.The first two order parameters measure the average orientation of the dipole and the uniaxial quadrupole moments, being equal to η (1,2) = 0 for the random orientation in the bulk.The third order parameter quantifies the total quadrupole moment and has a nonzero value in the bulk.These order parameters can be used to examine the local ordering of interstitial water molecules residing between the DNA molecular surfaces [73].We have computed all three order parameters along the axis of a given DNA pair for water molecules located in a cuboid (defined by DNA-DNA interaxial spacing, DNA diameter, and its long axis) between a given DNA pair.We found that upon the increase of the DNA concentration the local ordering of water molecules, as quantified by the order parameters η (1,2,3) , shows distinct features as well as a noticeable symmetry.In particular, the η (2,3) are symmetric w.r.t. the midpoint and indicate progressive orientational layering, while η (1) is antisymmetric and quantifies the ordering strength of the opposed DNA surfaces.η (2) ≥ 0 for vicinal water indicates that it is oriented preferentially perpendicular to the DNA surface, while η (2) ≤ 0 for the interstitial water indicates it exhibits an almost uniform orientation.These conclusions tally favorably also with water ordering and resulting hydration forces between lipid membranes and polar surfaces in general [81,135].We consequently conclude that the hydration forces contribute essentially also to the osmotic pressure at large DNA concentrations.This conclusion is further supported by comparing the total osmotic pressure in the system consisting of DNA and the aqueous solution, with the one obtained for DNA only, without any contribution of the solvent.Not only is the latter smaller, but can even show an opposite (negative) sign, which indicates that the aqueous solvent adds an essential contribution to the balance of forces in DNA arrays based on the water-mediated hydration interaction, as has been demonstrated in previous experiments on DNA arrays [4,10].

Conclusions and Perspectives
In this review, we presented an overview of computational approaches to simulate densely packed DNA arrays.We highlighted the structural characterization of DNA phases using orientational order parameters, the accurate determination of dielectric properties of solvent and DNA molecules, and modeling of the solvent in the grand-canonical ensemble fashion, as key ingredients to adequately study phenomena such a phase transitions between different high-density DNA mesophases, e.g., with hexagonal and orthorhombic local packings.We presented our attempts to tackle these issues using the multiscale MD method AdResS.
To perform a truly open grand-canonical molecular simulations one should take recourse to the open boundary molecular dynamics (OBMD) [57,116,[136][137][138] that permits simulations directly in the grand-canonical ensemble.The OBMD method combines features of AdResS and open MD [139,140], which enable the exchange of energy, momentum, and matter with the external environment through the imposed arbitrary time-dependent external pressure tensor.As discussed in this review, water can be considered as a dielectric continuum from the second DNA coordination shell onward, while the fluctuating cloud of ions surrounding the DNA extends to a wider region.Thus, to properly implement the OBMD simulations that include also the ionic component of the bathing solution, one would need to consider larger all-atom regions that would eventually undermine the original idea of an efficient simulation.To circumvent this limitation one could couple the explicit all-atom salt solution, used in the cylindrical layer surrounding a DNA molecule, and the implicit water with explicit ions, used in the outer part of the simulation box.This would allow for efficient molecular simulations of biomolecules solvated in bathing salt solutions at any ionic strength conditions.The work along these lines is currently underway.

Figure 1 .
Figure 1.An atomistic DNA molecule embedded in a dual-resolution salt solution.The proximal solvent around the DNA molecule is modeled using the atomistic level of detail whereas the distal water molecules are represented as grey beads, i.e., one bead corresponds to one water molecule.The salt ions are shown in green (Na + ) and blue (Cl − ).Reprinted (adapted) with permission from Ref. [70].Copyright (2015) American Chemical Society.

Figure 2 .
Figure 2. Schematic depiction of a simulation box with an atomistic DNA in the center surrounded by MARTINI salt solution.Two levels of resolution are used for solvent molecules.AT level of resolution is used for proximal solvent molecules, within a certain radius from the DNA's center of mass.The distal water molecules are represented as blue MARTINI CG beads, where one CG bead corresponds to four water molecules.Na + and Cl − ions are shown in orange and yellow, respectively.Adapted from Ref. [71].

Figure 3 .
Figure 3. Top-down view on the simulation box containing 16 DNA molecules arranged on a hexagonal lattice.The DNA molecules are at all times modeled with the high atomistic resolution.For the solvent, we use two levels of resolution, i.e., the atomistic resolution in the rhombic region containing the DNA molecules and the coarse-grained representation outside.Adapted from Ref. [73].