INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Water plays an important role in modulating the structural and functional properties of biological macromolecules, and its interaction with both proteins and nucleic acids is the subject of an intense investigation from both an experimental and a simulative approach. The surface hydration of proteins depends on the geometry and quality of surface groups: waters within polar cavities exchange more slowly and can be detected by x-ray diffraction, while nonpolar cavities do not usually contain long-lasting water molecules (Schwabe, 1997; Kuhn et al., 1992; Luise et al., 2000). In the case of nucleic acids it has been shown that their hydration is directly related to their conformation (Westhof, 1988, 1993; Rentzeperis et al., 1993; Chalikian et al., 1994; Jayaram and Beveridge, 1996; Cheatham and Kollman, 1997a). Information on the water position around a DNA sequence is provided by x-ray diffraction, which reveals the favored average position occupied by a water molecule. The observation of highly localized water molecules forming a "spine of hydration" in the minor groove of the B-DNA dodecamer (CGCGAATTCGCG) (Drew and Dickerson, 1981; Kopka et al., 1983; Subramanian et al., 1988, 1990) was successfully reproduced using static DNA models (Subramanian and Beveridge, 1989; Berman, 1994) and later by using a dynamic model for several DNA sequences (Duan et al., 1997; Young et al., 1997; Feig and Pettitt, 1999).

A network of water molecules bridging the side-chain protein atoms and the DNA bases have also been observed in the x-ray diffraction pattern of protein-DNA systems, such as the complex between trp repressor and trp operator or between the BamHI restriction enzyme and its specific DNA sequence (Shakked et al., 1994; Newman et al., 1994, 1995). These and other related works have raised the question of whether interfacial water molecules play an important role in mediating sequence-specific recognition and in enhancing the protein affinity for specific DNA sequences (Schwabe, 1997).

An accurate structural and dynamical description of DNA in aqueous solution, helping to understand the molecular basis of protein-DNA and drug-DNA recognition, can be obtained through molecular dynamics simulations. Realistic simulations of DNA molecules are currently possible because of three methodological improvements: the inclusion of hydration and ionic patterns in their environment, the introduction of Ewald summation technique (Cheatham et al., 1995, 1997; Young and Beveridge, 1998) to account for the undesired effects of the truncation method, and the use of an accurate force field. The dynamical properties of the water molecules around the trp operator have been recently studied by MD and NMR spectroscopy (Bonvin et al., 1998; Sunnerhagen et al., 1998). Both the simulative and experimental study indicated that water molecules observed by x-ray diffraction are not characterized by long residence time.

Here we make use of the most recent Amber force field (Cornell et al., 1995) and of the particle mesh method (Darden et al., 1993) to investigate through MD simulation the dynamical properties of water around the DNA duplex interacting with the BamHI endonuclease. The results indicate that hydration is characterized by relatively short water residence times at any DNA site. However, evidence is presented for a spine of hydration in the DNA minor groove, indicating that sites having a high probability of being hydrated do not necessarily trap water molecules for a long time.

	COMPUTATIONAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

System setup

The canonical B-DNA structure (Arnott et al., 1976) was generated for the DNA sequence d(TATGGATCCATA)₂ using the nucgen module included in the AMBER 5.0 software package. The B conformation was chosen because x-ray diffraction has shown that in the DNA-BamHI endonuclease complex (Newman et al., 1995) the DNA duplex retains a B-like conformation.

The DNA macromolecule was immersed in a rectangular box filled with TIP3P (Jorgensen et al., 1981) imposing a minimum solute-wall box distance of 10 Å. The system was neutralized with Na⁺ cations using the AMBER leap module. The total system is composed of 762 DNA atoms, 22 Na⁺ counterions, and 3376 water molecules, giving a total of 10,912 atoms.

Dynamics simulation protocol

The system was modeled using the AMBER95 all-atom force field (Cornell et al., 1995) with periodic boundary conditions. A cutoff radius of 9 Å was used for nonbonded interactions, updating the neighbor pair list every 10 steps. The electrostatic interactions were calculated with the Particle Mesh Ewald method (Darden et al., 1993; Essmann et al., 1995). The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain all bond lengths involving hydrogen. Optimization and relaxation of solvent and ions were initially performed while keeping the solute atoms constrained to their initial position with decreasing force constants of 500, 25, 15, and 5 kcal/mol Ų. Thereafter the system was minimized without any constraints, and warmed up. A 3-ns simulation was carried out at constant temperature (Berendsen et al., 1984) of 300 K and at a constant pressure of one atmosphere with a 2-fs time step. Pressure and temperature coupling constants were 0.5 ps.

Analysis

During the production run, from 0.5 to 3 ns, the atom coordinates were saved every 0.1 ps for analysis. The AMBER carnal module was used to check some structural properties (root-mean-square deviation (rmsd), hydrogen bond). The hydrogen bond criterion was a maximum donor-acceptor distance of 3.5 Å and a minimum donor-proton-acceptor angle of 120°. The mean and maximum water residence times and the coordination number were calculated using the algorithm described by Garcia and Stiller (1993), considering a radius sphere of 3.5 Å around DNA sites and a time resolution of 0.1 or 10 ps, respectively. The solvent-accessible surface area of DNA sites was evaluated with the program NACCESS (Hubbard and Thornton, 1993), using a default probe size of 1.4 Å. DNA geometry parameters were calculated with CURVES (Lavery and Sklenar, 1988, 1989).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Evaluation of the energetic contributions

A molecular dynamic simulation of the DNA oligonucleotide d(TATGGATCCATA)₂, entirely solvated and neutralized with Na⁺ ions, was performed for 3 ns using the AMBER95 all-atom force field. Fig. 1, A and B show the electrostatic and van der Waals energies versus simulation time of DNA-DNA, DNA-ions, and DNA-water interaction, respectively. Both the electrostatic and van der Waals potential energies for the DNA-DNA, DNA-water, and DNA-ions interactions were all fairly stable over the 3 ns of simulation. The equilibration of the system is dominated by the DNA-solvent and DNA-ions electrostatic interactions, which show relatively large fluctuations, greater than those observed in a previous MD simulation, carried out with the Gromos force field for a DNA sequence corresponding to the Trp operator (Bonvin et al., 1998). This may be due to the fact that the partial charges attributed to the atoms in the Gromos force field are lower than those attributed in AMBER, which was used in this work (Cornell et al., 1995). The DNA-DNA electrostatic interaction is almost negligible over all the trajectory. The van der Waals potential energies are negligible for what concerns the DNA-ion interaction, while the DNA-water and mainly the DNA-DNA interactions have an important role for the internal stability of the double helix.

View larger version (36K): [in this window] [in a new window]   FIGURE 1   Non-bonded energies as a function of simulation time of the DNA duplex interacting with BamHI. (A) Electrostatic energies; (B) van der Waals energies.

Because of the above consideration the system can be considered fully equilibrated after 300-400 ps. To have a reliable sampling of our system, the analysis concerning its structural and hydration properties have been carried out for the last 2.5 ns of our trajectory, i.e., the time range between 0.5 and 3.0 ns.

Structural analysis

The stability of our system during the MD trajectory has been checked by monitoring specific structural parameters as a function of time. The rmsd calculated for all atoms from the B-DNA, the A-DNA, the average molecular dynamics, and the crystallographic x-ray structure (PDB entry 2bamhi), excluding the first two and the last two basepairs of the DNA tract, are reported in Fig. 2, A and B, respectively. The deviation from B-DNA increases with time, whereas it decreases from A-DNA (Fig. 2 A), although over all the trajectory the simulated structure remains closer to that of B-DNA. The rmsd from the average molecular dynamics structure (Fig. 2 B, bottom curve) oscillates around 1 Å, indicating that the DNA settles in a well-defined and stable configuration during the simulation. The rmsd from the crystallographic x-ray structure of the DNA segment interacting with the BamHI endonuclease (Fig. 2 B, top curve) is less than that from the B-DNA structure. Because the starting structure of the simulation is the B-DNA, the system spontaneously evolves toward the crystallographic structure during the thermalization and oscillates around this configuration. This result is in agreement with the x-ray evidence that in the DNA-BamHI complex the DNA structure remains quite rigid, without any major bends or kinks, while the protein undergoes a series of conformational changes (Newman et al., 1995). Up to now the 3D structure of the DNA segment in absence of the BamHI endonuclease is not available, preventing any comparison with our simulation.

View larger version (30K): [in this window] [in a new window]   FIGURE 2   Root-mean-square atom positional deviations as a function of simulation time. (A) From A-DNA (top curve), from B-DNA (bottom curve); (B) from the average MD structure of the DNA in solution (bottom curve) and from the x-ray diffraction structure (top curve). The structures were superimposed and the rms deviations calculated for all atoms excluding the first two and the last two basepairs.

A further confirmation of the stability of the simulated structure comes from the analysis of the average number of the interstrand Watson-Crick hydrogen bonds (black bars) and of the whole hydrogen bonds (gray bars) plotted in Fig. 3. The number of genuine Watson-Crick hydrogen bonds is strictly close to two, and to three for the A-T and C-G basepairs, respectively. A lower number is only observed for the last two basepairs at one end of the helix, which leads to a subtle structural deformation at the tail of the DNA helix. Comparison of our data with the ones previously extracted from an MD simulation for a DNA duplex belonging to the trp operator fragment (Bonvin et al., 1998) indicates that our structure displays a higher ability to maintain the inter-strand hydrogen bonds. This may be due either to a higher intrinsic stability of our DNA segment or to the avoidance of truncation effect in the treatment of the electrostatic interactions because of the use in our simulation of the particle mesh Ewald method (Essman et al., 1995), which has been shown to be well suited for the treatment of electrostatics in the case of strongly charged systems like DNA.

View larger version (58K): [in this window] [in a new window]   FIGURE 3   Average number of inter-strand basepairing hydrogen bonds (gray bars) as a function of the DNA sequence calculated from the last 2.5-ns trajectory. The fractions of the Watson-Crick hydrogen bond are indicated in black. The criterion for the occurrence of a hydrogen bond is defined in Computational Methods. The box indicates central basepairs involved in BamHI recognition.

Hydration analysis

The trajectory has been analyzed to work out information on the degree of hydration around the various groups forming the DNA sequence. Before performing the analysis, the presence of time-dependent artifacts (Sagui and Darden, 1999; Harvey et al., 1998) due to the use of a relatively long nonbond update time, i.e., 20 fs, was checked, evaluating the water dynamics throughout the trajectory. The water kinetic energy was fairly constant during the overall simulation, its variation being <5%, allowing us to exclude the presence of any artifact. In the analysis we have evaluated the average and maximum residence time and the coordination number of water molecules around specific sites. These data are reported for the backbone and the major and minor groove, respectively, and are discussed in comparison with other structural an dynamical values obtained in previous work.

Backbone

The coordination number and the average and maximum residence times of water molecules around the phosphate and ribose oxygen atoms evaluated from the last 2.5 ns of simulation are reported in Table 1. The data indicate that of the four phosphate oxygens, the two partially charged phosphate oxygens, O1P and O2P, are highly hydrated, their average water coordination number ranging between 3.5 and 3.8. However, the coordination number around the ester oxygens, O3' and O5', oscillates around 1. An interesting observation is that the oxygen phosphate and the ester oxygen coordination number does not depend on the presence of specific bases because these values are fairly constant all along the DNA sequence. Such a finding is in agreement with the analysis carried out by Schneider et al. (1998) on 59 crystallographic structures of DNA. In that case it is reported that each charged phosphate oxygen (O1P and O2P) has three water-ordered molecules in its first coordination sphere. The lower coordination number observed by Schneider is likely due to the fact that in his work the selection has been made on water molecules hydrogen-bonded to the oxygens, so that the considered constraints (water-O1P and -O2P distance 2.8 Å and hydrogen bond angle >125°) are more strict than in our case, where we have considered all the water molecules having an oxygen phosphate-oxygen-water distance 3.5 Å. However, in both analyses it is found that water has a high affinity for the phosphate oxygens and that ester oxygens are hydrated only marginally.

The two types of oxygen atoms behave differently also as far as the water residence time is concerned. Phosphate oxygen atoms, in fact, have average and maximum residence times that are about twice those displayed by the ester oxygen atoms (see Table 1). This could be attributed to the difference in partial charge, which is greater for the phosphate (0.776) than for the ester (0.509) oxygen atoms, resulting in a long-lasting water interaction for the first class of oxygens. A comparable result is also found in an MD simulation of the trp operator (Bonvin et al., 1998) indicating that phosphate group hydration is independent of DNA sequence.

Major and minor groove

The average and maximum water residence times and the coordination number of water molecules around the various sites of the bases belonging to the major and minor groove are reported in Tables 2 and 3, respectively. In the tables the average solvent accessibility surface of the various sites evaluated from the MD trajectory is also reported in comparison with the solvent accessibility surface worked out from the DNA structure as observed in the x-ray diffraction of the BamHI-DNA complex (Newman et al., 1995), once the DNA sequence and the protein are separated. The two values are similar for all the major groove sites and some deviation is only observed in the minor groove around the Gua-5 base. This result is confirmed by the average values of the minor groove width taken over 500-3000 ps, reported in Fig. 4 A, in comparison to the values extracted from the DNA structure of the BamHI-DNA complex. The largest difference is observed in correspondence of the Gua-5 base. The average twist, obtained from the MD simulation, is found to have values lower than those observed in the x-ray structure (Fig. 4 B), suggesting that the amber force field underestimates this quantity, as reported in other MD simulations (Cheatham and Kollman, 1997b; Young et al., 1997). However, these results, together with the trend of the structural parameters of Fig. 2 confirm, that MD is sampling conformations quite close to the x-ray structure.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Minor groove average width (A) and average twist (B) for the six central basepairs calculated from the last 2.5-ns trajectory () compared with the values evaluated from the crystal structure (×). Error bars are single standard deviation about the mean. The two parameters were averaged from values extracted every 50 ps.

In the major groove each pyrimidine has one hydration site; in particular, at least one water molecule is present during all the trajectory in the first coordination sphere of the N4 of cytosine and O4 tymine, respectively (Table 2). The average water residence time is slightly shorter for water molecules residing in the proximity of a cytosine than of a tymine, however in both cases water is fast-exchanging with the bulk solvent, indicating that these sites, although highly hydrated, do not trap water molecules for a long time. In the case of purines, they both display two hydration sites in correspondence to the N7 and N6 atoms of the adenine and N7 and O6 atoms of the guanine base, respectively. The average coordination number is around 1 for both the N7 and N6 atoms of the adenine bases, while it is higher for the guanine atoms. The water molecules detected in the major groove are not shared between atoms belonging to different strands or between atoms belonging to the same base except for the Ade-6 and Ade-18, where one water molecule is shared for 2% and 20% of the simulation trajectory, respectively, between the N7 and N6 atoms. This effect slightly decreases the coordination number of the N7 and N6 atoms of the adenine bases and makes the C-G pairs more hydrated than the A-T basepairs, as may be observed by summing the average coordination number of the CG and the AT reported in Table 2. A higher hydration of the C-G basepairs when compared to the AT pairs was found in a recent long MD simulation (10 ns) of the DNA fragment d(C5T5)(A5G5) (Feig and Pettitt, 1999). This result is also in agreement with the hydration values evaluated by molar volume and compressibility measurement as a function of base composition (Chalikian et al., 1994) and by free energy calculations (Elcock and McCammon, 1995). The water residence times of the purines are also characterized by relatively short values (Table 2).

In the minor groove (Table 3) the average water coordination number for each base is close to one except for the guanine base, which displays a non-shared water molecule around both the N2 and N3 atoms, again confirming the higher degree of hydration of the C-G than the A-T basepairs. The average water residence time is short and comparable to that observed in the major groove. Some differences between major and minor grooves are observed at the level of the maximum water residence time, which is longer in the minor than in the major groove. The number of water molecules with maximum residence time greater than 80 ps, reported in Fig. 5 for the six central basepairs, is higher in the minor than in the major groove, indicating that water is trapped for a longer time in the minor groove. As a matter of fact, the minor groove shows the presence of localized water molecules. In fact, the 5'-A₆T₇C₈-3' strand shares water molecules with the opposite 5'-A₁₈T₁₉C₂₀-3' bases for a significant percentage of time. In Fig. 6 a scheme of the minor groove, where base atoms sharing water molecules are connected by a line, has been reported. There are pairs of atoms, such as C₈ O₂ and A₁₈ N3, T₇ O₂ and T₁₉ O₂, A₆ N3 and C₂₀ O₂, which share a water molecule for 24%, 44%, and 28% of the total simulation time, respectively. This implies that a water molecule has a high probability of being found in a position bridging two atoms belonging to the two strands despite the fact that the average water residence time is short. These positions are local minima for water molecules and create a spine of hydration, as already observed in the minor groove of other DNA sequences (Duan et al., 1997; Feig and Pettitt, 1999). Fig. 5 also shows that there are water molecules in the minor groove that may be shared for a different percentage of time by various atoms, sometimes even by four sites, although for a relatively short percentage of the simulation time, i.e., 5%. Such shared water molecules are not observed in the major groove, indicating that water localization occurs only in the minor groove. However, despite this localization, the average water residence times are relatively short also in the minor groove, although slightly longer than in the major groove (see Tables 2 and 3), confirming that the static and dynamic properties of water are not necessarily correlated (Levitt and Park, 1993). In fact, a localized water molecule, i.e., a water molecule revealed by x-ray diffraction, has a favored average position where the water is almost always present, although likely fast-exchanging with bulk water. This behavior is better clarified in Fig. 7, where the time history of water molecules in the first coordination sphere of T₇ O₂ and T₁₉ O₂ are represented. Each water molecule is represented by a different color and the fast changes of colors indicate that the two sites are almost continuously hydrated, but by many different water molecules. The figure also represents as a function of time the water molecules that are shared by the above-mentioned atoms and it confirms that although water molecules have a high probability of bridging the two sites during the trajectory, each water molecule resides in this position for a relatively short time.

View larger version (48K): [in this window] [in a new window]   FIGURE 5   Number of water molecules residing on each base having a maximum residence time >80 ps in the major (gray bars) and minor (black bars) grooves, respectively.

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Connectivity of the hydration sites in the six central basepairs of the DNA duplex recognized by the BamHI endonuclease. Atoms sharing a water molecule, having a distance <3.5 Å from each atom, are connected by lines. The reported value represents the percentage of time over the whole simulation time in which water molecules are shared by the connected atoms. Lines are solid, dashed, and dotted for water molecules shared by two, three, and four sites, respectively. The values represent the percentage of simulation time in which water molecules are shared between the various sites.

View larger version (32K): [in this window] [in a new window]   FIGURE 7   Hydration time history of the O2 T7 and O2 T19 atoms. Each color represents a different water molecule at a distance <3.5 Å from the O2 T7 atom (top line) and from the O2 T19 atom (bottom line). The middle line represents the fraction of time in which water molecules are shared between these two atoms.

The average residence times of water molecules within 3.5 Å of nonexchangeable DNA protons were also calculated using a time resolution of 10 ps, i.e., counting water molecules leaving and entering the 3.5 Å radius within 10 ps as being continuously within the distance cutoff. The values shown in Table 4 are three to four times lower than those reported in a previous MD study carried out on a different DNA sequence using the same resolution time (Bonvin et al., 1998). Such a difference can, at least in part, be due to the use in our study of the TIP3P water model, which has a diffusion constant higher (Van der Spoel et al., 1998) than the SPC model used in the Bonvin work. Decreasing the time resolution to 0.1 ps reduces the evaluated residence time by one order of magnitude (data not shown) indicating that resolution time is the parameter that must be calibrated when comparing experimental and calculated water residence times.

Upon considering the importance of the choice of the resolution time the calculated residence time values are comparable to the ones experimentally measured through NMR on various B-DNA sequences (Sunnerhagen et al., 1998; Phan et al., 1999). In Phan's work it is also found that the residence times of the major groove were shorter than those in the minor groove as found in our simulation, suggesting that this trend is a general feature of the grooves that likely depends on their geometrical properties and not on the specific DNA sequence.

The x-ray diffraction of the BamHI-DNA complex have indicated that Gua-5 interacts with the protein through water-mediated hydrogen bonds (Newman et al., 1995). We don't find any long-lasting water molecule in this or any other position; however, we find that Gua-5, together with Gua-17, is the most hydrated base. A similar result has been obtained in the MD simulation of the trp operator (Bonvin et al., 1998) where no specifically long residence time values were observed in any site corresponding to the water crystallographic sites.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The structure and hydration of the DNA sequence interacting with BamHI has been investigated by means of a 3.0-ns molecular dynamics simulation in water. The use of the particle mesh Ewald method allowed an accurate treatment of the electrostatics and the generation of a stable MD trajectory. The solvated structure in the MD simulation reaches a stable configuration close to the configuration observed through x-ray diffraction in its complex with the BamHI endonuclease. This is in line with the evidence that, in the DNA-BamHI complex, the DNA remains quite rigid, while the protein undergoes conformational change to fit with the DNA sequence (Newman et al., 1995)

Hydration analysis indicates a preferential hydration for phosphate oxygen atoms than for ester oxygen atoms, as already observed in other static and dynamical analysis (Falk et al., 1962; De Oliviera Neto, 1986; Schneider et al., 1998; Bonvin et al., 1998). Water residence times are relatively short for both the backbone and the major and minor groove sites, being slightly shorter in the major than in the minor groove, as also experimentally reported in other DNA duplexes (Phan et al., 1999). A spine of hydration is found in the minor groove along the central CTA sequence confirming that it is not a unique feature of the AATT sequence (Duan et al., 1997; Feig and Pettitt, 1999).

No special long water residence time values were observed for sites where crystallographic waters have been detected, as also found in the case of the trp operator. These results indicate that the sites where water molecules are observed by x-ray diffraction correspond to sites where water molecules have a high probability of being found, but these molecules are usually fast-exchanging with the bulk solvent.