J . Phys. Chem. 1990, 94, 4660-4665
4660
would occur because of direct interaction with the tip and substrate (Q and 7R in our model). If the consequent broadening is not too large, it might be possible to obtain chemical information on a nanometer scale. Acknowledgment. We are grateful to George Wolf, Roland Hanson, and Tom Moore for the loan of equipment, and to John
Page, Tom Thundat, Larry Nagahara, Rick Oden, and Lin Oliver for useful discussions and advice. Roger Wartell encouraged this work by pointing out the contrast variations in STM images of DNA. Pui S. Ho and Leigh Clarke shared their work with us prior to its publication. We received financial support from the NSF (BBS8615653) the O N R (N00014-87-K-0478) and the office of the Vice President for Research at ASU.
Molecular Dynamics Simulation of a Deoxydinucleoside-Drug Intercalation Complex: dCpG/Proflavin S. Swaminathan, D. L. Beveridge,* Chemistry Department, Wesleyan University, Middletown. Connecticut 06457
and H. M. Berman Institute f o r Cancer Research, Fox Chase Cancer Center, 7701 Burholme Avenue, Philadelphia, Pennsylvania 19111 (Received: August 11, 1989; In Final Form: January 2, 1990)
A molecular dynamics (MD) study of the dCpG/proflavin crystal hydrate based on the GROMOS force field and simulation
protocol is reported. The objectives of the study are (a) to determine the extent to which the theoretical calculations account for the dynamical motion of the complex and the organization of the water network and (b) to investigate the nature of the intermolecular hydrogen bonding network and determine the extent of interconnections on a static and dynamical basis. The results show the GROMOS energy function to have a well-defined energy minimum with a root-mean-square deviation of 0.19 A for the complex and 0.38 A for the water. The hydrogen bond network in the energy-minimized water structure is not as fully connected as was indicated by the 3.2-A crystallographic contacts. MD calculations were performed on a system of two unit cells (eight asymmetric units) for 25 ps at a temperature of 300 K. The dynamical behavior of the complex was reasonable, even to the point of reproducing mixed sugar puckers. The water structure was examined in terms of a hydrogen bond density representation. Here, in a dynamical sense and symmetry averaged over asymmetric units, the water network was found to be fully hydrogen bonded and consistent with 3.2 8, contact representation. The contact representation is thus a good indication of the time-averaged hydrogen bond network.
Introduction
The crystal structure of 2:2 complex of the dinucleoside dCpG and the simple intercalator proflavin reported in a series of papers by Neidle, Beman, and Shieh’-*provided the first example of a nucleic acid system with highly structured hydration. The dCpG/proflavin crystal has subsequently proved to be of considerable interest from the point of view of both structural studies and theoretical calculations. The system provided considerable detail on the conformational flexibility of the nucleic acid at the intercalation site and dynamics of the intercalator. However, the nature of the water network proposed by crystallographers has been questioned3 with regard to whether the virtal bonds drawn at “contact” distances of 3.2 A accurately represent the underlying intermolecular hydrogen bond network or not. The dCpG/proflavin crystal, as a well-determined, intermediate-sized system, provides a prototypical case for evaluating the performance of molecular simulations on nucleic acids and a basis for investigating the sensitivity of results to the choice of force field, potential function truncations, and other characteristics and assumptions in the calculations as well as assumed system size. We describe herein further theoretical studies on the dCpG/ proflavin crystal hydrate based on both energy minimization (EM) calculations and molecular dynamics (MD) simulations and investigate further the issues raised above. A comparison of the 0 K EM structure with the instantaneous (snapshot) and diffu( I ) Shieh, H. S.;Berman, H. M.; Dabrow, M.; Neidle, S. Nucleic Acids Res. 1980, 8, 85. (2) Neidle, S . ; Berman, H. M.; Shieh, S. Nature 1980, 288, 129. (3) Savage, H. In Water Science Reuiews; Franks, F., Ed.; Cambridge University Press: Cambridge, U.K., 1986; Vol. 1.
0022-3654/90/2094-4660$02.50/0
sionally averaged structures, termed by Eisenberg and Kauzmann4 the “I” and “D” structures, respectively, permits the investigation of the dynamical as well as static aspects of structure in this system. This differentiation has proved to be crucial in understanding the nature of the hydrogen bond network in the structured water.
Background The dCpG/proflavin complex is shown in Figure 1. The crystal consists of stacked intercalated duplexes. The general organization of water in the crystal is shown in Figure 2. The observed hydration network is shown in Figures 3 and 4. The major groove cavity is filled with water molecules arranged in pentagonal arrays, and the minor groove features a flat polygonal disk of waters. In the arrays there are five edge-linked pentagons, four of which consist solely of waters while the fifth consists of four water molecules and a phosphate oxygen. These pentagons are anchored to the complex via the heteroatoms on the bases and the proflavines. Previous theoretical studies on nucleic acid crystal hydrates based on molecular simulation have been carried out on dCpG/proflavin and related systems. Mezei et aLs carried out a Monte Carlo (MC) simulation on the water in the crystallographic unit cell, with the nucleic acid adduct atoms fixed in their X-ray positions. Here 108 water molecules were included, a number chosen on the basis of independent density measurements (4) Eisenberg, D.; Kanzmann, W. Structure and Properties of Water; Oxford University Press: New York, 1969. ( 5 ) Mezei, M.: Beveridge, D. L.; Berman, H. M.; Goodfellow, j. m.;Finney, J . L.; Neidle, S . J . Biomol. Srruct. Dyn. 1983, 1 , 287.
0 1990 American Chemical Society
Molecular Dynamics Simulation of dCpG/Proflavin
Figure 1. Molecular structure of the 2:2 dCpG/proflavin adduct in the crystalline solid. Top: view parallel to the plane of the proflavin. Bottom: view perpendicular to the plane of the proflavin.
The Journal of Physical Chemistry, Vol. 94, No. 1 I , I990 4661
Figure 3. Oxygen atoms of the ordered water network, viewed along the crystallographic b axis, for the dCpG/proflavin crystal hydrate. A portion of the network extending over three unit cells is shown. Virtual bonds are drawn between oxygen atoms within a hydrogen bond contact distance of 3.3 A.
9 \
P
\
I
/
Y
Figure 2. General organization of water in the dCpG/proflavin crystal. Water molecules are represented by 3.2-A spheres centered on the observed water oxygen positions.
of the crystaL6 The result after some 2000K configurations of sampling was partially but not completely successful in accounting for the observed ordered water positions, with some 64%of the crystallographic sites found to be within 0.6 A of a calculated “generic” solvation site.’ Clementi and co-workerss-I0 reported simultaneously an MC simulation on the water of the dCpG/proflavin crystal hydrate, treating a system of 12 unit cells so that the long-range ionic potentials due to the phosphate anions and proflavin cation would not have to be truncated unrealistically. They also varied the number of waters in the crystal and, on the basis of energetics and positional analysis9 as well as solvation sites,I0 predicted that the optimum number of waters in the crystal should be 122-1 32 per unit cell. With this number, extremely good agreement on the subset of 100 waters closest to the crystallographic positions ( 6 ) Neidle, S.; Goodfellow, J. M. Private communication. (7) Mezei, M.; Beveridge, D. L. J . Compur. Chem. 1984, 5, 523. (8) Kim, K. S.; Corongiu, G.; Clementi, E. J. Biomol. Srrucr Dyn. 1983, I , 263. (9) Kim, K. S.; Clementi, E. J . Am. Chem. SOC.1985, 107, 227. (IO) Kim, K . S.; Clementi, E. J . Phys. Chem. 1985,89, 3655.
Figure 4. Structure of an asymmetric unit of dCpG/proflavin hydrated crystal.
was obtained. However, the density arrived at by this approach is considerably larger than that expected by the crystallographers and estimated from experimental density measurements. The agreement with observed water positions is bound to improve at higher densities, since the space is more tightly packed with waters and some will obviously end up close to the observed positions. The authors discuss possible sources of errors in their calculation in terms of potential functions and the neglect of cooperative effects. Additional Monte Carlo simulations on nucleic acid crystal hydrates have been reported by Elliot and Goodfellow.”J2 The small crystal hydrate systems relevant to nuclei acid solvation that have been studied by MD simulation to date are ~ - ~cyclodextrins ~ had earlier the a- and ~ - c y c l ~ d e x t r i n s . ~The (11) Elliot, R. J.; Goodfellow, J. M. J . Theor. Biol. 1987, 127, 403. (12) Elliot, R. J.; Goodfellow, J. M J . Theor. Biol. 1987, 128, 121. ( I 3) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. F.Eur. Biophys.
J .. . 1987. 15., 197. .- . , . .. . . (14) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W.
J . 1987, 15, 21 I .
F.Eur. Biophys.
4662
The Journal of Physical Chemistry, Vol. 94, No. 11, 1990
Swaminathan et al.
been the subject of neutron diffraction studies, and the waters of hydration were well determined with respect to both position and orientation. Koehler et al. carried out MD calculations using GROMOS on the a-cyclodextrin hexahydrate crystal13 and the 0-cyclodextrin dodecahydrate cry~ta1.l~ They obtained generally excellent agreement with experiment on crystals and also compared the crystal results with simulations on aqueous hydration. The calculations generally support the idea of flip-flop hydrogen bond dynamicsi5 and also provide evidence for occurrence of novel three-center hydrogen bond~’~J’ in the crystalline solid and aqueous solution.
Calculations The calculations in this study are based on empirical energy functions for the intra- and intermolecular interactions in the dCpG/proflavin crystal hydrate and involve locating the local minimum with a geometry closest to the crystal structure on the energy surface and then probing the energetic and structural characteristics of this minimum using molecular dynamics. Comparisons of the calculated results with the crystal structure permit an objective assessment of the extent to which the calculations account for the observed results. Detailed analysis of the calculations leads to information on the dynamics of the intercalation site and the nature of the hydrogen bond network in the hydration complex. The calculations described herein are based on the GROMOS system of energy functions and the GROMOS86 molecular simulation routines.i6 All parameters for the dCpG duplex are in the standard GROMOS prescription. The parameters for the proflavin were obtained from the standard GROMOS atom types with the N and H charges in the NH+ unit assigned to be +0.6 and +0.4, respectively. The object of theoretical study was a system of two completely intact unit cells consisting of 8 intact 2:2 dCpG/proflavin adducts (1056 atoms) and 216 SPC water molecules (648 atoms), a total of 1704 atoms and 5 109 degrees of freedom. The system for simulation was constructed by doubling the crystallographic elementary cell in the c direction, thus providing a simulation cell of dimensions a = 32.991 A, b = 21.995 A, and c’= 2c or 27.018 8, where c = 13.509 A. This permits a minimum truncation of potentials at b/2 or 10.998 8,. A switching function was implemented to feather the potentials smoothly to zero over the range from 9.5 to 10.5 8, and applied on a group by group rather than atom by atom basis, with PO4- and NH+ being the only groups that carry a net charge. Bond length constraints were applied to all covalent bonds to hydrogen atoms by using SHAKE,I9 and the MD time step was 0.001 ps. The number of water molecules was assumed to be 27 per asymmetric unit, obtained by converting each of the two fractional occupancies to unity by expansion across the symmetry boundary and adding an additional water molecule at a vacancy located by a crude Monte Carlo calculation on the distribution of unoccupied space in the crystal structure. As in the previous s t ~ d y this ,~ number of waters is chosen to be consistent with crystallographers’ estimates and density measurements on the system6 Energy Minimizations The EM studies were carried out to locate the local minimum on the energy surface closest to the crystal geometry. The intrinsic characteristics of the EM structures are also of interest. EM structures were obtained by three different protocols, hereafter referred to as EMI, EM2, and EM3. All minimizations were initiated from the crystal structure coordinates. The index of (IS) Koehler. J. E. H.; Saenger, W.; van Gunsteren, W. F. Eur. Biophys. J . 1988, 16, 153. ( 1 6 ) Koehler. J. E. H.; Saenger, W.; van Gunsteren, W. F. J . Mol. Bioi. 1989, 203, 241. ( 1 7 ) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. F. J . Biomol. Struct. Dyn. 1988, 6, 181. (18) van Gunsteren, W.F.; Berendsen, H. J. C. GROMOS: Groningen Molecular Simulation Program, University of Groningen, 1986. (19) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J . Comput. Phys. 1977, 23, 327.
Figure 5. Comparison of three energy-minimized structures, EMI, EM2, and EM3 (thin lines) with the molecular structure of dCpG/proflavin (thick lines) in the crystal hydrate.
quality in the EM’S was measured as the temperature achieved in one 0.001 -ps displacement in time, departing from the minimized structure. This index is essentially an inverse mass-weighted root-mean-square (rms) force, converted into the equivalent temperature. EM1 involved simply 500 steps of free minimization using the conjugate gradient (CG) method.*O The final temperature of the complete system was 70 K. The rms deviation of the EM structure from the crystal geometry was 0.35 8, for the dCpG/ proflavin complex and 0.64 8, for the water. The molecular dynamics calculations described in the following section were initiated from EM 1. In parallel, some refinements of the EM studies were explored. EM2 involved an initial optimization of the orientation of water molecules using 1000 passes of a Monte Carlo Metropolis procedure,2i with the water oxygens restrained to their crystallographic coordinates with a harmonic restraint function with a force constant of 25 kcal/(mol.8,). This was followed by 50 steps of CG free minimization, which achieved a system temperature of 7.5 K with an rms deviation of 0.21 8, for the crystal structure complex and 0.56 8, for the water. EM3 involved a protocol for introducing asymmetric unit symmetry into the minimization. Here an initial structure for the system based on the crystal geometry was subjected to 100 passes of Metropolis Monte Carlo simulation. The energy of the eight asymmetric units was examined, and the one corresponding to the lowest energy was selected. The most energetically favorable asymmetric unit structure was then used to form a new complete system of eight in a simulation cell. The process was iterated 10 times over. The minimization process was completed by 50 steps of the CG calculation. By use of this protocol, a final system temperature of 1 K was achieved. The minimized structure obtained from EM3 deviated by 0.19 8, from the crystal structure of the complex and 0.38 8, from that of the water. A comparison of the minimized structures for a typical asymmetric unit from EMl, EM2, and EM3 with the crystal structure is shown in Figure 5. Small but distinct deviations are seen for EMI, but EM2 and EM3 are virtually coincident with the crystal structure in this diagram. Thus the GROMOS energy surface is seen to have a well-defined local minimum in the vicinity of the crystal structure and behaves in a quite satisfactory manner in this regard. The calculated water structure for the crystal obtained from the EM3 study is shown in Figure 6 , displayed in the convention established by crystallography in Figure 3 and averaged with respect to asymmetric units. Here virtual bonds are drawn be(20) Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling, W. J. Numerical Recipes: The Art of Scientific Computing Cambridge University
Press: Cambridge, 1986. (21) Metropolis, N.;Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087.
Molecular Dynamics Simulation of dCpG/Proflavin
The Journal of Physical Chemistry, Vol. 94, No. I I, 1990 4663
Figure 8. Dynamical range of motion for dCpG/proflavin atoms in the MD simulation. The spheres on each atom are an isotropic analysis of the motion, with radii given by (l/2)(A3)’’z. The factor of 1/2 is introduced for clarity in the visualization of the structure.
/ \
I+c 1
Figure 6. Water network in the dCpG/proflavin crystal hydrate as determined by the energy minimization EM13. Virtual bonds are drawn only between atoms at distance C3.5 A with a pair interaction energy < I .5 kcal mol.
0
2
C
I
0
2
c 1 0
2
C
I
0
2
Figure 7. Comparison of three dynamical I structures at t = 0, 12.5, and 25 ps with the observed molecular structure (thicker lines) of dCpG/
proflavin complex in the crystal hydrate. tween water oxygens 53.50 A of each other, and with a pair interaction energy 5-1.5 kcal/mol. Thus only “hydrogen-bonded” contacts are portrayed in Figure 6 . The polygon disk and the pentagon network are recognizable but not fully connected in Figure 6, the “0 K” energy-minimized structure. The apparent discrepancy at the apex of the polygon disk is due to our particular resolution of the fractional occupancy problem.
Molecular Dynamics The molecular dynamics calculations were carried out by using standard GROMOS protocols and a time step of 0.001 ps. The pair interaction neighbor list was based on a cutoff of 10.9 A and updated at 0.05-ps intervals in the calculation. The calculation assumed as starting point the EM1 minimized form of the crystal and was first heated to 300 K and equilibrated. The heating phase involved random assignment of velocities from a Gaussian distribution corresponding to T = 50 K, followed by an isothermal MD step of 0.5 ps duration. The velocities were reassigned to provide and increment of 50 K in temperature and the MD step was repeated. This heating cycle was continued for six iterations, at which point the temperature of the system was 300 K. The heating was followed by 7 ps of Gaussian equilibrium and 15 ps of free equilibration. The production phase of the calculations was commenced at this point and extended for a duration of 25 ps. The vital signs of the calculation were all found to be stable. The dynamical structure of the dCpG/proflavin adduct in the crystal hydrate determined from the MD simulation is depicted in Figures 7 and 8. A comparison of the observed molecular
Figure 9. Conformational analysis of the left-hand strand of asymmetric units I , 3, 5 , and 7 in the system.
structure of the adduct from the crystal hydrated with the I structures (snapshots) from the MD at t = 0, 12.5, and 25 ps for an arbitrarily selected symmetric unit is shown in Figure 7 and provides a good idea of the dynamic range of the molecular motion. The preferential disorder observed as mixed sugar puckering is correlated with atomic motions in the calculations. In Figure 8 the dynamical behavior of the adduct averaged over all asymmetric units is represented in terms of spheres proportioned to (1/2)(Ar2)’I2for each atom, with the 1/2 factor introduced arbitrarily for better visibility. A fully detailed dynamical conformational analysis of the simulations has been carried out by use of the method of “dials and windows”22based on the “CU~V~S’’ analysis of the s t r ~ c t u r e . ~ - ~ ~ The dial systems for each of the strands, left- and right-handed, for each of the IUPAC conformational parameters (a, 0, 7,6, c, [, x) and pseudorotation phase angle 6 of each of the eight dCpG/proflavin complexes in the two-unit-cell system are shown in Figures 9-12. The conformational dials are defined on a conventional Klyne-Prelog system with the origin at the top (N) and the conformational angle increasing in a clockwise direction. (22) Ravishanker, G.;Swaminathan, S.; Beveridge, D. L.; Lavery, R.; Sklenar, H. J . Biomol. Srr. Dyn.1989, 6,669. (23) Lavery, R.; Sklenar, H. J . Biomol. Srrucr. Dyn.1988, 6,63. (24) Lavery, R.; Sklenar, H. J . Biomol. Srrucr. Dyn.1989, 6,655.
4664
s 0
Q
e
C
a
P
u
x
G
4
E
