Molecular Dynamics Laboratory

Counterion - DNA interaction in water solution.
Atom-atom potential functions elaborated for Na+ and K+ ions.

Alexander V. Teplukhin
Institute of Mathematical Problems of Biology
Russian Academy of Sciences
Pushchino, Moscow region, 142292, Russia

George G. Malenkov
Institute of Physical Chemistry
Russian Academy of Sciences
31 Leninsky prosp., Moscow, 117915, Russia

Valery I. Poltev
Institute of Theoretical and Experimental Biophysics
Russian Academy of Sciences
Pushchino, Moscow region, 142292, Russia


In modern times, computer simulation occupies a prominent place in the variety of approaches and methods used to study molecular and atomic processes. However, all the programs that has enjoyed the widest application (AMBER(1), GROMOS(2), CHARMM(3) are best suited to investigation of a complex molecular systems (a numerous set of potential functions (PF) and their parameters is required to take properly all possible atomic pairwise combinations into account). Thus, the PF set adjusted for special group of molecules (e.g., nucleic acid fragments) seems to be more reliable and useful. By this strategy water-water and neutralized DNA-water potentials designed by us (4) in 1984 are successively employed to date (5,6). Recently this set of PF’s was revised to include new subsystems: charged phosphate and counterions. Below we represent the parameters of AAPF specially adapted to model ion-water, ion-base, ion-phosphate and water-phosphate interactions.

Counterion-DNA interactions both in water solution and in vivo are one of the important factors, determining local and global structure of DNA double helix. Therefore, potential functions using in computer simulations should be carefully tested for conformity to available experimental data. As in previous papers (4,5) we use two-stage routine for each type of newly included interactions. Firstly, bimolecular systems (one ion and one water molecule, one ion and one base, etc.) were studied using energy minimization and Monte Carlo techniques. This gives us an opportunity to monitor the interaction energy and interatomic distances while matching them to experimental ones. Secondly, the Monte Carlo simulations of ion hydration and, separately, phosphate hydration were carried out to evaluate the ion coordination number, hydration energies, hydrogen bond (H-bond) probabilities, etc. In addition, the ion lability in the neighborhood of the bases or phosphate anion was examined.


Starting data for counterion’s potential functions. PF (7) parameters for carbonyl O were adopted as initial for Na+ ion and those for aromatic C as initial ones for K+ ion. In the subsequent refinement only PF parameters describing the counterion-O and counterion-pyridine N interactions were changed.

Molecular models

Geometry of the base (9met-adenine, 9met-guanine, 1met-thymine, 1met-cytosine) and water molecules as well as their atomic charges and AAPF coefficients were as they had been described in the previous works (4,5). Dimethyl phosphate (DMP-) anion was chosen as a computationally convenient prototype for the study of the charged phosphodiester group of polynucleotides (8). For this study, DMP- is modeled by a rigid molecule in gg conformation (description and comments in Refs. 9,10). The particular geometry of DMP- anion was similar to that in corresponding fragment of B-DNA sugar-phosphate backbone studied in our previous work (6). The values of phosphate atomic charges (11) and Lennard-Jones parameters, published in Ref.7, were chosen as the initial ones to be refined in this work (see below).

Computational details

Except for the water-counterion pairs, an ordinary method for the potential energy of bimolecular system minimization was used. An azimuthal scanning of the potential energy surfaces was completed in a cylindrical reference system. Its origin coincides with geometrical center fixed nucleic acid base or DMP- and the polar axis is orthogonal to base or O-PO- plane. For this purpose, the energy minimization in 2D (for counterion) or 5D (for rigid water molecule) space was performed for a sequence (with 10° increment) of fixed azimuthal angles. In addition, a search for optimal counterion (or water) positions in the regions of local energy minima was made including the azimuthal variable. Water-counterion pairs and small water-counterion clusters were studied by Monte Carlo method.

Monte Carlo simulations of a single counterion (or DMP- anion) hydration in dilute solution at temperature 300° K were performed using the same procedure (Metropolis sampling, NVT ensemble) as in our earlier works (4,5). There were 256 or 400 water molecules included into counterionic or DMP- unit cells, respectively. The volume of a particular unit cell was assumed to be a sum of partial volumes of all constituents at standard conditions to provide the correct system density. The partial volumes for various components of our systems are as follows: water 30 A3, Na+ ion - 11A3, K+ ion 5 A3, DMP- anion 100 A3. Data for cations are taken from Ref.12 and were slightly changed by us (reported values, recalculated by us to A3, were - 11.69 A3 and 5.5 A3, respectively). The same value (- 11 A3) for Na+ ion was used in Ref.13. The partial volume of DMP- was found as a difference between those for Na+DMP- ion pair (53.1 mL/mol (8)) and Na+ ion. The ion under consideration was fixed in the center of cubic simulation cell. Periodic boundary conditions (14) were applied and the interactions were calculated according to the minimum image prescription (14). To achieve equilibrium the initial segment of Markov’s chain, consisting of no less than 50000 trials per water molecule, was excluded from consideration. The averaging was carried out over the chains up to 2 million (cation-256 waters) or 500000 (DMP--400 waters) trials per water molecule in length.

Dependencies of various energetical characteristics of systems under consideration on water-ion or solute-ion separation distances (Figures 1, 3, 5,7) were plotted using mean values obtained via MC averaging over short (0.02 A for water-ion systems and 0.05 A for all others) incremental intervals. This method is conceptually identical with the one used for calculation of radial distribution functions.

In addition, something of the methodological aspects, associated with a particular system, are given in correspondent paragraphs of Results section.


Water-counterion systems

Several variants of atom-atom potential functions and coefficients were tested. Those of them, finally adopted by us, are described below. The potential function employed for Na+-water O interaction is of 1-12-10 format: Uij=qiqj/rij+Bij/rij12-Aij/rij10 and all other (K+-water O, ion-water H) are of 1-12-6 format: Uij=qiqj/rij+Bij/rij12-Aij/rij6.Here qi is the charge associated with a particle i, and B and A define the van der Waals interaction. The values for A and B coefficients are listed in Table V.

Water-ion pairs and small water-ion clusters. These simple systems, having all water molecules directly contacting with ion, are excellent objects for primary testing of ion-water potentials for conformity to experimental data. A selection of energetical and structural characteristics calculated by us using Monte Carlo method (Metropolis sampling at 300° K in cluster approximation) as well as experimental (15) and ab initio calculations (16) data are given in Table I. As it evident from the Table, our results are in rather good agreement with experimental and other calculated data, even though the particular PF’s were elaborated for the water molecule model with fixed geometry and charge distribution specially designed (4) for a bulk medium simulation.

Table I
Characteristics of water ion clusters M+(H2O)n

n ΔE ΔH ΔH298exp RMO RMOab initio
1 -20.93 -21.53 -24.0 2.39 2.21
2 -19.85 -20.45 -19.8 2.39 2.25
3 -18.10 -18.70 -15.8 2.40 2.29
4 -16.17 -16.77 -13.8 2.41 2.31
5 -13.30 -13.90 -12.3 2.42 2.32
6 -10.80 -11.40 -10.7 2.43 2.42
1 -14.51 -15.11 -17.9 2.77 2.58
2 -13.69 -14.29 -16.1 2.77 2.68
3 -12.54 -13.14 -13.2 2.78 2.70
4 -11.09 -11.69 -11.8 2.80 2.72
5 -9.72 -10.32 -10.7 2.81 2.73
6 -9.44 -10.04 -10.0 2.82 2.89

Incremental energy Δ E and enthalpy Δ H changes (kcal/mol) at 300° K for the reaction M+(H2O)n-1 +H2O → M+(H2O)n. Δ H values were calculated by subtraction of Δ A=pΔ V@ 0.6 kcal/mol (ideal work of expansion of 1 mol of water molecule gas at 300° > K) from the corresponding Δ E values. Δ H298exp taken from Ref.15. RMO and RMOab initio are the average ion-water O distances (in A) estimated, respectively, from Monte Carlo simulation and from ab initio calculation data (16).

Ion hydration in dilute solution. Systems containing 256 water molecules and one ion were confined in cubic unit cells with the edge size 19.72 A (Na+) or 19.734 A (K+). The most important data, obtained from the Monte Carlo simulation are presented in Figure 1 and Table II. Ion-water radial distribution functions (Figure 1a) both Na+ and K+ ions have well-pronounced first and second peaks, and the same is true for corresponding minima of Etotal dependencies (Figure 1b). First maxima these RDF’s are positioned at 2.42 A and 2.84 A for Na+ and K+ ions, respectively. These are in good agreement with values 2.38 A and 2.8 A (17) obtained for ion-water O distances distributions calculated from structural data on crystal hydrates of Na+ and K+ salts. However, owing to the interactions between the particular crystal ingredients, the ion-water distances are distributed in a wide range (for example, from 2.3 A to more than 2.55 A for Na+). So, the average ion-water distances for Na+ and K+ ions in these crystals are 2.42 A and 2.87 A (17), respectively. The information on water-ion distances as well as calculated dependencies of total potential energy of individual water molecule (Etotal in Figure 1b) on these distances were used to work out the effective way to determine coordinational number of ion (CN). We treated CN of ion as an average number of water molecules located around the ion in the sphere of radius RCN. Practically, we used RCN 3.1 A for Na+ and 3.4 A for K+, which are equal to water-ion distances corresponding to the first maxima of Etotal dependencies. Calculated CN values (Table II) are in good agreement with experimental data obtained using the similar ideology (17,18).

Experimentally measured values for the ion hydration enthalpy are typically fall in rather wide range (see Table II) and are strongly dependent on particular method of solving the sharing problem. The ion hydration energy Δ Ehyd presented in Table II is found as a difference the total potential energies of water systems calculated by Monte Carlo method with and without ion. As seen from the Table II, Δ Ehyd of K+ ion is in excellent agreement with experimental, but Na+ hydration energy seems somewhat below the experimental values.

Table II
Hydration characteristics of Na+ and K+ ions

Ion Etotal Ewi Δ Ehyd Δ Hhyd CN
Na+ -2610 -188 -109 -101—106 6.68
K+ -2584 -151 -83 -81—86 7.67

Presented are the values of total potential energy Etotal and water-ion contribution Ewi to it as well as ion hydration energies Δ Ehyd (all in kcal/mol) and the ion coordination numbers CN. Δ Hhyd (kcal/mol) - experimental enthalpy of hydration at 298° K (19,18). An ideal work 0.6 kcal/mol of 1 mol gas expansion should be added to &Delat; Hhyd while comparing with Δ Ehyd.

Worthy of mention are the dependencies of Etotal on Rio distances. As it seen from Figure 1b, the total potential energies of individual water molecules forming the first hydration shell around Na+ ion (first minimum of Etotal curve) are by 2.5 kcal lower than the energies in pure water. At the same time, water molecules in first hydration shell of K+ ion have only a slight gain in comparison with bulk water. Such a difference in hydrational characteristics of Na+ and K+ ions is in accordance with the Samoilov’s concept (20) on positive (for Na+) and negative (for K+) hydration.

Water-phosphate interaction

Negative DMP- charges are slightly rearranged between ester and anionic oxygens in comparison with original (11). This was made to diminish the phosphate oxygen charge asymmetry. The DMP- atomic charge distribution finally adopted by us is as follow: P +0.802, anionic O –0.730, –0.313, C –0.039, H +0.060 - all values in |e| units. The PF’s employed for water H-phosphate O interactions are of 1-12-10 format and all other are of 1-12-6 format. The parameters of ester oxygens were the same as for hydroxyl O in Ref.5 and those for anionic O were newly developed. Final data on the PF coefficients are in Table V.

Water-DMP- bimolecular complex. Global minimum of the interaction energy (-19.04 kcal/mol) corresponds to water bridge formation between anionic oxygens (position 1 in Figure 2; ROwO- is 2.72 A). Other two local minima corresponding to the bridging between ester and anionic oxygens (positions 2 and 3 in Figure2; ROwO- is about 2.68 A and ROwOest is about 2.84 A) are by 3 kcal/mol higher. At the same time, the shortest DMP-- water distances corresponds to the formation of linear (slightly bent) H-bond with phosphate oxygens (positions 4 and 5 in Figure 2; Hw…O- –1.76 A, Ow…O- –2.71 A, and 1.95 A and 2.79 A for ester O; respective energies are -15.5 and -13.2 kcal/mol). This is in accordance with other calculation results of other authors (21,22), see Table III.

Table III
Characteristics of DMP- -water interactions in vacuo
Site This work

Ab initio

Mol. mech.

Mol. mech.

  E RHwO- RHwOest E RHwO- E RHwO- E**
1 -19.04 1.90 - -22.19 2.10 -19.40 1.87 -19.0
2 -16.05 1.83 2.06 nd nd -16.0
3 -16.00 1.83 2.06 nd nd -17.0
4 -15.51 1.76 - -18.82* 1.75 -14.55 1.72 nd
5 -13.22 - 1.95 nd nd nd

Site numbers as in Figure 2. E DMP-- water interaction energy (kcal/mol); RHwO- and RHwOest - distance (in A) between water H and anionic or ether O of DMP- (data for water O are in text of the paper); nd - no data; (*) - average for two values (-19.68 and -17.95) presented in Ref. 21; (**) data for H2PO4.

DMP- hydration in dilute solution. System containing 400 water molecules and DMP- anion was confined in a cubic cell with the edge size of 22.956 A. Partial Ow-Ophos and Hw-Ophos RDF’s are depicted in Figure 3 (a, b). Those, corresponding to anionic O have clear-cut first maxima at 2.75 A and 1.79 A for Ow…O- and Hw…O- distances. The number of H-bonds between waters and anionic O calculated by Monte Carlo method using geometrical criterion (4) is, on average 3.3. As to H-bonds between waters and ester O they form much more rarely (on average 0.9-1.2 per oxygen) and their geometrical characteristics are distributed over a wide range (Ow…O and Hw…O maxima are near 2.85 A and 1.92 A, see Figure 3a,b). This is in excellent agreement with experimental data on water rich crystal hydrates of nucleotides and their analogs. For example, as has been observed in the crystal structure of 5'-dCMP•Na2•7H2O (23) two of three anionic O’s form three H-bonds and third forms two H-bonds with water molecules. In this crystal mean Ow…O- distance was 2.77± 0.06 A and Hw…O- was about 1.84±0.1 A.

Multilayer structure of water, surrounding DMP- , is clearly seen in Figure 3c, where dependency of local water density (referenced to bulk water) as well as its running sum (Figure 3d) on distance between water oxygen and P atom of DMP- are presented. First shell consists of 10 water molecules confined in a sphere 4.45 A in radius, but only about 8.7 of them are H-bonded to phosphate O. Dependency of the total potential energy of individual water molecule on Ow…P distance (Figure 3c) has well pronounced first minimum (3.4 kcal/mol below the bulk water level) corresponding to water-phosphate H-bonds formation. So, DMP- might be treated as the local water structure stabilizing (positively hydrated (20)) anion.

Main structural features of DMP- water shell are the water bridges forming between phosphate oxygens. The most probable of them are as follows: two-water bridges O1-…O2- (found in 57% configurations of Markov’s chain) and Oester…O- (25%); three-water bridges - 100% (this value means that anionic oxygens might be connected via two different three-water bridges simultaneously) and 37%, respectively. Single-water bridges as well as bifurcated H-bonds occur very rarely (probability less then 3%).

Calculated hydration energy (the difference between the total interaction energies of systems with and without anion: -4001.3 and -3908 kcal/mol, respectively; water -- DMP- interaction energy is -181.7 kcal/mol) of DMP- is -93.3 kcal/mol. It is in good agreement with value -92.5 kcal/mol obtained from Monte Carlo calculation of DMP- hydration in NPT ensemble (9).

Counterion-phosphate interaction

The PF’s employed for Na+-phosphate O interactions are of 1-12-10 format and all others (including those for K+-phosphate O) are of 1-12-6 format. Only parameters of counterion-oxygen PF’s were refined. Final data on the PF coefficients are in Table V.

Counterion-DMP -pairs in vacuo. Global minimum of the interaction energy (-139.7 kcal/mol for Na+ and -120.3 kcal/mol for K+) corresponds to O-PO- bisector position (site 1, Figure 4), while the ion 'bridges' between ester and anionic oxygens (sites 3,4, Figure 4) are unstable (energy values obtained by azimuthal scanning are higher the minimal by about 26 kcal/mol for Na+ and 19 kcal/mol for K+). It should be noted that in crystal hydrates these sites have only a slight occupancy (especially for Na+). Anionic O-counterion distances in the global minima are 2.39 A and 2.79 A for Na+ and K+, respectively. Closest distances between anionic O and ion are 2.31 A for Na+ and 2.6 A for K+ and correspond to the ion position in cone coaxial with PO- bond (site 2, Figure 4). The distances, corresponding to the ion bridging between ester and anionic O (sites 3,4, Figure 4) were obtained by azimuthal scanning. The closest are in range 2.38-2.52 A for Na+ (so, if ester O…Na+ distance is 2.38 A, respective value for anionic O…Na+ is 2.52 A and reversibly) and 2.76 A for K+. Calculated geometrical parameters are in excellent agreement with those from crystal hydrates data. For example, distances between anionic O and ion are distributed in region 2.30—2.46 A for Na+ (24-27) and 2.7—3.2 A for K+ (28,29), those for ester O are 2.5—2.58 A for Na+ (26) and 2.84 A for K+ (30).

Counterion-DMP-pairs in dilute solution. System containing 400 water molecules, DMP- anion and counterion was confined in a cubic cell with the edge size 22.95 A for Na+ and 22.96 A for K+. The position of DMP- was fixed in the cube center and O-PO-plane was set parallel to the cube face. The counterion maximal shift value 0.065 A was chosen to obtain a 48 % acceptance rate for Monte Carlo trials (for comparison, maximal shift and rotational angle of water molecules were about 0.122 A or radians, respectively). Periodical boundary conditions with minimum image prescriptions were applied. Ion-ion interactions between the unit cell neighbors were neglected.

Initial position of counterion was chosen near the global minimum interaction energy of ion-DMP- pair in vacuo. First, MC run of 50000 trials per water molecule was completed for the system with fixed geometry of ion-DMP- pair. After this, counterion movement was included in MC procedure. MC run was stopped when ion-phosphorus atom separation becomes more than 10.5 A. Such a dissociation occurs quickly for K+ ion (150-200 thousand trials per water molecule) in contrast to Na+ where 2 or 3 times more configuration is needed. Described routine was repeated for each counterion (starting from newly generated initial configuration of water molecules) to receive more reliable data.

Space of the simulation cell was divided by 3D grid with elementary cell volume of 0.5*0.5*1.5 A3 (x, y, z, respectively; Z-axis was orthogonal to O-PO- plane and X-axis coincides with -PO- bisector line) to trace dissociation pathway of counterion. So doing, PO2-group was situated in the middle of central Z-slice. Then probability of counterion location in each of these elementary cells was estimated (as ratio Ni/N, where Ni is the number of MC configurations in which ion was situated in i-th cell and N is the total number of MC configurations for the particular MC run). As it follows from the analysis of these probability maps, dissociation pathways have many common features. So, starting from initial position near site 1 (see Figure 4) counterion gradually (several dozen thousands of MC trials) moves out of -PO- plane by about 1.5 A and shifts by about 1 A along Y-axis breaking van der Waals contact with one of anionic O. This configuration (site 2 in Figure 4) is relatively stable, because it takes 100 or 150 thousand MC trials to either continue dissociation or (that more probable in case of Na+) to jump to other anionic O. Further, dissociation path is confined in a slightly curved pipe (diameter is about 1.5 A for Na+ and up to 2.5 A for K+) passing parallel to O-PO- bisector line. As it follows from the probability map data, there are at least two attractive regions near each anionic O. First one corresponds to van der Waals contact of counterion with anionic oxygen. Second, less pronounced, is located at distance about 4.5 A from anionic O and might be associated with water bridged arrangement of counterion-DMP- pair. Nevertheless, there are no any remarkable peculiarities of energetical characteristics (see Figure 5 and Ref. 31) corresponding to second region.

As it was found from the calculations, the total energy of dissociated complex is approximately by 25 kcal/mol lower than that of associated one. The dependencies of Etotal on P-ion distance Rpi as well as Eidmp, Ewdmp and Eiw contributions are depicted in Figure 5. As it is evident from this Figure, the dependencies of Ewdmp and Eiw on Rpi become practically constant for separation distances more than 9.5 A, where its values are the same as for independently hydrated DMP- and cation (see previous paragraphs). Standard deviation of instant values of Etotal is about 32-35 kcal/mol and might be directly associated with thermal fluctuations and correspond to constant volume heat capacity about 20-23 cal/mol/deg (experimental value 21 cal/mol/deg (18)). Nevertheless, standard deviation of Etotal mean values obtained by averaging over Rpi microintervals (0.05 A) and plotted in Figure 5 is about 10 kcal/mol (linear fit).

Counterion-nucleic acid base interaction.

The PF’s employed for Na+-carbonyl O and Na+-pyridine N interactions are of 1-12-10 format and all the others are of 1-12-6 format. Only parameters of counterion-O and counterion-pyridine N PF’s were refined. Final data on the PF coefficient are in Table V.

Na+-base and K+-base interactions in vacuo. Positions of local minima of the ion-base interaction energy are shown in Figure 6. The energy values and respective ion-atom distances are given in Table IV. Calculated geometrical characteristics are in good agreement with experimental distances obtained from X-ray data on crystal hydrates: Na+…carbonyl O - 2.43 (2.34-2.56) A (24,26), Na+…pyridine N - 2.415 and 2.613 A (32), K+…carbonyl O - 2.77-2.79 A (29), K+…pyridine N - 3.15 A (28). At the same time, binding energies obtained via ab initio calculations (33,34) are systematically lower than ours values and correspond to shorter ion - base O or N distances: 2.0-2.15 A (32) and Na+…N - 2.3-2.41 A, Na+…O - 2.26 A, K+…N - 2.79-2.86 A, K+…O - 2.67 A (34). This is not surprising, because our PF’s were elaborated to simulate interactions into bulk media, where such short ion base contacts are practically absent. Moreover, several of ion-base complexes (sites near NH2-group) studied by us are non-planar: see Z in Table IV. It should be noted, what K+ - base atom distances, calculated via ab initio (34), are in good agreement with both x-ray and our data.

Table IV
Characteristics of base-cation interactions in vacuo


9m-adenine 9m-guanine 1m-cytosine 1m-thymine


1 2 3 1 2 1 1 2
E -18.02 -15.16 -14.28 -45.30 -8.64 -46.36 -19.95 -14.94
E1ai -26.4 -24.0 -21.3 -53.9 -14.6 -51.7 -32.9 -28.7
E2ai nd nd -26.5* -59.7* nd nd nd nd
RIO -- -- -- 2.42 -- 2.40 2.37 2.39
RIN 2.46 2.47 2.46 2.46 2.46 2.47 -- --
Z 0.00 0.05 0.37 0.00 0.38 0.00 0.00 0.00
E -12.66 -10.24 -8.64 -35.09 -3.51 -35.65 -14.78 -10.49
E2ai nd nd -15.0* -42.95* nd nd nd nd
RIO -- -- -- 2.76 -- 2.71 2.73 2.77
RIN 2.85 2.87 2.84 2.81 2.89 2.84 -- --
Z 0.12 0.01 0.69 0.00 1.93 0.00 0.00 0.00

Site numbers are as in Figure 6. E - base-ion interaction energy (kcal/mol) and E1ai - data obtained via ab initio calculation of Na+ interaction with adenine, guanine, cytosine and uracil (taken from Ref. 33) and E2ai - ab initio energy of ion interaction with adenine and guanine (taken from Ref. 34); nd - no data; (*) - mean value for two variants presented in Ref. 34. RIO and RIN - distance (in A) between ion and base O or N (see Figure 6, dotted lines), respectively. Z - distance (in A) between ion and base plane.

Ion-base pairs in water solution. System, containing 400 water molecules, base and ion, was confined in cubic cell with edge size 23.015 A for Na+ or 23.017 A for K+. Position of base was fixed in the cube center and base plane was set parallel to pair of cube faces. So, z-axis was set orthogonal to base plane and x-axis was set parallel to long axis of base. Periodical boundary conditions with minimum image prescriptions were applied. Ion-ion interactions between the unit cell neighbors were neglected.

Molecule of 9met-guanine was chosen to investigate ion-base pair dissociation pathway and ion mobility near the base, because this ion-base complex is the most stable in vacuo (respective site of cytosine is inaccessible in DNA duplex). Initial position of ion was chosen in the vicinity of the global minimum interaction energy of ion-base pair in vacuo (site 1, Figure 6, G). MC run of 50000 trials per water molecule was completed on system with fixed geometry of ion-base pair. After this, counterion movement was included in MC procedure. MC run was terminated when ion-N7 guanine atom separation becomes more than 10.5 A.

The same, as in case of counterion-DMP- study, routine was applied to trace ion dissociation pathway. In such a way, the base was situated in the middle of central z-slice. Again, as in case of counterion-DMP- dissociation, ion firstly moves out of base plane by about 1.5 A breaking van der Waals contact with O6 atom of guanine. Estimated mean residence time in this region is about 100000 and 50000 MC trials for Na+ and K+ respectively. On breaking van der Waals contact with N7 (ion - base N atom distance within 3.6 A) the ion usually shifts to O6 and makes contact with oxygen atom. Residence time estimated for this site is in equivalence to one or two dozen thousand MC trials per ion/water. After two or three such jumps the ion moves to secondary less pronounced attractive region at distance about 4.5 A from N7 or O6 (2-4 A out of base plane). Then ion randomly walks along N7-O6 edge of guanine. It takes one (K+) or two (Na+) hundred thousand MC cycles to leave this region.

There are no statistically reliable data on final portion of dissociation path, yet. We carried out three MC runs: one for Na+ ion and two for K+ ion. In two cases (one for Na+ ion - totally 1000000 MC trials per ion/water and one for K+ ion - totally 300000 MC trials per ion/water) ion dissociation movement was in direction from site 1 (Figure 6, G) to C8 side of guanine and in one case (for K+ ion - totally 500000 MC trials per ion/water) in direction to N1 side of guanine. In last case K+ moved out of base plane by 6 A (in contrast to 2-4 A in first two cases). Finally, owing to periodical boundary conditions, all three MC runs finished in the same region (±2 A) which is situated near the middle of line (collinear to Watson-Crick base pair long axis) connecting C8 atoms of guanines in adjacent unit cells. It should be noted, that final position of anion might be strongly dependent on particular choice of simulation cell geometry as well as PF set. Thus, extensive study is needed to obtain more data.

As it was found from the calculations, the total energy of dissociated complex (including water-water contribution), Etotal, is roughly by 5 (for K+ ) or 10 (for Na+ ) kcal/mol lower than of associated one, see Figure 7. This effect (estimated from linear fit) is rather slight in comparison with standard deviation (about 10 kcal/mol) of Etotal mean values obtained over 0.05 A RN7i microintervals (especially for K+). However, it may illustrate dissociation tendency of guanine-ion pair hydration.

The dependencies of Etotal on distance RN7I between the ion and N7 atom of 9met-guanine as well as ion-base (Eig), water-ion (Eiw) and water-guanine (Ewg) contributions are presented in Figure 7. A step change in Eig near 4.5 A is due to strengthening of ion-O6 atom interactions while ion moves from N7 atom to N1 side of guanine (this occurs then the ion switches van der Waals contact from N7 to O6). The dependency of water-ion interaction energy Eiw on N7-ion separation distance becomes practically constant after 8.5 A. The same is true for water-guanine interaction energy Ewg.

Table V
Coefficients of potential functions

  Na+ K+
Atoms A B A B
H1 121 17800 126 27300
H2 121 17800 126 61700
H3 121 42200 126 81600
C1 305 349000 316 601000
C2 385 406000 400 704000
N1 327 317000 334 544000
N2 391 366000 400 630000
N3 190000* 1095000 2147 1153000
O1 129000* 688000 1674 701000
O2 129000* 688000 1674 701000
O3 169000* 954000 3220 1730000
P 1182 1352000 1156 2195000
OW 150840* 830200 1954 954000

  Ow Hw(H1)
  A B A B
O3 326 531400 6225* 15900

Atom types: H1, H2, and H3 are hydrogen atoms attached to N or O, to aromatic C, and to aliphatic C, respectively; C1 and C2 are aliphatic and aromatic carbons; N1 is amino-group nitrogen, N2 and N3 are pyrol and pyridine nitrogens, respectively; O1 is ester oxygen, O2 and O3 are carbonyl and anionic oxygens, respectively; P is phosphorous atom; Ow is water oxygen. Values marked with (*) are utilized in PF’s of 1-10-12 format.


  1. Weiner, S.J., Kollman, P.A., Nguyen, D.T. and Case, D.A., J. Comput. Chem. 7, 230-252 (1986).

  2. van Gunsteren, W.F. and Berendsen, H.J.C., GROMOS86 Groningen Molecular Simulation Systems; Groningen University, Groningen, 1986.

  3. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S. and Karplus, M., J. Comput. Chem. 4, 187-217 (1983).

  4. Poltev, V.I., Grokhlina, T.I. and Malenkov, G.G., J. Biomol. Struct. Dyn. 2, 413-429 (1984).

  5. Poltev, V.I., Malenkov, G.G., Gonzalez, E.J., Teplukhin, A.V., Rein, R., Shibata, M. and Miller, J.H., J. Biomol. Struct. Dyn. 13, 717-725 (1996).

  6. Teplukhin, A.V., Zhurkin, V.B., Jernigan, R. and Poltev, V.I., Mol. Biol. (Moscow) 30, 121-135 (1996), Engl. Transl. Molecular Biology 30, pt.2, 75-84 (1996).

  7. Zhurkin, V.B., Poltev, V.I. and Florentiev, V.L., Mol. Biol. (Moscow) 14, 1116-1129 (1980), Engl. Transl. Molecular Biology 14, 882-895 (1980).

  8. Jayaram, B., Mezei, M. and Beveridge, D.L., J. Comput. Chem. 8, 917-941 (1987).

  9. Alagona, G., Ghio, C. and Kollman, P.A., J. Am. Chem. Soc. 107, 2229-2239 (1985).

  10. Huston, Sh.E. and Rossky, P.J., J. Phys. Chem. 93, 7888-7895 (1989).

  11. Sasisekharan, V. and Lakshminarayanan, A.V., Biopolymers 8, 505-514 (1969).

  12. Desnoyers, J.E. and Jolicoeur, C., In Modern aspects of electrochemistry, J. O’M. Bockris and B.E. Conway Eds., Plenum Press, N.Y., 1969.

  13. Chandrasekhar, J. and Jorgensen, W.L., J. Chem. Phys. 77, 5080-5089 (1982).

  14. Allen, M.P. and Tildesley, D.J., Computer simulation of liquids, Oxford University Press, N.Y., 1987.

  15. Dzidic, I. and Kebarle, P., J. Phys. Chem. 74, 1466-1480 (1970).

  16. Feller, D., Glendening, E.D., Woon, D.E. and Feyereisen, M.W., J. Chem. Phys. 103, 3526-3542 (1995).

  17. Malenkov, G.G., In The chemical physics of solvation, R. Dogonadze et al. Eds., pt. A, Elsevier, Amsterdam etc., pp.355-389 (1985).

  18. Friedman, H.L. and Krishnan, C.V., In Water: A comprehensive treatise. F. Franks Ed., Plenum Press, N.Y., v.3, pp.1-118 (1973).

  19. Vasil’ev, V.P., Zolotarev, E.K., Kapustinskii, A.F., Mischenko, K.P., Podgornaya, E.A. and Yatsimirskii, K.B., Zh. Fiz. Khim. (Moscow) 34, 1763-1767 (1960). Russian J. Phys. Chem.

  20. Samoilov, O.Ya., Structure of electrolyte solutions and the hydration of ions, (Engl. Transl.), Consultants Bureau Enterprises, Inc., N.Y., 1965.

  21. Alagona, G., Ghio, C. and Kollman, P.A., J. Am. Chem. Soc. 105, 5226-5230 (1983).

  22. di Oliveira, M., J. Comput. Chem. 7, 617-628 (1986).

  23. Pandit, J., Seshadri, T.P. and Viswamitra, M.A., Acta Cryst. C39, 342-345 (1983).

  24. Young, D.W., Tollin, P., and Wilson, H.R., Acta Cryst. B30, 2012-2018 (1974).

  25. Camerman, N., Fawcett, J.K. and Camerman, A., J. Mol. Biol. 107, 601-621 (1976).

  26. Seeman, N.C., Rosenberg, J.M., Suddath, F.L., Kim, J.J.P. and Rich, A., J. Mol. Biol. 104, 109-144 (1976).

  27. Rosenberg, J.M., Seeman, N.C., Day, R.O. and Rich, A., J. Mol. Biol. 104, 145-167 (1976).

  28. Adamiak, D.A. and Saenger, W., Acta Cryst. B36, 2585-2589 (1980).

  29. Emerson, J. and Sundaralingam, M., Acta Cryst. B36, 537-543 (1980).

  30. Viswamitra, M.A., Hosur, M.V. and Katti, S.K., In Conformation in Biology, the Festschrift celebrating the sixtieth birthday of G.N. Ramachandran F.R.S. R. Srinivasan and R.H. Sarma Eds., Adenine Press, N.Y., pp.439-444 (1983).

  31. Friedman, R.A. and Mezei, M., J. Chem. Phys. 102, 419-426 (1995).

  32. Katti, S.K., Seshadri, T.P. and Viswamitra, M.A., Curr. Sci. 49, 533-535 (1980).

  33. Perahia, D., Pullman, A. and Pullman, B., Theoret. Chim. Acta (Berl.) 43, 215-225 (1977).

  34. Burda, J.V., Sponer, J. and Hobza, P., J. Phys. Chem. 100, 7250-7255 (1996).