Ion distributions around left- and right-handed DNA and RNA duplexes: a comparative study

The ion atmosphere around nucleic acids is an integral part of their solvated structure. However, detailed aspects of the ionic distribution are difficult to probe experimentally, and comparative studies for different structures of the same sequence are almost non-existent. Here, we have used large-scale molecular dynamics simulations to perform a comparative study of the ion distribution around (5′-CGCGCGCGCGCG-3′)2 dodecamers in solution in B-DNA, A-RNA, Z-DNA and Z-RNA forms. The CG sequence is very sensitive to ionic strength and it allows the comparison with the rare but important left-handed forms. The ions investigated include Na+, K+ and Mg2 +, with various concentrations of their chloride salts. Our results quantitatively describe the characteristics of the ionic distributions for different structures at varying ionic strengths, tracing these differences to nucleic acid structure and ion type. Several binding pockets with rather long ion residence times are described, both for the monovalent ions and for the hexahydrated Mg[(H2O)6]2+ ion. The conformations of these binding pockets include direct binding through desolvated ion bridges in the GpC steps in B-DNA and A-RNA; direct binding to backbone oxygens; binding of Mg[(H2O)6]2+ to distant phosphates, resulting in acute bending of A-RNA; tight ‘ion traps’ in Z-RNA between C-O2 and the C-O2′ atoms in GpC steps; and others.

The long-range Coulomb interaction was evaluated by means of the Particle-Mesh Ewald (PME) method [8] with a 9Å cutoff and an Ewald coefficient of 0.30768. Similarly, the van der Waals interaction were calculated by means of a 9Å atom-based nonbonded list, with a continuous correction applied to the long-range part of the interaction. The production runs were generated using the leap-frog algorithm with a 1 fs timestep with Langevin dynamics with a collision frequency of 1 ps −1 . Data was saved every picosecond of the simulation. The SHAKE algorithm was applied to all bonds involving hydrogen atoms. The nucleic acid structures generated by means of these constant temperature (T = 300K) and pressure (1 atm) runs were analyzed by means of the 3DNA (v2.1) package [9] and the PTRAJ package in AMBERTOOLS [1].
The initial coordinates of the nucleic acid structures were generated as follows. The usual right-handed B-DNA and A-RNA structures were generated using the NAB package in AMBERTOOLS (ver.12) in Arnott's conformation [10,11], while the left-handed Z-DNA structure was available from our previous investigations [12]. For the Z-RNA structure, we first generated the (CG) 3 structure of Z-DNA, and then modified it to generate the corresponding RNA structure. The resulting structure was then equilibrated for 10 ns in an explicit solvent environment with a neutralizing number of Na + ions. The average structure was then analyzed using 3DNA and compared to the NMR-based Z-RNA structure found in the PDB data bank (PDB id: 1T4X) [13]. Since the local base-pair and step parameters of this equilibrated structure were found to be quite close to the experimental structure, we simply replicated the (CG) 3 structure to generate the initial (CG) 6 Z-RNA duplex. The duplexes were first equilibrated in a smaller water box (with neutralizing ions) and then placed in a larger cubic box of ∼ 82Å side, filled with a suitable number of water molecules. Such a large box is necessary, since cylindral distribution functions need to be calculated to a distance of at least 30Å. To achieve electroneutrality, 22 water molecules were replaced by Na + cations. For simulations involving high salt concentrations, the number of water molecules replaced by salt anions and cations was dictated by the system size and density. In order to minimize any correlation effects on the initial distribution of the ions, these were placed at random at a distance of 10Å away from each other and the duplex structure. molarity is used to indicate excess salt (over the neutralizing ions). Thus, we will refer to the case (i) above as "zero salt", although it has 0.06M Na + , due to the neutralizing ions.
The energy of each final system of oligonucleotide, waters and ions was then minimized with the nucleic acid and ions fixed, followed by further minimization where atomic motion was allowed. Subsequently, the temperature was gradually raised at constant volume from zero to 300 K over a period of 50 ps, with the DNA/RNA structure and ions restrained with a harmonic constant of 50 kcal/mol. The resulting structure was then further equilibrated at 300 K for 100 ps. The restraints on the duplex atoms and ions were then gradually reduced by decreasing the restraining harmonic constant in 5 steps during equal intervals of time. The final configuration was then used for a 1 ns constant pressure simulation which was completely unrestrained, thereby adjusting the box size so as to achieve a system density of approximately 1 g/ml. This was then followed by 120 ns of constant pressure and temperature production runs. Equilibrium data was collected from only the last 100 ns of these runs, since a minimum of 20 ns is required in order to stabilize the ion distribution.
Our analysis was focused on calculating the distribution of the mobile ions around the nucleic acid structures. To that end, we calculated the diffusion coefficients for the ions, the cylindrical and radial distribution functions, and carried out an analysis of the efficacy of different sites to localize the ions.
The diffusion coefficient (D) for the ions was obtained by conventional means: i.e., D is obtained by calculating the slope of the mean-square displacement as a function of time. We calculated this quantity over a 10 ns time interval, taking data every 10 ps, and averaging the result over all the ions. Turning to the radial distribution function (RDF), this is defined in terms of the distance r of an atom B from a central atom A: where ρ B is the bulk density of B and ρ B (r | r A = 0) is the conditional distance-dependent density of atoms of type B a distance r from central site A. In practice the RDF is calculated as: where N B (r, Δ) is the average number of atoms of type B located between distances r − 1 2 Δ and r + 1 2 Δ from central atom A, and V B (r, Δ) is the volume of the spherical slice between r − 1 2 Δ and r + 1 2 Δ. RDFs are particularly useful for investigating the binding properties between specific atoms such as negatively charged nitrogens and oxygens on the duplexes and positively charged ions. In these RDFs the location of the peaks reveals the typical binding distances. In the results presented here, the RDFs between two atom types were computed as an average over all the equivalent atom pairs along the duplex. An ion can be attracted to more than one atom in the duplex, and thus it can contribute to the RDFs of different atom types in the duplex. A particularly useful variant of RDFs involves the cylindrical distribution functions, which track the ion concentration measured radially outwards from the central axis of the duplex. The calculation of these functions is the same as for the RDFs, except that r is measured on the plane perpendicular to the duplex axis and V B (r, Δ) now represents the volume of a cylindrical slice between r−Δ/2 and r+Δ/2 from the chosen central axis. In our analysis, the global z-axis was calculated by the package 3DNA, as in SCHNAaP [14]. The vectors used to define the axis are a combination of C1' and G-N9/C-N4 atoms along the same strand, as developed by Rosenberg et al. [15]. Calculating the cylindrical distribution function not only gives insight into the properties of the ion distribution, but also provides for a good check for the convergence of the simulation.

Results
Diffusion Constants. Table 2 summarizes the calculated diffusion coefficients for cations in the presence of the different duplexes. Generally speaking, these results are in good agreement with values obtained in B-DNA studies [16]. First, we consider the results for the Na + ion. At low ion concentration, the diffusion coefficients for the right-handed structures are higher than for the corresponding left-handed structures. For comparison, the calculated diffusion coefficient for 22 free Na + ions in TIP3P waters is 1.56 × 10 −9 m 2 /s, which is higher than the experimental value of 1.33 × 10 −9 m 2 /s for ions in pure water [17].    [18], and theories based on various versions of the Poisson-Boltzmann equation. At short distances near the duplexes, these theories are expected to fail because the atomic details matter [19].
However at larger distances, continuum theories can work. To probe this point, we considered a simple model consisting of two concentric cylinders (inner cylinder radius r=a; outer cylinder, r=b), with a constant charge density embedded on the inner cylinder and the electric field vanishing on the outer cylinder. The potential generated by such a model may be solved analytically using the linearized Poisson-Boltzmann equation [20]. The application of this model to our systems gives excellent results at larger distances. For example, SI Fig.11 shows results for 0.4 M Na + (radii a=12.0Å and b=29.5Å ) and shows that this simple continuum model well represents the cylindrical ion distribution for distances larger than about 12.0Å (RNA) and 14.0Å (DNA).    [20]. As may be expected, good fits are obtained away from the central axis.