Molecular Dynamics Simulation of Methylene Blue−Guanine Complex

---

J. Phys. Chem. B 2000, 104, 1073-1077

1073

Molecular Dynamics Simulation of Methylene Blue-Guanine Complex in Water: The Role of Solvent in Stacking Mironel Enescu,*,‡ Bernard Levy,‡ and Victor Gheorghe§ National Institute of Laser, Plasma and Radiation Physics, P.O. Box MG-36, 76900-Magurele, Romania, Laboratoire de Physico-Chimie des Rayonnements (UMR no.8610), UniVersite Paris Sud, bat. 350, 91405 Orsay Cedex, France, and Biophysics Department, UniVersity of Bucharest, 76900-Magurele, Romania ReceiVed: July 20, 1999; In Final Form: October 14, 1999

The conformation of the methylene blue (MB)-guanine (GUA) complex in water was investigated by molecular dynamics simulation. It was found that the T-shaped vacuo stable conformation becomes unstable in water and turns into a well-defined stacked conformation. The free energy of the complex formation was determined by the thermodynamic integration method and found equal to -7.2 ( 0.2 kcal/mol. The variation of the hydration energy upon complexation was found to be positive (about 14 kcal/mol), being compensated by the solute-solute interaction energy (about -21 kcal/mol). The main contribution to this solute-solute interaction energy arises from the van der Waals forces (about -14 kcal/mol).

1. Introduction The role of water in the stacking of the planar aromatic molecules is largely debated.1-8 Two processes are usually considered: (i) the organization of the complex by the solvent (the hydrophilic domains of the complexed molecules are oriented toward the solvent while the hydrophobic ones are pushed one against the other2,9) and (ii) the reorganization of the solvent structure around the two molecules upon complexation.3,10,11 From an experimental point of view it is difficult to establish the relative importance in stacking of the solute-solute dispersion forces and of the solvent reorganization energy. Moreover, in many cases it is considered that the stacking conformation results from specific molecular interactions such as charge-transfer4,12 or electrostatic interactions13-15 rather than by hydrophobic interaction. The recent development of the molecular simulation method provides a new and interesting tool in the present field. The simulation method has been successfully applied in the study of ionic hydration,16-18 dimerization of small hydrophobic molecules,19 and complexation of large organic molecules.7,20 The molecular simulation studies of benzene dimer in water21-23 have demonstrated the importance of the solute-solute forces in the dimer conformation. It was also suggested that the stacking is not a general property of aromatic molecular aggregates. On the other hand, the heterocyclic molecules seems always to form a stacked complex in water,3,15,24,25 their molecular association being important in many areas of biological interest. For instance, it is involved in the photosensitizing activity of planar organic dyes causing nucleic acid damage.26-28 Among these dyes methylene blue (MB) represents a model molecule, and its interaction with the nucleic acid bases is of particular interest. In a previously reported molecular dynamics (MD) simulation, the conformation of the MB-guanine (GUA) complex in vacuo was found to be T-shaped, nonstacked.29 Spectroscopic †

National Institute of Laser, Plasma and Radiation Physics. Universite Paris Sud. § University of Bucharest. ‡

data29 as well as the flat form of the two molecules strongly suggest that the conformation of the MB-GUA complex in water is a stacked one. In the present paper we try to answer the question whether the MD simulation predicts a stacked structure when the solvent is included. We also calculate the free energy of complexation and analyze the contributions arising from different processes: solute-solute interaction and solvent reorganization. A detailed description of this particular molecular association process is obtained. It reveals the importance of the solvent in the conformation of the complex. 2. Method In a MD simulation, a trajectory in the phase space of the molecular system is generated by integrating Newton’s equations of motion. The set of molecular configurations obtained along the trajectory has to form a representative statistical ensemble which may be used to evaluate the statistical distribution of different variables describing the system. In particular, the most probable conformation of a molecular complex may be identified. In the present work the MD simulations have been performed with the CHARMM program30 using CHARMM 22 force field parameters on a HP 725/80 graphic station. The molecular parameters of GUA (partial atomic charges, geometry, nonbonded and bonded energy parameters) were taken from ref 31. The partial atomic charges and the geometry of MB were derived by an ab initio molecular orbital calculation already reported,32 while the van der Waals parameters were transferred from the CHARMM 22 atomic parameters file. Since the intramolecular bonded energy is not important in the present calculation, the corresponding bonding parameters of MB were adjusted in order to preserve the molecular geometry during MD simulation. The SHAKE routine was used to constrain all covalent bonds including hydrogen atoms. The electrostatic and van der Waals potentials were progressively shifted (by SHIFT and VSHIFT options) between 10 and 12 Å in order to obtain a smooth transition to zero. GUA and MB were included in a water sphere of 15 Å radius. Water molecules were represented by the TIP 3 model.33 A static force field acting on the water molecules was added in order to

10.1021/jp992486f CCC: $19.00 © 2000 American Chemical Society Published on Web 01/11/2000

1074 J. Phys. Chem. B, Vol. 104, No. 5, 2000

Figure 1. The position and orientation of the reference frames attached to MB and GUA, respectively.

simulate the action of the solvent molecules outside the sphere.34 This boundary condition preserves the density into the sphere during simulation. The water sphere was equilibrated at 300 K, its density being 0.99 g/cm3. The time step chosen for the integration of motion equations was 1 fs. Prior the simulation, the system was heated to 300 K in 6 ps then equilibrated for 60 ps at the same temperature. The simulation was carried out at constant temperature by coupling the system to a Berendsen thermostat with a relaxation time of 0.1 ps. The boundary potential used in the present simulation preserves the mean water density during simulation and indirectly the mean pressure. The statistical ensemble of configurations used to determine the molecular complex conformation was obtained by taking snapshots of the trajectory. The spatial conformation of the complex can be described by means of six parameters. A reference frame was attached to each molecule as shown in Figure 1. The first three parameters are the spherical coordinates of the origin O′ of the frame attached to the guanine molecule with respect to the reference frame attached to MB. The remaining three parameters describe the relative orientation of the axis of the two frames. The variation interval of the configuration parameters was reduced in order to take into account the molecular symmetry. The free energy for MB-GUA association was derived by the thermodynamic integration method16 performing two artificial deletions of GUA in the solvent sphere (as suggested in ref 19): first for GUA alone, then for GUA complexed to MB. The energy terms involving GUA were multiplied by a λ2 factor, where λ is the reaction parameter varying from 1 to 0. The λ value was diminished by a constant amount at each dynamics step (1 fs). The thermodynamic integration was performed during 200 ps. The deletion of guanine would cause a density variation in the system not greater than 2%, as expected from the variation of the total atom number in the system. In fact, this variation should be significantly reduced by the action of the boundary potential. The remaining effects of density variation in the two deletion processes were supposed to compensate each other and were neglected in the computation of the free energy. 3. Results The Conformation of the MB-GUA Complex in Water. As was already suggested,20 the simulation of a molecular association process in water should start from molecules separately diffusing in the solvent sphere and continue until the molecules meet one another and form a stable configuration. This approach requires extensive simulation without evidence that the resulting conformation corresponds to an absolute free energy minimum. In the present work we are mainly interested by the solvent effect on the vacuo stable conformation of the

Enescu et al. complex: in particular we want to answer the question whether a T-shaped, nonstacked conformation (Figure 2a) is still stable or if it is modified to a face-to-face conformation (Figure 2b). Hence, we started the MD simulation from a MB-GUA conformation close to the vacuo stable conformation. The simulation was extended on a time interval of 500 ps. A stable conformation was reached after about 100 ps. The derived statistical distributions of the six conformation parameters are presented in Figure 3. The distribution of the distances between the “centers” of the two molecules (Figure 3a) shows a maximum at 3.5 Å as expected for face-to-face complexes of planar molecules. The width at half-maximum is only about 0.5 A, suggesting a well-defined structure. However, the angular coordinates of the guanine “center” are rather dispersed (Figures 3b and 3c). Concerning the relative orientation of the two molecules in the complex, the most interesting result is the statistical distribution of the angle γ between the molecular planes (Figure 3f). A sharp maximum centered at 10° indicates a defined faceto-face structure. The present result is in clear contrast with the corresponding conformation obtained for vacuo interaction32 showing a flat distribution ranging from 20° to 60°. We note that the results reported in ref 32 were derived using the same set of force field parameters as in the present paper. The remaining two parameters, R and β, describing the relative orientation of the guanine axis with respect to the MB axis, exhibit also narrow distributions (Figure 3 d, e), meaning that the relative orientation of the two molecules is stable, which is not the case for the complex in vacuo. We note also that the transition between the starting T-shaped conformation and the face-to-face conformation is discontinuous as illustrated in Figure 4. This is a characteristic feature of processes controlled by solvent fluctuations. To test the reproducibility of these results we performed additional simulations starting from different T-shaped initial conformations with a few different values of the R, β, γ parameters. In all of these cases no stable T-shaped structure was found, the only stable structures being of face-to-face type. In some cases the molecules separated one from another without complex formation. Free Energy Computation. In a first step, the solvent properties have been tested by performing a 200 ps simulation at 300 K on the 15 Å radius sphere of pure solvent containing 464 molecules. The derived potential energy, -4160 kcal/mol, corresponds to a heat of vaporization of 9.0 kcal/mol, in very good agreement with the value reported for TIP model (8.9 kcal/ mol).33 As a further test, the free energy of hydration of Na+ was computed by deleting the solvated ion in a 200 ps thermodynamic integration process. A hydration energy of -105.4 kcal/mol was found, close to the experimental value -100.2 kcal/mol.16 In the thermodynamic integration method the variation of the free energy of the system when the reaction parameter λ is modified from λA to λB is calculated as

FBA )

dλ ∫AB 〈∂H ∂λ 〉λ

(1)

where the brackets represent an ensemble average obtained using the Hamiltonian H(λ). The difference FBA was calculated in two cases: (i) deletion of GUA from the MB-GUA complex in water and (ii) deletion of GUA alone in pure water. The difference between the free energies obtained in cases (i) and (ii) gives the complexation free energy. It was found equal to

Methylene Blue-Guanine Complex in Water

J. Phys. Chem. B, Vol. 104, No. 5, 2000 1075

Figure 2. Possible conformation of MB-GUA complex: T-shaped, nonstacked (a) and face-to-face (b).

Figure 3. The geometry of the MB-GUA complex in water. (i)The distribution of the spherical coordinates of the GUA center O′ with respect to the reference frame attached to MB centered in O: (a) the OO′ distance, (b) the spherical angle θ, and (c) the spherical angle φ. The angle θ is formed by the OO′ vector and the O z-axis, while φ is the angle between the projection of the OO′ vector on the xOy plane and the O x-axis. (ii)The frequency distribution describing the relative orientation of the two molecules: (d) R, the angle between the axes Ox and O′x′, (e) β, the angle between the axes Oy and O′y′, and (f) γ, the dihedral angle between the molecular planes.

-7.2 kcal/mol. The integration was started with λ ) 1 (corresponding to GUA present in the system) and stopped at λ ) 0 (corresponding to GUA completely removed). The values

of the integral FBA corresponding to integration intervals equal to 0.05 taken on the reaction trajectory are given in Figure 6 for both cases.

1076 J. Phys. Chem. B, Vol. 104, No. 5, 2000

Enescu et al.

Figure 4. Variation of the dihedral angle γ during the first 100 ps of the molecular dynamics simulation illustrating the transition between the T-shaped conformation and face-to-face conformation.

Figure 6. Dependence on λ (reaction parameter) of the average interaction energy between MB and GUA, respectively, between GUA and solvent molecules.

By substituting this equation into eq 1 one finds for small variation of λ

∆FBA ) 〈EG〉λ + 〈UG-S〉λ + 〈UG-MB〉λ 2λA∆λ

Figure 5. Variation of the free energy of the system upon the deletion of GUA in the solvent sphere (open circles) and deletion of GUA complexed to MB in the solvent sphere (solid circles). The differences ∆F have been computed for variation intervals of the reaction parameter λ of 0.05.

As the reaction parameter is varied the system undergoes successive transformations. The conformation of the complex begins to be affected at λ ) 0.6. At λ ) 0.5, the two molecules are separated and the solvent fills the cavity previously occupied by GUA. The guanine is slowly diffusing to the solvent boundary. It leaves the solvent at λ ) 0.3. During this process the solvent orientation around the two molecules also is affected. To illustrate the evolution of the system it is interesting to evaluate the average solute-solute and solute-solvent interaction energy during the simulation. When the derivative of the Hamiltonian is calculated (eq 1) only the explicit dependence on λ should be considered. In our case, the Hamiltonian part explicitly dependent on λ is

U(λ) ) λ2(EG + UG-S + UG-MB)

(2)

with EG being the internal energy of GUA, UG-S the interaction potential between GUA and solvent, and UG-MB the interaction potential between GUA and MB. After derivation with respect to λ and statistical averaging one obtains

〈∂U∂λ 〉 ) 2λ(〈E 〉 + 〈U λ

G λ

G-S〉λ

+ 〈UG-MB〉λ)

(3)

(4)

Equation 4 allows the evaluation of the sum in the right member based on the values represented in Figure 5. The contribution of the internal energy of GUA was then eliminated by performing a separate deletion of GUA in vacuo. The variation of the two remaining terms with respect to λ representing the change of the solvent orientation around GUA and the dissociation of the complex is represented in Figure 6. Usually, the estimation of the error in the thermodynamic integration method can be obtained by performing a reverse simulation, starting from the product: the difference between the forward and reverse free energy variations gives the computational accuracy. In our case the reverse simulation is not possible: indeed, when the two molecules are separated the rebinding of the complex cannot be realized in 200 ps. Instead, the precision of the result was tested by repeating the guanine deletion in the complex in two other separated simulations. Finally a accuracy of (0.2 kcal/mol was estimated. We note also that the variation of the free energy (Figure 5) as well as that of the intermolecular potential energy (Figure 6) are smooth, showing that the disintegration of the complex and the diffusion of GUA are slow enough with respect to the variation of λ. Otherwise, sudden variations of the free energy should be detected resulting from a non continuous molecular reorientation. Analysis of Different Contributions to the Free Energy of Complexation. The free energy of complexation in water may by viewed as the balance between the solute-solute interaction energy and the variation of the hydration energy upon complexation. In the case of the MB-GUA complex we have

∆F ) ES(C) - ES(MB) - ES(G) + EI(MB-G)

(5)

where ES(C), ES(G), and ES(MB) are, respectively, the hydration energy of the complex, GUA, and MB and EI(MB-G) is the interaction energy between MB and GUA. Given the fact that the conformation of the complex is relatively stable, the contribution of the interaction energy may be evaluated without

Methylene Blue-Guanine Complex in Water extensive molecular simulation. Hence we have randomly extracted 10 configurations from the dynamic trajectory previously obtained and calculated the direct MB-GUA interaction energy. The mean value was -21.2 kcal with a dispersion of (2.5 kcal/mol. Then, using the final value of ∆F (-7.2 kcal/ mol) and eq 5 one obtains the variation of the hydration energy upon complexation (sum of the three first terms of the righthand side of eq 5) as being positive and equal to 14 kcal/mol. It is to be noted that the dominant contribution to the solutesolute interaction energy arises from the van der Waals interaction (about -14 kcal/mol). This is a different situation as compared to a vacuo complex where the dominant contribution is electrostatic. 4. Discussion According to the present results, the following model for the methylene blue-guanine complexation in water emerges. The solvent destabilizes the vacuo stable T-shaped conformation which is based on the electrostatic interaction between the partial atomic charges of the two molecules. Under favorable solvent fluctuation a single cavity holding the two molecules in a faceto-face configuration is formed. There are three factors involved in the MB-GUA association in water. The first factor is solvent reorganization resulting in a negative contribution to the free energy of hydration. This factor has been considered by some authors as the most important one in the hydrophobic hydration.10 However, the solvent reorganization energy computed for the transition of benzene dimer from the T-shaped to the stacked conformation was reported to be not greater than -2 kcal/mol.22 A small contribution of the solvent reorganization to the free energy of methane dimer formation in water has also been reported.8 In fact, a recent study35 based on Monte Carlo simulations pointed out that the present contribution should not be very important at room temperature because the coordination number of a solvent molecule around the cavity slowly decreases when the solvent-accessible surface increases, but with little energy loss. In addition, it results from a recent model of hydrophobic hydration that the contribution of the solvent reorganization to the solute-solute association should be always small because the entropic and the enthalpic contributions to the free energy compensate each other.36 The second factor is the variation of the solvent-solute interaction. In the case of benzene dimer, the variation of the solvent-solute interaction energy upon transition from the T-shaped to the stacked structure was found to be about 2 kcal/ mol.22 In the present case the variation of the solvent-solute interaction energy upon molecular association should be important given the charge polarity of MB and GUA molecules. Indeed, the present factor and the previous one give here a net positive contribution to the free energy of the process of about 14 kcal/mol which may be viewed as a desolvation energy penalty. The third factor is the solute-solute interaction energy. A net contribution of about -21 kcal/mol compensates the desolvation energy. The most important contribution to the solute-solute interaction arises in the present case from the van der Waals forces. The complexation is thus possible because the van der Waals interaction between the solvent and the hydrophobic regions of the solute molecules is considerably smaller than the corresponding solute-solute interaction. The role of the solute-solute interaction in the conformation of aromatic and heterocyclic dimers has already been demonstrated. The T-shaped conformation of the benzene dimer is favored

J. Phys. Chem. B, Vol. 104, No. 5, 2000 1077 by the quadrupole interaction of monomers22 while the stacked conformation of adenine dimer in water has been attributed to the electrostatic interaction between the adenine molecules.15 The van der Waals interaction has been proved to be important in the case of toluene dimer,23 molecular association between porphyrins and quinones,2 and the molecular structure of (2,2)-bis(indol-1-y1-methyl) accetate.7 The fact that van der Waals interaction is the main driving force in the methylene blue-guanine association in water is a consequence of the solvent action: the polar parts of the two molecules are prevented from being in the most favorable relative position by the solvent molecules, while the apolar ones are promoted into close contact. It is by this mechanism that the solvent contributes indirectly to the stacked conformation. This mechanism could be also active in the case of many other heterocyclic complexes of biological importance. References and Notes (1) Dewey, T. G.; Wilson, P. S.; Turner, D. H. J. Am. Chem. Soc. 1978, 100, 4550. (2) Kano, K.; Hayakawa, T.; Hashimoto, S. Bull. Chem. Jpn. 1991, 64, 778. (3) Mukerjee, P.; Ghosh, A. K. J. Am. Chem. Soc. 1970, 92, 6419. (4) Dahl, T. Acta Chem. Scand. 1994, 48, 95. (5) Valdes-Aguilera, O.; Neckers, D. C. Acc. Chem. Res. 1989, 22, 171. (6) Lee, B.; Graziano, G. J. Am. Chem. Soc. 1996, 118, 5163. (7) Yuan, P. P.; Miller, J. L.; Kollman, P. A. J. Am. Chem. Soc. 1999, 121, 1717. (8) Lu¨deman, S.; Absher, R.; Schreiber, H.; Steinhauser, O. J. Am. Chem. Soc. 1997, 119, 4206. (9) Spencer, W.; Sutter, J. R. J. Phys. Chem. 1979, 83, 1573. (10) Sinanoglu, O. Int. J. Quantum Chem. 1981, 18, 381. (11) Marglit, R.; Rotenberg, M. Biochem. J. 1984, 219, 445. (12) Smith, M. L.; McHale, J. L. J. Phys. Chem. 1985, 89, 4002. (13) Mecozzi, S.; West, A. P., Jr.; Dougherty, D. A. J. Am. Chem. Soc. 1996, 118, 2307. (14) Dougherty, D. A. Science 1996, 271, 163. (15) Newcomb, L.; Haque, T. S.; Gellman, S. H. J. Am. Chem. Soc. 1994, 116, 4993. (16) Straatsma, T. P.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 5876. (17) Obst, S.; Bradaczek, H. J. Phys. Chem. 1996, 100, 15677. (18) Garcia-Tarres, L.; Guardia, E. J. Phys. Chem. B 1998, 102, 7448. (19) Jorgensen, W. L.; Buckner, J. K.; Boudon S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, 3742. (20) Mordasini Denti, T. Z., van Gusteren, W. F.; Diederich, F. J. Am. Chem. Soc. 1996, 118, 6044. (21) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768. (22) Linse, P. J. Am. Chem. Soc. 1992, 114, 4366. (23) Chipot, C.; Jaffe, R.; Maigret, B.; Pearlman, D. A.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 11217. (24) Mutai, K.; Gruber, B. A.; Leonard, N. J. J. Am. Chem. Soc. 1975, 97, 4095. (25) Dunn, D. A.; Lin, V. H.; Kochevar, I. E. Photochem. Photobiol. 1991, 47, 53. (26) Kelly, J. M.; van der Putten, W. J. M.; McConnell, D. J. Photochem. Photobiol. 1995, 62, 55. (27) Enescu, M.; Lindqvist, L. Photochem. Photobiol. 1995, 62, 55. (28) Tuite, E.; Norden, B. J. Am. Chem. Soc. 1994, 116, 7548. (29) Enescu, M.; Gheorghe, V. J. Mol. Struct. (THEOCHEM) 1998, 423, 213. (30) Brooks, R. B.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (31) MacKerell, A. D., Jr.; Wiorkiewicz-Kuczera, J.; Karplus, M. J. Am. Chem. Soc. 1995, 117, 11946. (32) Levy, B.; Enescu, M. J. Mol. Struct. (THEOCHEM) 1998, 432, 235. (33) Jorgensen, W. L. J. Am. Chem. Soc. 1981, 103, 335. (34) Brunger, A.; Brooks, C. L., III; Karplus, M. Chem. Phys. Lett. 1984, 105, 495. (35) Ikeguchi, M.; Shimazu, S.; Nakamura, S.; Shimazu, K. J. Phys. Chem. B, 1998, 102, 5891. (36) Grunwald, E.; Steel, C. J. Am. Chem. Soc. 1995, 117, 5687.