Effect of initial ion positions on the interactions of monovalent and divalent ions with a DNA duplex as revealed with atomistic molecular dynamics simulations

Monovalent (Na+) and divalent (Mg2+) ion distributions around the Dickerson-Drew dodecamer were studied by atomistic molecular dynamics (MD) simulations with AMBER molecular modeling software. Different initial placements of ions were tried and the resulting effects on the ion distributions around DNA were investigated. For monovalent ions, results were found to be nearly independent of initial cation coordinates. However, Mg2+ ions demonstrated a strong initial coordinate dependent behavior. While some divalent ions initially placed near the DNA formed essentially permanent direct coordination complexes with electronegative DNA atoms, Mg2+ ions initially placed further away from the duplex formed a full, nonexchanging, octahedral first solvation shell. These fully solvated cations were still capable of binding with DNA with events lasting up to 20 ns, and in comparison were bound much longer than Na+ ions. Force field parameters were also investigated with modest and little differences arising from ion (ions94 and ions08) and nucleic acid description (ff99, ff99bsc0, and ff10), respectively. Based on known Mg2+ ion solvation structure, we conclude that in most cases Mg2+ ions retain their first solvation shell, making only solvent-mediated contacts with DNA duplex. The proper way to simulate Mg2+ ions around DNA duplex, therefore, should begin with ions placed in the bulk water.


Introduction
Ion interactions with nucleic acids (including both DNA and RNA) are an important subject that has received extensive investigations since the early 1960s (Elson, 1965). Positively charged cations may interact with highly negatively charged nucleic acids via simple electrostatic interactions to screen the electrostatic repulsion along the nucleic acids and assist their folding and/or compaction. Cations may also bind at specific sites and become integral parts of the structures, possibly playing important enzymatic roles (Frederiksen & Piccirilli, 2009;Freisinger & Sigel, 2007;Muller, 2010). For example, function of the hammerhead ribozyme is typically lost without specific divalent cations, but may be restored at molar concentrations of monovalent salts (Murray, Seyhan, Walter, Burke, & Scott, 1998;O'Rear et al., 2001). Recently, variants of the 8-17 DNAzymes have been designed to have enzyme activity dependent on the binding of a specific divalent cation and such properties have been exploited for use as sensitive metal ion sensors (Xiang, Tong, & Lu, 2009;Xiang, Wang, Xing, Wong, & Lu, 2010). Knowing precisely where and how cations interact with nucleic acids will help elucidate the struc-ture of nucleic acids which is increasingly recognized as playing an important role in many cellular functions.
Fully atomistic molecular dynamics (MD) simulations can provide a wealth of microscopic information about the ion interactions with DNA and RNA. Extensive allatom MD simulations of monovalent ions interacting with DNA have been reported that have provided important insight (Ponomarev et al., 2004;Rueda et al., 2004;Savelyev & Papoian, 2006;Young et al., 1997). In comparison, atomistic simulations of divalent ions like Mg 2+ interacting with the nucleic acid are fewer. As more and more structures of nucleic acids with bound divalent Mg 2 + ions become available, ) MD simulations of nucleic acids with divalent ions are routinely performed (Auffinger, Bielecki, & Westhof, 2004;Kirmizialtin & Elber, 2010;Kirmizialtin, Pabit, Meisburger, Pollack, & Elber, 2012;Kirmizialtin, Silalahi, Elber, & Fenley, 2012;Lee et al., 2008). While it is recognized that representing Mg 2+ as a point charge without accounting for polarization or charge-transfer effects leads to an inadequate description of cation behavior (Sponer et al., 2000), this type of MD simulations with explicit divalent ions will still be very useful for elucidating the structure and functions of nucleic acids. Hence, it is desirable to have a better understanding on how to properly simulate the divalent ions interacting with nucleic acids given the known limitations of force fields.
In the current study, we performed atomistic MD simulations of the Dickerson-Drew dodecamer of DNA immersed in monovalent (Na + ) and divalent (Mg 2+ ) salts, respectively. We examine the limit of current state-of-art atomistic MD techniques to model the interaction of cations with nucleic acids, especially the magnesium ions. We have tested the effect of initial cation positions on the simulated ion atmosphere surrounding the nucleic acids. Special attention was paid to analyze the solvation shell structures around the metal ions and how that may be disrupted upon binding with the nucleic acids.

Methods
Atomistic MD simulations of the DNA in explicit water and ions were carried out mostly using the parm99 force field and partly with parmbsc0 and parm10 with the AMBER simulation package (Case et al., 2008;Perez et al., 2007;Wang, Cieplak, & Kollman, 2000;Zgarbova et al., 2011). A B-form DNA structure of sequence d(CGCGAATTCGCG) 2 was created with the nucgen module of AMBER. The DNA duplex was solvated with approximately 25,000 TIP3P water molecules, forming a $100 Â 100 Â 100 Å 3 cubic box (Jorgensen, Chandrasekhar, Madura, Impey, & Klein 1983). Then, a specific number of cations (Na + or Mg 2+ ) were added to each system using two different methods of establishing their initial conditions. In type A simulations, we used the default parameters of the Leap module in AMBER, which places cations in the simulation box based on the calculated electrostatic grid, resulting in many cations being placed very near the DNA duplex. In type B simulations, the cations, after an initial placement by Leap, were randomly moved so that they were no closer than 35 Å from the DNA and 3 Å from each other. For both simulation types, the number of anions (Cl À ) required to achieve system neutrality was placed in the box with Leap after placement of the cations; for type B simulations, the placement of the anions was performed before cation randomization. As we will show later, the dynamics and distributions of the cations were found to be dependent on their initial placements.
We investigated two sets of ion parameters for Na + and Mg 2+ ions. The older versions of AMBER used the ions94 parameter set adapted from the works of Åqvist and Smith and Dang (Åqvist, 1990;Smith & Dang, 1994). This version of ion parameters has been used in a number of studies (Feig & Pettitt, 1999;Ponomarev et al., 2004;Rueda et al., 2004;Young et al., 1997), but it became known that the use of ions94 parameters led to the formation of salt crystals in long time scale molecular dynamic simulations (10 ns) at a concentration well below the known experimental saturation concentration (Auffinger, Cheatham, & Vaiana, 2007;Chen & Pappu, 2007;Savelyev & Papoian, 2006;Savelyev & Papoian, 2007). This problem led to the development of new ion parameters for monovalent ions, by Joung and Cheatham (ions08), that are currently the default choice in AMBER 12 (Joung & Cheatham, 2008). The parameters for Mg 2+ were not updated in AMBER, but Allnér and Villa recently developed Mg 2+ parameters that better reproduce experimental data on water exchange kinetics (Allnér, Nilsson, & Villa, 2012). They started with parameters for Mg 2+ available in CHARMM27 and adjusted the repulsive LJ term to reproduce experimental data on water exchange kinetics. The new parameters they developed for Mg 2+ ions have a weaker attractive interaction and a larger radius compared with the old parameters. Table 1 summarizes the Lennard-Jones parameters for different models of ions investigated in the current study, where ɛ and r ii are the depth and the position of the minimum of the LJ potential. Most of the simulations were performed with ions94 parameter, unless otherwise stated.
All systems were first minimized to remove bad initial contacts with harmonic restraints on the DNA for 1000 steps, followed by 5000 steps of unrestrained minimization. Afterward a 10 kcal/(mol⁄Å 2 ) positional restraint was used throughout the equilibration and production runs to hold the DNA fixed. With this restraint force, the heavy atoms on DNA will move less than 0.2 Å. The equilibration stage includes warming of the system temperature from 0 K to 300 K over a span of 20 ps with constant volume, followed by explicit NPT ensemble simulation at 300 K for 10 ns. The production run was also performed in NPT ensemble simulation with pressure of 1 bar using Langevin thermostat with a collision frequency of 1.0 ps À1 and pressure relaxation time 2.0 ps. A time step of 2 fs was used throughout with SHAKE constraints on all bonds involving hydrogen atoms (Miyamoto & Kollman, 1992;Ryckaert, Ciccotti, & Berendsen, 1977). Long-range electrostatics were treated with the particle mesh Ewald summation with a 10 Å direct space cut-off (Darden, York, & Pedersen, 1993). The NPT simulations were performed for 50 ns for nearly all systems and the production runs were taken from the last 40 ns. This length of production run is judged to be sufficient based on convergence tests using either number of bound ions or ion density distribution around DNA duplex, when except there is a strong dependence on the initial placements of divalent ions. We will discuss this effect shortly. Trajectory analysis was carried out using the Ptraj program distributed with the Amber-Tools package as well as some scripts and Fortran codes developed by us. System trajectories were visualized with Visual MD (VMD) 1.9.1 (Humphrey, Dalke, & Schulten, 1996). The DNA's RMSD during the production runs from the initial minimized structure was less than 0.3 Å during the entire simulations. In the current study, DNA can be considered as essentially fixed since we are not interested in the dynamics of the DNA, but only interested in how ions may be bound with DNA. We note that several recent MD simulations of ion interactions with nucleic acids also used fixed DNA/RNA structures (Kirmizialtin & Elber, 2010;Kirmizialtin et al., 2012aKirmizialtin et al., , 2012b. Table 2 summarizes the systems investigated in the current study. The DNA structures used in simulations are nearly the same. Simulations performed with parm99 or parmbsc0 led to very little difference in the DNA structure since the DNA is restrained. The main focus of the current paper is on the initial starting coordinates of ions.

Effect of initial positions of cations on the number of bound ions in MD simulations
We are interested in the comparison of number of bound ions around a DNA duplex revealed in different MD simulations. A B-form DNA duplex is a highly negatively charged polyelectrolyte, and the counterions near the DNA feel a strong attraction because of the electrostatic interactions. This concept of bound ions was proposed by Manning and is very useful to understand the properties of polyelectrolytes (Manning, 1979(Manning, , 1980(Manning, , 2007. If one can determine the electrostatic potential, ϕ (r), according to the Poisson-Boltzmann equation, then any cations that reside in a location, r, where the product of its charge, q, and potential, ϕ(r), preserves this relationship |qϕ|>k B T (k B is the Boltzmann constant and T is the temperature) may be considered bound since the electrostatic attraction between the cation and the DNA duplex is larger than the thermal energy of the ions. However, there is no easy way to monitor the number of bound ions in MD simulations using this definition. Hence, we used a slightly different approach: we monitored the number of cations within 10 Å from the closest C1′ atom on the DNA duplex. This quantity can be calculated during the MD simulations at any time. The C1′ atom is on the sugar ring where the base connects, and as each base possesses only one, the C1′ is quite convenient for this analysis. The C1′ atom is approximately 5 Å away from the phosphate groups which define the outer surface of the DNA duplex. To further verify that using this cut-off distance would yield the bound ions, we made use of a software suite, APBS, to numerically solve the Poisson-Boltzmann equation around a DNA duplex. The details of the APBS calculation will be presented in a future report, but here we want to show that ions within this cut-off lie within the potential surface that defines |qϕ| > k B T. Figure 1(a) renders multiple frames of Na + ions within the cut-off distance in the MD simulations with the calculated electrostatic potential surface of 1 k B T, whereas Figure 1(b) renders the simulations with Mg 2+ ions within the cut-off distance with the electrostatic potential surface of 0.5 k B T. Only in the case of Mg 2+ ions, a few ions within the cut-off distance were found to lie outside the electric potential surface, but such events are rare. We used this number of bound ions as a way to check whether the simulations starting with different initial cation positions would yield equilibrium values. Figure 2(a) and (b) present the number of cations within the cut-off distance determined from MD trajectories for Na + and Mg 2+ ions with two different starting positions. For the Na + simulations, after the initial equilibration stages, the number of ions within 10 Å from any of C1′ atoms reached a plateau value and was independent of the initial positions of the ions. In fact, in all Na + ion simulations the ion atmosphere converged within the 40 ns production window reflected by only slight deviation in the number of cations and anions within the cutoff distance as simulation time progressed. However, the number of Mg 2+ ions within the cut-off distance depended on the initial placement of the cation; simulations with ions placed by Leap (type A) had more ions bound to the DNA than those starting with Mg 2+ ions initially placed away from DNA structures (type B). Simulations of types A and B never reached the same average numbers of ions for divalent Mg 2+ ions even in simulations as long as 50 ns. Although MD simulations with two different initial starting positions did not converge to the same value, the number of bound ions determined for either type A or type B simulations converged within the production time. This is reflected in quick initial drop in type A or rising in type B in the number of ions, which indicates a relaxation time of no more than 5 ns. As will be discussed below, the difference in the number of bound Mg 2+ ions in these two types of simulations is due to the solvation structure of Mg 2+ ions. In type A simulations, Mg 2+ ions were initially placed too close to DNA duplex and formed permanently bound ions which then could not be equilibrated.

Anomalous behavior for Mg 2+ ions is related to their solvation shell structure
The inability to reach equilibrium between two simulation types for Mg 2+ ions led us to investigate the origin of this behavior. We examined the interactions of Mg 2+ ions with the nucleic acids. The analysis revealed that there was a fraction of Mg 2+ ions in type A simulations that were permanently bound to nucleic acids over the entire 50 ns of trajectories ( Table 2). The distances observed for permanently bound Mg 2+ ions were equal or less than 2.0 Å, which is the distance from the Mg 2+ ion to an oxygen atom in its first solvation shell. Essentially, these Mg 2+ ions are in direct contact with these atoms on nucleic acids and have become partially desolvated. Most of these sites are on the phosphate groups of oxygen atoms (O1P/O2P) with a few cases in the major groove on the O6 of guanine or the minor groove of O2 atoms on cytosine and thymine.
It is known experimentally that Mg 2+ forms an octahedral first solvation shell with residence times in the microsecond range (Bleuzen, Pittet, Helm, & Merbach, 1997;Neely & Connick, 1970). By using a smaller simulation system with just Mg 2+ ions immersed in TIP3P water, we tested and confirmed that Mg 2+ ions immersed in TIP3P water bath formed an octahedral first solvation shell with six bound waters in the AMBER parm99 force field with ions94 parameters. These solvation shell water molecules never exchanged with other water molecules in solution over 300 ns of MD simulation at 300 K. We calculated radial distribution functions (RDF) of the oxygen atoms of the bulk water around Mg 2+ ions. The radial distribution function shows a first peak at a distance of 2.0 Å and the second peak from 3.5 to 5 Å. The first peak is related to the first solvation shell and integration of the first peak intensity yields exactly six water molecules permanently bound to the Mg 2+ ions. The second peak is the second solvation shell and water molecules in the second shell were more loosely bound and number of water molecules reaches up to 30.
We then examined the number of water molecules within the first solvation shell of the Mg 2+ ions in type A and type B MD simulations. Figure 3 presents the number of cations having different numbers of water molecules in the first solvation shell for both type A and type B simulations using a cut-off distance of 3.4 Å for the first solvation shell. In type B simulations, all Mg 2+ ions were found to have exactly six water molecules within the 3.4 Å cut-off distance. However, in type A simulations, most of permanently bound Mg 2+ ions had only four water molecules in the first solvation shell. A few of them even had zero water molecules and these ions were found to be entirely buried in groove region of DNA. For the rest of nonpermanently bound Mg 2+ ions in type A simulations (i.e., the mobile Mg 2+ ions), most of them had only five water molecules instead of six. A closer examination revealed that these partially desolvated, mobile Mg 2+ ions were paired with one or two permanently bound Cl À ions with a coordination distance of 2.55 Å. The difference in the first solvation shell for the mobile Mg 2+ ions in two types of simulations can be traced to how the ions were initially introduced in the simulations. Recall that in type A simulations, both Mg 2+ and Cl À ions were placed by Leap according to an electrostatic grid. These ions were placed fairly close to the nucleic acids and were close to each other. Once the MD simulation began, the nearby Mg 2+ and Cl À ions formed an ion pair that never dissociated during the MD simulations. As a result, the number of water molecules within the first solvation shell became less than six even when the ion pair drifted away from the nucleic acid. In type B simulations, the initial locations of Mg 2+ ion were swapped and placed at least 35 Å away from the nucleic acid. The Cl À ions were left at their original positions placed by Leap, near the nucleic acid. Hence, in type B simulations, Mg 2+ ions were far away from the nucleic acids and the Cl À ions and these Mg 2+ ions were surrounded by water molecules. After the simulation began, surrounding water molecules quickly formed octahedral first solvation shells with the Mg 2+ ions and this first shell never exchanged with other solvent molecules. These analyses suggest that one should be extremely careful when deciding on how to initially place Mg 2+ ions and Cl À ions in MD simulations. More discussion on how to properly simulate Mg 2 + ions in explicit solvent water system will be presented later on.

Monovalent Na + ions can be equilibrated regardless their initial positions
Compared with Mg 2+ , data in Figure 2 suggest that simulations with monovalent Na + can reach equilibrium regardless of the ions' initial positions. This again can be understood from the solvation structure around the Na + ion. The first solvation shell of the Na + ion is more loosely bound. The first shell peak distance of the radial distribution function of oxygen atoms to the Na + ion in the bulk water was at 2.4 Å, further away than first solvation shell around Mg 2+ . We found that within 3.0 Å cut-off, the mean number of water molecules surrounding the Na + ion was 5.83 agreeing well with previously reported values (Koneshan, Rasaiah, Lynden-Bell, & Lee, 1998). The average residence time of water molecules in the first solvation shell was reported to be $10 ps to 20 ps (Impey, Madden, & Mcdonald, 1983;Koneshan et al., 1998). Within our simulation protocols, we found similar residence times for the water in the first solvation shell. During the MD simulations, the water molecules surrounding the Na + exchanged frequently with the bulk water molecules. Additionally, close contacts between Na + and Cl À ion pairs were observed with an average distance around 2.78 Å that lasted about 250 ps. So, unlike Mg 2+ and Cl À ion pairs, Na + and Cl À ion pairs were able to dissociate and reform in nanosecond-long MD simulations. Therefore, the ion distribution of monovalent salt around the DNA was independent of the initial starting positions of the ions, but the ion distribution of divalent Mg 2+ obtained in MD simulations was sensitive to their initial starting configurations.

Mode of interactions of Na + ions with DNA
Although how Na + ion interacts with DNA has been investigated extensively by others, we performed these analysis for two purposes: 1) to compare against the behavior for Mg 2+ ions; and 2) to validate our simulation and analysis methods against the reported literature. We first computed the RDF of the closest ions around any atom of the DNA. The RDF was calculated by treating DNA atoms as solute and the cations as solvent molecules and only the closest distances are counted in the calculation. The obtained RDF was not normalized by the ion density since we are only interested in the peak of the RDF which reflects the closest contact distances. Figure 4 presents the RDF obtained for Na_31 system for both simulation types, which are nearly indistinguishable. The closest contact between Na + ions with any atoms on DNA was $2.4 Å, in agreement with other studies (Cheatham & Young, 2000). This was a direct contact between Na + ions and DNA atoms (not water mediated). Guided by the RDF's first peak, a 3 Å cut-off was used to generate Na-DNA contact data using the hbond command in Ptraj. Table 4 summarizes prominent contact sites on the DNA for each ion. Prominent contacts are defined as the ion-DNA contact with the highest occupancy during the MD simulation for each cation. Most contacts are with the phosphate O1P/O2P atoms and some with N7 on the purines (guanine or adenine) in the major groove. These contacts are also long lived, lasting hundreds of picoseconds. Additional contacts at minor grove of O2 on cytosine and thymine and a small fraction of contacts on the O4' in the sugar ring were also observed.
We further analyzed sequence dependence of binding events on DNA duplex where binding is defined as any contacts made by a Na + ion with any atoms except hydrogen on DNA within 3.0 Å. This was done by calculating the percentage of occupied frames on each individ-  ual residue based on all the contacts generated with hbond command with a 3.0 Å cut-off. Supplementary Figure S1 presents the results of the analysis. We have considered two types of contacts: contacts made with each residue with any DNA atoms and contacts made with only oxygen atoms on the phosphate group (O1P/ O2P). We find that contacts made with O1P/O2P constitute a good fraction of all contacts. The simulations of type A and type B showed a similar pattern of sequence-dependent enhancement for guanine and adenine, although the extent of enhancement varied for the two simulation types. The preference for these sites (guanine and adenine residues) in the grooves correlates well with the N7-cation binding behavior at the AT-tract discussed previously and has been observed in other MD simulations and X-ray crystal structures (Egli, 2002;Sethaphong, Singh, Marlowe, & Yingling, 2010;Singh, Sethaphong, & Yingling, 2011). Conversely, the contacts made with phosphate groups have much reduced sequence dependence. Na + ions interact with the DNA phosphate groups fairly uniformly. Figure 5 presents the atomic radial distribution function for the Mg 2+ ions around any DNA atom. In type A simulations, the RDF was dominated by the first peak at 2 Å which was due to the aforementioned permanently bound Mg 2+ ions. In type B simulations, the closest contact was at 4 Å, due to solvent-mediated interactions between Mg 2+ ions with DNA. For high concentration type A simulations, a solvent-mediated peak can also be observed upon focusing on that region indicating solvent-mediated contacts also exist (see insert Figure 5(a)).

Mode of interactions of Mg 2+ ions with DNA
Binding sites for type A simulations were discussed previously and were summarized in Table 3. A cut-off of 5 Å was used to generate contact data in accordance with the first peak in RDF of type B simulations. The major interaction sites in type B simulations are summarized in Table 5. The most abundant type is at phosphate O1P/ O2P sites, which on average lasts in the range of hundreds of picoseconds. The next abundant type is at the O6 or N7 atoms on purines in the major groove and at O2 on pyrimidines in the minor groove. The longest binding event observed reached as high as $20 ns at the minor groove. The sequence dependence of the binding event for Mg 2+ ions in type B simulations was analyzed following the same protocol used for Na + ions. Supplementary Figure S2 presents the data. Sequence enhanced binding was observed for residues 4, 8, 16, and 20 for the 91Mg_B system. The preferred residues were identical due to the palindromic sequence; residues 4 and 16 were guanines at the border of the AT-tract and 8 and 20 were thymines. The binding at the phosphate groups, however, was largely uniform, with no significant sequence dependence.
Binding events for the divalent Mg 2+ ions were on average much longer than for Na + ions, with the longest exceeding 20 ns. These Mg 2+ ions all having fully solvated first shells were unable to form direct coordination complexes; they only formed solvent-mediated contacts with both minor and major groove atoms as well as the phosphate groups. In order for these stable water-mediated contacts to form, the second solvation shell must be perturbed. The nanosecond-scale contacts required the ejection of $5 second shell waters. Figure 6(a) and (c) illustrate this correlation between stable binding events and ejection of second shell waters. The cation can closely approach the duplex with a full second shell, but cannot remain there for long without ejection of second shell waters. Figure 6(b) and (d) present representative binding motifs at the major and minor grooves, where the fully solvated Mg 2+ ion made three contacts with the nucleic acid via first solvation shell waters. Binding  events with the O1P/O2P were significantly short lived (average duration time $200 ps) when compared to the some of the binding events at grooves. This would seem counterintuitive based on a solely electrostatic viewpoint due to the greater charge on O1P/O2P. A closer examination reveals that a fully solvated Mg 2+ ion can only make two water-mediated contacts with nucleic acid atoms at the phosphate groups, and consequently the binding remains short lived when compared against the binding at the grooves.

Sensitivity of simulations to force fields and ion parameters
Lastly we investigate sensitivity of simulation results on the force fields used for the nucleic acids and the ions. Extensive discussions on the influence of monovalent ion parameters and force fields for nucleic acids used in MD simulations on the resulting DNA-ion interactions have been reported in the literature (Chen & Pappu, 2007;Noy, Soteras, Javier Luque, & Orozco, 2009;Pérez, Lankas, Luque, & Orozco, 2008;Rueda et al., 2004;Varnai & Zak-rzewska, 2004). We focus here on the number of bound ions and roles of solvation shell structure during the binding. The use of parm99 or parmbsc0 did not have much effect on these properties since in our simulations the DNA duplex is restrained. The difference in parm99 and parmbsc0 in AMBER simulation affects the alpha/ gamma transition in the phosphate backbone which has been linked with conformational stability of B-DNA vs. A-DNA. (Perez et al., 2007) In our simulations, the initial structure of DNA duplex, created using nucgen model, adopts the B-form. The RMSD of DNA duplex during the entire MD simulation was less than 0.3 Å, regardless of whether we have used parm99 or par-mbsc0 force fields. The ion parameters had some impact on the observed interactions between ions and nucleic acids. For the monovalent Na + ions, the LJ potential in ions08 for Na + has a deeper well at a smaller radius (see Table 1) than in ions94. Consequently, we observed that the first peak of RDF between the closest Na + ions with any DNA atom shifted to a slightly smaller radial distance with a higher intensity (see Figure 7). Analysis of contact between Na + and DNA atoms reveal that the average time in which a  Na + ion binds with DNA atom is longer in simulations with ions08 than with ions94. However, the number of bound ions around the DNA duplex, as defined in section (A), had little dependence on the ions' parameters. The two simulations yield the same number of bound ions as shown in Figure 8. The number of bound ions is determined primarily by the long-range electrostatic interactions between ions with DNA, and hence it is not sensitive to the parameters in the short-range LJ interactions. For the Mg 2+ ions, similar problems of equilibration with respect to different initial starting positions were observed with the use of Allnér and Villa's parameter. This is not surprising since Mg 2+ ions form a tightly bound first water shell in both parameter sets. The exchange of water molecules in the first solvation shell is too slow to occur in these standard MD simulations. Initial improper positions of the Mg 2+ ions will lead to the same artifacts in the simulations. However, when we allow the Mg 2+ ions to first form water solvation shell in the bulk and then proceed with the MD simulations, the number of bound Mg 2+ again is found to be independent of the ion parameters. Figure 6. Panel (a) depicts the binding events of one Mg 2+ ion in MD simulations of 31Mg_B, monitored by the minimum distance between Mg 2+ ion and C1′ atom on the DNA duplex (blue circles) and total number of waters in the second solvation shell (red squares). Whenever distance between Mg 2+ and C1′ atom drops below 10 Å, the ion can be considered as bound. Notice how stable water-mediated binding does not occur until there is significant perturbation of the second shell. Panel (b) illustrates the binding of a fully solvated Mg 2+ ion making 3 water-mediated contacts with three different nucleic acid residues (spheres). Noninteracting first shell waters and nucleic acid atoms are represented with ball and stick. Image corresponds to binding events during the last 5000 ps in (a) and a prominent contact with the major groove. Panels (c) and (d) are similar to Panels (a) and (b), except it is for 91Mg_B. Image in panel (d) corresponds to binding events during 10,000-30,000 ps in panel (c) with a prominent contact at the minor groove.

Conclusion
We have investigated how to properly simulate the ion atmosphere formed by Na + and Mg 2+ around nucleic acids using the state-of-art MD simulations. We focused on the initial placement of ions in the simulations: type A simulations, using the initial locations of the ions determined by an electrostatic grid and type B, in which cations were randomly placed at least 35 Å away from nucleic acids. For Na + , the initial simulation conditions had little effect on the obtained number of bound ions around the duplex. The number of bound Na + ions reached equilibrium within 10 ns independent of the initial positions. For the Mg 2+ ion simulations, drastically different results were obtained depending on the initial locations of the cations. In type A simulations, many Mg 2+ ions were permanently bound to the DNA because they were placed initially too close to the nucleic acid and these ions did not form a fully solvated inner coordi-nation shell. On the other hand, in type B simulations when Mg 2+ ions were initially placed in the bulk water, the ions formed octahedral first solvation shell and this solvation shell was not disrupted throughout any of the simulations. Fully solvated Mg 2+ ions bind with DNA with events lasting as long as 20 ns. This observation alerts us to the importance of considering how to introduce the Mg 2+ ions into MD simulations. Since experiments have suggested that the residence time of water molecules in the first solvation shell of Mg 2+ ion lasts around a microsecond and removal of first solvating shell water molecules is energetically unfavorable, the direct ion coordination between Mg 2+ and DNA seen in MD simulation type A is likely an artifact of MD simulations. A recent computational study using replicaexchange MD simulations also reported that binding of "bare" magnesium ions to RNA atoms never occurred. (Kirmizialtin & Elber, 2010) In fact, they have found that crystallographically observed binding sites of magnesium ions, which were modeled from X-ray scattering data without first solvation shell, are well reproduced by their MD simulations results with magnesium ions that retain the first solvation shell. In addition, Petrov et al. used high-level quantum mechanics to calculate the potential binding modes of Mg 2+ ions with dimethyl phosphate and found that Mg 2+ ion prefers outer-sphere coordination with the phosphate moiety. (Petrov, Funseth-Smotzer, & Pack, 2005) Therefore, we conclude that the proper way to simulate Mg 2+ ions should begin with solvated Mg 2+ ions. Replica exchange simulations, such as those used by Kirmizialtin and Elber, are useful to obtain proper exchange between solvated Mg 2+ and desolvated Mg 2+ ions. (Kirmizialtin & Elber, 2010) However, since the exchange is very rare, regular MD simulations with solvated Mg 2+ ions will be sufficient to reproduce real experimental observations at room temperature.
Our MD simulations also show that the binding between Na + ions with DNA is principally through direct contact and these binding events last around several hundred picoseconds. These observations agree with earlier reports. For Mg 2+ ions, the water-mediated contact can be long lived. The longest binding event found in the simulation lasts up to 20 ns. These long-lasting contacts occurred in the groove region where solvated Mg 2+ ions make three water-mediated contacts with electronegative atoms. Groove binding was correlated to a strong sequence specific enhancement for both cations. Although, Mg 2+ ions did bind to the phosphate backbone with events lasting around 1-2 ns; these events were substantially shorter than groove binding and demonstrated no sequence dependence. The stable binding of solvated Mg 2+ ions with nucleic acids is correlated with the ejection of second shell water molecules. The more water- Figure 7. RDF of Na + to any closest atoms on DNA (excluding hydrogens) for Na_31A and Na31_08 MD systems. mediated contacts a solvated Mg 2+ ion can make with nucleic acids, the longer the binding event.