J. Phys. Chem. 1996, 100, 2535-2545
2535
Quantum-Classical Molecular Dynamics Simulations of Proton Transfer Processes in Molecular Complexes and in Enzymes P. Bała,†,§ P. Grochowski,‡,§ B. Lesyng,*,‡,§ and J. A. McCammon*,⊥ Institute of Physics, N. Copernicus UniVersity, Grudzia¸ dzka 5/7, 87-100 Torun´ , Poland; Department of Biophysics, Warsaw UniVersity, Zwirki i Wigury 93, 02-089 Warsaw, Poland; Interdisciplinary Centre for Mathematical and Computational Modelling, Warsaw UniVersity, Banacha 2, 02-097 Warsaw, Poland; and Department of Chemistry and Biochemistry, and Department of Pharmacology, UniVersity of California at San Diego, La Jolla, California 92093-0365 ReceiVed: September 11, 1995; In Final Form: NoVember 22, 1995X
A quantum-classical molecular dynamics model (QCMD) designed for simulations of proton or electron transfer processes in molecular systems is described and applied to several model problems. The primary goal of this work is the elucidation of enzymatic reactions. For example, using the QCMD model, the dynamics of key protons in an enzyme’s active site might be described by the time-dependent Schroedinger equation while the dynamics of the remaining atoms are described using MD. The coupling between the quantum proton(s) and the classical atoms is accomplished Via extended Hellmann-Feynman forces, as well as the time dependence of the potential energy function in the Schroedinger equation. The potential energy function is either parametrized prior to the simulations or can be computed using a parametrized valence bond (VB) method (QCMD/VB model). The QCMD method was used to simulate proton transfer in a proton bound ammonia-ammonia dimer as well as to simulate dissociation of a Xe-HI complex in its electronic excited state. The simulation results are compared with data obtained using a quantum-classical time-dependent self-consistent field method (Q/C TDSCF) and with results of fully quantum-dynamical simulations. Finally QCMD/VB simulations of a hydrolytic process catalyzed by phospholipase A2, including quantum-dynamical dissociation of a water molecule in the active site, are reported. To the best of our knowledge, these are the first simulations that explicitly use the time-dependent Schroedinger equation to describe enzyme catalytic activity.
1. Introduction During the past few years, there has been increased interest in the development of quantum-dynamical methods and their application to chemical and physical problems at the molecular level. These methods are based on a variety of approximations. Some of them have a heuristic background, and others are derived from variational principles1,2 or are exact analytical or numerical solutions of the time-dependent Schroedinger equation.3-7 In most practical problems, the models contain both classical and quantum degrees of freedom. The separation of the phase space into quantum and classical domains will be dictated by the problem under study as well as by the required accuracy of the solutions. Electronic degrees of freedom are usually treated as quantum variables, whereas positions and momenta of heavy atoms are treated as classical ones. In timeindependent theories, the Born-Oppenheimer approximation is a common route to the separation of quantum and classical degrees of freedom; see ref 8. In general, the problem is more complicated. The separation process reduces the global dynamical problem into a system of coupled, lower-dimensional equations of motion. Effective approaches for constructing potential energy surfaces for small molecular systems and incorporating quantum effects in classical trajectory simulations are presented in ref 9. In the reduction process, one can distinguish two phases. The purpose of the first one is to determine an effective Hamiltonian * To whom correspondence should be addressed. † Copernicus University. ‡ Warsaw University, Zwirki i Wigury. § Warsaw University, Banacha. ⊥ University of California at San Diego. X Abstract published in AdVance ACS Abstracts, January 15, 1996.
0022-3654/96/20100-2535$12.00/0
for the system under the study. During the second one, one separates the quantum and classical degrees of freedom by splitting the resulting dynamical problem into the classical and quantum equations of motion. Note that treating a subset of dynamical variables as a quantum one, one modifies once more the potential energy surface for the classical motion of the atomic nuclei. This can be made more specific by considering a proton transfer process. The original, unreduced molecular Hamiltonian contains electronic, protonic and heavy-nuclei degrees of freedom. Suppose that we want to describe the quantum proton dynamics in the environment of heavier, classical nuclei. In the simplest approximation, one assumes that the electronic charge distribution adjusts itself instantaneously to the current positions of all other degrees of freedom, and in this way generates an effective potential energy function for the protonic and nuclear degrees of freedom. One gets this potential by integration of the original Hamiltonian over the electronic coordinates, invoking the time-independent Born-Oppenheimer (BO) approximation, or another adiabatic approximation, see, e.g., refs 8, 10, and 11. In this way, one gets an effective Hamiltonian for the atomic nuclei, including protons. The numerically-obtained potential energy function from the BO approximation can be fit to analytical functions, which allows for the explicit definition of the effective Hamiltonian; see, e.g., refs 12 and 13 for practical solutions of the fitting procedures. It should be stressed that the quality of the potential energy surface is critical for the quality of further studies. For example, relatively small errors in the curvature of the fit can produce significant differences in the dynamics on this surface.14,15 In the second phase, either one can run the dynamics treating all © 1996 American Chemical Society
2536 J. Phys. Chem., Vol. 100, No. 7, 1996 atomic degrees of freedom classically by integration of the Newtonian equations of motion (see e.g. refs 16-18) or one can select the proton coordinates as the quantum-dynamical degrees of freedom. In the last case, one solves the timedependent Schroedinger equation for the protonic degrees of freedom, and the classical equations of motion for the heavy nuclei. Enzyme reactions have been studied using, among others, classical molecular dynamics (MD) simulations with empirical potential energy functions (EP). This class of the methods will be denoted as MD/EP. Using such methods, fragments of renzymatic eaction paths,19 aqueous polarization effects,20,21 and hydrogen bond dynamics in enzymes22 were studied. However, conventional MD methods with molecular mechanics force fields are limited; in particular, they cannot be applied for the description of chemical or biochemical reactions in which bondbreaking or bond-forming processes occur. To overcome this problem quantum-mechanical (QM) methods have been used for calculations of the potential energy functions and their gradients (forces). Below we briefly mention some of these methods. Time-independent semiempirical (SEM), ab initio (ABI), and density functional (DFT) methods have been applied to generate effective potential energy functions, and the forces (gradients of the potential functions) have been used in MD simulations. Models which belong to this class will be denoted as MD/QM (more specifically, MD/SEM, MD/DFT, or MD/ABI, respectively). In particular, molecular dynamics combined with the semiempirical methods AM1 and MNDO (MD/AM1 and MD/ MNDO, respectively) were used for simulations of polyatomic molecules22-24 and chemical reactions in vacuo25 and on surfaces.26,27 Examples of ab initio molecular dynamics calculations MD/ABI, called often first principles molecular dynamics, include studies of the dynamics of a silicon cluster28 and vibrations of C60.29 The MD/DFT methods are widely used in solid state physics and various molecular systems.30-35 During the past several years, a first principles variational model developed by Car and Parinello36 has become quite popular for molecular dynamics simulations. In this approach a fictitious kinetic energy term is used to obtain a system of coupled equations for the wave function and the nuclei, which permits an adiabatic propagation of the electronic state in parallel with the nuclear motions.37,38 This is a more efficient procedure than trying to find at each geometry a new optimal quantum state with the lowest energy value. The most important limitations of this method are related to the lack of quantum excitations in the system, which leads to an underestimation of tunneling phenomena. Recently a rigorous QD model was implemented in the field the time-dependent density functional theory (TDDF).39,40 This model contains a physical kinetic energy of the electrons (quantum particles). Furthermore, because of the rigorous relation to the Schroedinger equation for the full system, the TDDF equations conserve the total momentum and the total energy quite well. One of the effective quantum and quantum-classical methods is the time-dependent self-consistent field model (TDSCF).41 In this model, the coupling between the quantum and classical particles goes through an average potential energy surface. TDSCF was applied at a semiempirical level for simulations of the dynamics of small molecules.23 The theory of the electronnuclei dynamics (END)42 can be considered as an extension of the TDSCF methods in such sense that electrons and nuclei are allowed to interact without any restrictions. The electron-nuclei dynamics takes place in a generalized phase space. The END
Bała et al. method has been applied to study ion-atom scattering problem43 as well as intramolecular proton transfer processes.42 However, since the ab initio and density functional methods are computationally expensive, their applications in quantum and quantumclassical simulations are limited to small molecular systems. Recently a density matrix evolution (DME) method44 was developed to simulate the dynamics of quantum systems embedded in a classical environment. The DME model was applied to study the inelastic collisions of a classical particle with a quantum harmonic oscillator,45 and proton tunneling in an intramolecular hydrogen bond in aqueous solution.46,47 We have developed a quantum-classical molecular dynamics model (QCMD), which allows one to simulate dynamical processes in quantum-classical systems.48-51 In particular, this model can be used to simulate the motion of quantum particles, protons or electrons, in molecular systems and in the solid state. The ultimate targets of our studies are enzymatic reactions. In the QCMD model, the proton dynamics in enzyme active sites are described by the time-dependent Schroedinger equation, and the dynamics of the remaining atoms is performed using classical molecular dynamics. Coupling between the quantum proton(s) and the classical atoms is accomplished Via extended Hellmann-Feynman forces50 as well as the time dependence of the potential energy function in the Schroedinger equation. We have also evaluated corrections to the wave function resulting from fast motions of the classical particles.50 The potential energy function was parametrized prior to the simulations.12,13,52 At the testing stage, the QCMD method was applied to simulate quantum proton hopping in model proton-bound dimers48-51 and to simulate a tautomeric rearrangement in malonaldehyde.48 Recently we formulated and parametrized an approximate valence bond (VB) method53 for proton transfer and hydroxyanion nucleophilic reactions in enzyme catalytic processes. The parametrization of the VB method is based on density functional as well as conventional ab initio calculations, and calibrated with respect to experimental data in the gas phase. Influence of an electrostatic field generated at the classical protein environment on the active site is mediated Via electrostatic fields. In this model, denoted as QCMD/VB, the VB method is used for the calculations of the potential energy surface of some of the classical nuclei as well as the quantum protons. The method was applied to enzyme reactions.53,54 Hydrolytic enzymes like phospholipases, serine proteases, and acetylcholinesterases are the primary targets of our studies. Figure 1 presents a classification scheme of the quantum, classical and quantum-classical models, used mostly in simulations of proton transfer processes in molecular systems. 2. Quantum-Classical Molecular Dynamics Models Used in This Study A. QCMD Model. In the QCMD model, selected degrees of freedom are treated quantum dynamically, and the others classically or stochastically. In general, the model consists of a system of coupled equations which govern the dynamics of the quantum and classical subsystems (Figure 2). The dynamics of the quantum subsystem is described by the time-dependent Schroedinger equation, while the rest of the system is described by the Newtonian equations of motion. The coupling between the quantum and classical domains is described by the time-dependent potential function V ) V(x,{RR(t)}) in the Schroedinger equation, where V is a function of the instantaneous atomic positions RR, and by extended Hellmann-Feynman forces50 modifying the classical forces in the Newtonian equations of motion (eq 2). These coupling terms are important for the energy transfer between the quantum
PT Processes in Molecular Complexes and Enzymes
J. Phys. Chem., Vol. 100, No. 7, 1996 2537 Corrections to the wave function which account for fast motions of the classical nuclei have also been evaluated,49,50 but in practical applications for proton transfer processes, including the results presented in this study, can be neglected. The forces between classical atoms, FR are computed using standard potential energy functions. The forces between quantum particles and classical atoms located sufficiently far from the quantum subsystem are calculated in the classical limit; in practice, when the distance to a remote atom is much larger than the hydrogen bond length, the classical forces can be used instead of the rigorous extended Hellmann-Feynman forces. B. TDSCF Model. In practical applications, one of the most successful quantum dynamical models is the time-dependent self-consistent field (TDSCF) method developed by Alimi at al.41,55 For any single coordinate, the effective potential in the TDSCF equations is time-dependent and represents an average field due to the motion of all other degrees of freedom. In our implementation of the model, some of the coordinates are described quantum mechanically while other modes are treated classically. This scheme is usually referred as the quantum/ classical TDSCF (Q/C TDSCF).41 With this approximation the time-dependent wave function of the system is assumed to be separable:
Figure 1. Schematic classification of the quantum, classical, and quantum-classical models.
|ψ(x,R,t)〉 ≈ |Ξ(R,t)〉|ψ(x,t)〉
(3)
Substituting this form into the time-dependent Schroedinger equation, multiplying by Ξ* or ψ*, and integrating over R and x, respectively, one gets the quantum dynamical TDSCF equations:
(
)
∂ψ(x,t) p2 ip ) - ∆x + Veff(x,t) ψ(x,t) ∂t 2m ∂Ξ(R,t) ip
∂t
(
)
p2 ∆R + Veff ) -∑ c (R,t) Ξ(R,t) R 2MR
(4)
(5)
where
Figure 2. Scheme of the QCMD model. R and P denote positions and momenta of the classical particles. ψ is the wave function of the quantum particle(s). VC and VQ are the potential energy functions of the classical and quantum regions, respectively (see Figure 1). FC is the classical force and FQ is the force resulting from the extended Hellmann-Feynman theorem.
and classical domains as well as for energy dissipation processes which influence the time evolution of the whole system.49,50 Assuming that the positions of the classical nuclei change slowly in time, one can treat them as parameters in the quantum Hamiltonian, and in the simplest representation the QCMD equations are the following:
(
)
∂ψ(x,t;{RR}) p ip ) - ∆x + V(x;{RR}) ψ(x,t;{RR}) ∂t 2m 2
(〈|
MRR¨ R ) FR - ∑ bk Φk 2
k
|Φ 〉 + E ∂R )
∂V(x;{RR}) ∂RR
k
(6)
Veff c (R,t) ) 〈ψ(x,t)|V(x,R)|ψ(x,t)〉x
(7)
The time dependence of the effective potentials allows for the energy flow between the subsystems. The mixed quantum/ classical scheme is obtained by taking the classical limit of the equation for the R coordinates,
MRR¨ R ) -∂Veff c (R,t)/∂RR
(8)
The mean field acting on the quantum particle becomes
Veff(x,t) )
1
N
∑V(n)(x,R) Nn)1
(9)
(1)
where n labels different initial values sampling the different initial conditions for the classical trajectories R(n).
(2)
3. Potential Energy Surface
∂bk2
k
Veff(x,t) ) 〈Ξ(R,t)|V(x,R)|Ξ(R,t)〉R
R
where Φk and Ek are instantaneous eigenfunctions and eigenvalues of the time-independent Hamiltonian H ) -(p2/2m)∆x + V(x;{RR}). bk are the expansion coefficients of the wave function in the basis of eigenfunctions of H: ψ(x,t;{RR}) ) ∑k)1bke-(i/p)EktΦk(x).
In the case of proton transfer, the potential energy function V(x;{RR}) is a Born-Oppenheimer potential energy surface which is obtained numerically by solving the time-independent electronic Schroedinger equation, varying the proton position and computing the energy of the system. For this purpose conventional ab initio or semiempirical molecular orbital methods as well as valence-bond orbital or density functional
2538 J. Phys. Chem., Vol. 100, No. 7, 1996
Bała et al.
Figure 3. Coordinate system for the colinear Xe-HI cluster.
TABLE 1: Parameters Used in the Numerical Integration of the QCMD Equations of Motion parameter
Xe-HI
[H3N-H-NH3]+
phospholipase A2
no. of grid points ∆x, Å time step ∆t, fs no. of Chebychev polynomials
256 × 2048 0.028 × 0.0065 0.1 33
128 × 128 0.029 × 0.0375 0.05 51
256 0.029 0.25 60
methods can be applied. As mentioned, in order to study the proton transfer processes in molecular complexes and in the active site of phospholipase A2 we have developed a valencebond (VB) orbital method.53 Several resonance, valence- and ionic-type structures are included in the calculations, and the VB matrix elements are parametrized to reproduce ab initio MP2/6-31G* and/or density functional calculations for these sites. The influence of remote classical atoms on the active site is effected by explicit inclusion of electrostatic interactions with remote atomic charges in the VB Hamiltonian’s diagonal elements. The discretized energy surface changes its shape in time according to the actual atomic positions and actual value of the electrostatic field generated in the active site region by the remote atoms. Changes of the shape of the potential energy surface play an important role in the proton quantum dynamics and lead to the proton transfer. The proper treatment of these coupling terms is important for the accurate description of the energy transfer between the quantum and classical domains, as well as the energy dissipation processes which influence the time evolution of the whole system.49,51,54
Figure 4. Simulation results for the Xe-HI cluster. 〈r〉 is the mean interatomic HI separation (a); R is the distance between the Xe atom and the center of mass of the HI molecule (b).
4. Simulation Results and Discussion The QCMD and Q/C TDSCF methods were applied to study photodissociation of the Xe-HI cluster as well as to study proton transfer in a model proton-bound ammonia-ammonia dimer and in the active site of phospholipase A2. In this active site, due to a strong electrostatic field, dissociation of a water molecule takes place; this is the first step in the catalyzed reaction. In the case of the first two systems, the quantumclassical simulation results will be compared with full QD simulations. A. Photodissociation of Xe-HI. The triatomic linear system Xe-HI (see Figure 3) is described in Jacobi coordinates r and R, where R is the distance between the Xe atom and the center of mass of the HI molecule and r is the interatomic HI separation. In the QCMD model, the hydrogen atom is moving along the quantum coordinate r while the heavy particles follow the classical equations of motion. The potential energy surface for the excited state of the linear Xe-HI cluster is approximated as a sum of pairwise interactions between the atoms,41 a potential developed based on experimental spectroscopic data. The initial R and r distances are taken as equilibrium values in the ground electronic, state of the Xe-HI complex with the proton wavefunction of a Gaussian shape. In the Q/C TDSCF model the R coordinate is represented by five initial conditions chosen at Req, Req ( ∆R, and Req ( 2∆R, with ∆R ) 0.1 Å.41 The grid parameters and the time step are listed in Table 1. Additionally, full QD simulations were performed where both
Figure 5. Conservation of the total energy for the Xe-HI cluster using different models.
R and r coordinates were treated quantum mechanically. The initial wave function was assumed to be a product of two Gaussians centered at the ground state equilibrium values of R and r, respectively. The mean values of the H-I distance obtained in both the QCMD and Q/C TDSCF schemes tend to have smaller amplitudes of oscillation than in the full quantum case (Figure 4a). In all simulations the Xe-(HI) distance is practically the same and the trajectories overlap (Figure 4b). The total energy of the system is also reasonably conserved (Figure 5), although the Q/C TDSCF simulations show small oscillations. The total energy in the exact quantum dynamical simulations is higher by 0.07 kcal/mol due to the zero-point energy of the R mode. A slightly higher energy in the Q/C TDSCF model, compared to the QCMD scheme, comes from the average of the potential energy over the classical trajectories. Another quantity of interest related directly to the wave function is the square of the autocorrelation function a(t) ) |〈ψ(x,0)|ψ(x,t)〉|2. As
PT Processes in Molecular Complexes and Enzymes
J. Phys. Chem., Vol. 100, No. 7, 1996 2539
Figure 6. Autocorrelation function obtained using different models.
Figure 7. Coordinate system for the proton-bound ammonia dimer.
presented in Figure 6, the QCMD and TDSCF results are very close and agree qualitatively with the exact QD calculations. B. The Proton-Bound Ammonia-Ammonia Dimer, [H3NH-NH3]+. In the proton-bound ammonia-ammonia dimer, the proton is described by the time-dependent Schroedinger equation while the ammonia molecules are represented as classical particles with effective masses equal to 17mp (Figure 7). The potential energy surface for this system was fitted to an analytical function13 using ab initio MP2/6-31G* calculations.12 At distances R > 2.5 Å the potential for the proton transfer is bistable with the barrier increasing with R, and at distances R < 2.5 Å the barrier disappears. The equilibrium distance is R ) 2.72 Å. The initial ammonia-ammonia distance was stretched by 0.13 Å from its equilibrium value, and an initial Gaussian wavepacket describing the proton state was localized in one of the wells, at its lowest energy level. 1 ps simulations were performed with time step ∆t ) 0.1 fs (cf. Table 1). The oscillations of the classical distance R are shown in Figure 8, and the mean proton position 〈x〉 in Figure 9. Due to the initial localization of the wave function in one well of the potential, the system is in a nonstationary state and the quantum calculations show periodic oscillations of the proton wave function between wells. The changes of the amplitude of the proton motion come from the fact that the potential energy surface cannot be represented as sum of terms which depend exclusively on one coordinate. The R mode and the proton density distribution are coupled. As a result, two very close tunneling frequencies are visible in the full QD model. The proton motion within the QCMD scheme is close to the exact one, although only one tunneling frequency is observed. In Q/C TDSCF, the proton trajectory oscillates with a higher frequency (twice the main tunneling frequency), which suggests that coupling between the quantum and classical subystems is overestimated in this case. The classical oscillations of the R mode in the QCMD and Q/C TDSCF models are compared with the full quantum results; see Figure 8. Due to the coupling with the x coordinate small oscillations of the mean value of R and its rms deviation are observed in the full QD calculations. In all cases the equilibrium distance is similar. The QCMD trajectory exhibits periodic oscillations of the ammonia-ammonia distance in the region bounded by the root mean square deviation of the quantum wave
Figure 8. Distance between the classical atoms in the proton-bound ammonia-ammonia dimer obtained within the QCMD and Q/C TDSCF models (solid line). For comparison, the mean position of the R mode obtained using the exact QD model is presented (dashed line). The root mean square deviation around the mean ammonia-ammonia distance obtained in QD is also plotted (dotted lines).
function. The constant amplitude of these oscillations comes from balanced energy partition between the R and x coordinates and reflects proper coupling between these modes. The Q/C TDSCF model shows the amplitude of the R oscillations increasing with time, which is the result of energy flow from the quantum mode to the classical one. The total energy of the system in all models is well conserved (Figure 10). The QCMD simulation shows the total energy to be smaller than the exact QD value by 0.62 kcal/mol, which is equal to the zero-point energy of the R mode. The Q/C TDSCF model gives a relatively high energy of the system (6.11 kcal/ mol compared to the exact value of 3.59 kcal/mol). This is the result of the averaging of the potential energy surface over the classical trajectories, which is strongly dependent on R. The average potential for the proton motion in the Q/C TDSCF model is steeper and higher total energy as well as higher proton tunneling frequency are obtained. The trajectories obtained within the QCMD and QD models are quite similar. One observes a difference in the amplitude of the R oscillationssthe QCMD trajectory oscillates in the region localized around the mean position of the quantum wave function with an amplitude close to the mean-square deviation of the wave function, while in the Q/C TDSCF model the amplitude tends to increase above the region defined by the exact quantum calculations. C. Phospholipase A2. The QCMD method was designed for the modeling of enzymatic reactions. Initial simulations were performed for phospholipase A2, an enzyme which consists of about 110 residues and which hydrolyzes phospholipids. 1. Phospholipase A2: Description of the Enzymatic System. Several X-ray structures of phospholipase A2 are known; they have a high level of homology. In our studies, crystal structures
2540 J. Phys. Chem., Vol. 100, No. 7, 1996
Bała et al.
Figure 11. Structure of phospholipase A2 including the substrate molecule.
Figure 12. A substrate molecule for phospholipase A2. The arrow points out the ester bond to be hydrolyzed. Figure 9. Mean proton position 〈x〉 in the proton-bound ammoniaammonia dimer obtained within the QCMD and Q/C TDSCF models. For comparison, the mean proton position obtained in the exact QD calculations is presented (thin line).
Figure 13. First two steps of the enzymatic reaction catalyzed by phospholipase A2. Figure 10. Total energy of the proton-bound ammonia dimer.
of both native56 and inhibited57,58 enzymes were used from a cobra venom. The secondary structure of the enzyme molecule consists of three connected R-helices. The active site is located between two parallel R-helices. The active site comprises a Ca2+ coordinated with polar functional groups of protein and water molecules. Histidine (His47), located 6 Å away from the Ca2+ ion, plays an important role in hydrolysis. The crystal structure shows that this Histidine 47 is bonded to an aspartic acid residue (Asp99) and to structural water with hydrogen bonds (Figure 11). This is similar to the situation in serine proteases or acetylcholinesterases, in which the role of the water molecule is played by a serine residue. The crystal structure of the enzyme molecule with an inhibitor was also determined.57-59 The inhibitor replaces all the water molecules in the active site. The substrate molecule (Figure 12) replaces all water molecules but one. This particular water molecule occupies the space between His47 and the calcium ion. Based on various crystallographic and biochemical data,
the following mechanism of hydrolysis was proposed. At the first stage of the enzymatic reaction, a proton is transferred from the water molecule to the histidine residue, leaving the OHgroup prepared for the nucleophilic attack at the ester bond of the substrate; see Figure 13 for a description of the first two steps of the enzymatic reaction. The intermediate product of the attack is the tetrahedral substrate-hydroxy anion adduct stabilized by the calcium ion.60 During the catalysis, the oxyanion, calcium ion and the attacking nucleophile are buried in the enzyme’s interior (Figure 11). Access of solvent to the active site is limited and the nearest water molecule has been found to be at least 5 Å from the attacking group. The calcium ion is essential both for the binding of substrate and for the catalysis.60 It lowers the activation barrier of the transition state and controls the substrate binding.57,58 In the uninhibited form of phospholipase A2 the calcium coordination sites are occupied by oxygens of two water molecules.56 Until now, theoretical tools were not able to provide any detailed information on the quantum dynamics of the catalytic processes. Difficulties in the modeling of full enzymatic reactions result
PT Processes in Molecular Complexes and Enzymes
J. Phys. Chem., Vol. 100, No. 7, 1996 2541
mostly from large changes of the electronic charge distribution, as well as quantum dynamical effects, including proton or electron tunneling61-63 during the catalytic process. 2. Potential Energy Surface for the Enzymatic Reaction: A Fast Valence Bond Method. A VB method was formulated and parametrized for the QCMD/VB simulations of the enzymatic reaction. The method was used to compute the potential energy surface for the atomic and proton motions in the enzymes active site. The QCMD/VB approach allows for direct microscopic simulation of the catalytic reaction and takes into account the interactions between atoms in the active site and atoms belonging to the classical protein environment, which is described using empirical potential (EP) functions. The method was parametrized based on DFT calculations performed for small molecules (water, imidazole, substrate) and a molecular complex modeling the enzymatic active site. The modeled enzyme was divided into two regions. The first one is called the classical region and contains the protein environment for the active site. The second one, called the quantum region, is equivalent to the active site and comprises the molecular fragments directly involved in the enzymatic reaction: the imidazole ring of His47, a water molecule, and a segment of the substrate (-CH2-CO-O-CH2-) with the ester bond to be hydrolyzed (see Figures 12 and 13 for the description of the active site and the substrate molecule). The total potential energy surface used in the QCMD/VB simulations is
VTOT ) VEP + EVB
(10)
The VEP potential is calculated using empirical potential energy functions (GROMOS force field in this case) and describes interactions within the classical region as well as the bonding and Lennard-Jones interactions between the classical and quantum regions. The VB component is the total energy resulting from the VB calculations for the quantum region, including also the electrostatic interactions with the classical region and the polarization effects inside the active site.53,54,64 It reproduces the ground-state energy of the isolated quantum region as well as the perturbation energy of this site in the presence of the external electric field generated by the classical protein environment. The potential energy surface for the proton motion is bistable, with a deep minimum located 0.98 Å away from the water oxygen. The second minimum, near the imidazole nitrogen, is 45 kcal/mol above the global one in the model of the isolated active site. The influence of the enzymatic environment reduces this difference to 20 kcal/mol.51,53 One should note that the barrier height is modified by changes of the position of the atoms both in the active site region and in the enzymatic environment.54 3. Force Fields. The classical interaction parameters used in the simulations are based on the GROMOS65 parametrization with implicit solvent. Several modifications of the electrostatic charges at the atoms in the neighbourhood of the Ca2+ cation were introduced as reported in our previous studies.50,51 Atomic charges on the water molecule, imidazole ring of His47, and the substrate directly involved in hydrolysis are obtained within the VB method and recalculated at the each time step. The charges on the substrate PO4- group and on the terminal NH3 group are fixed in the dynamics. Their values were obtained from density functional calculations using DMol.66 Because the VB method is used to describe the potential energy surface for the proton motion in the hydrogen bond formed by the water molecule and the imidazole ring of histidine, some of the bonded interactions are directly included in the quantum VB potential and must be removed from the
classical force field. In particular, the HOH bond angle and the HO stretching interactions are directly included into the VB potential energy surface. Similarly, all bonded interactions involving the most important substrate atoms are removed from the classical force field, and are included in the VB method. 4. Energy Minimization and Thermalization. The initial structure of the phospholipase A2-substrate complex was generated using the crystallographic data for the enzymeinhibitor complex.57,58 The inhibitor molecule was changed into the substrate by replacing part of the phosphonate group of the transition state analog by a carbon atom and one water molecule. The dynamics of the proton in the hydrogen bond formed by the His 47 and the water molecule was treated classically during the minimization and thermalization procedures. The potential energy surface for the atoms located in the active site are calculated at the each step using the VB model. The forces acting on the atoms were calculated as the potential energy gradients. Prior to the MD calculations, relaxation of the structure of the system was performed. 1000 steps of steepest-descent minimization were done followed by 10 0.1 ps NVT dynamics runs at a temperature of 5 K, with the time step of 0.5 fs. The final coordinates were chosen as a starting point for further energy minimization. 1000 one-step runs at a temperature of 5 K were performed, reassigning the velocities at each time step. In this way the excess kinetic energy was successfully removed from the system. The structure from the last step was used as a starting point in the standard thermalization procedure. A 20 ps dynamic simulation was performed with a time step of 0.5 fs. Every 0.2 ps, new velocities were generated from a Maxwellian distribution. Every 5 ps, the reference temperature was increased to approach a final value of 300 K. Then another 25 ps dynamics with velocity scaling characterized by a short relaxation time (0.1 ps) were performed at the temperature of 300 K. The last stage of the thermalization procedure consisted of 55 ps dynamics in the NVT ensemble with the relaxation time of 1 ps. The overall rms deviation from the crystallographic structure for the thermalized structure with the VB potential energy surface is smaller than that obtained exclusively with the GROMOS force field (2.272 vs 2.659 Å). The most important differences between the thermalized and crystallographic structures were found in both cases for flexible residues, such as lysine, at the surface of the enzyme. 5. Dynamics of the Enzyme-Substrate Complex. A 100 ps molecular dynamics simulation in the NVT ensemble was performed for the phospholipase A2-substrate complex using both the VB potential energy surface (MD/VB simulation) and the GROMOS force field (MD simulation) for the whole complex. The time step was 0.5 fs, and the temperature was set to 300 K with the relaxation time of the thermal bath equal to 1 ps. Proton transfer did not occur in either the MD/VB or the MD simulation, although a strong stretching of the imidazolewater hydrogen bond was observed in the first case. The mean N-O distance in the imidazole-water hydrogen bond, averaged over 100 ps, is shorter in comparison to the case when the VB potential surface for the proton dynamics is used (Table 2). The water molecule involved in the hydrogen bond with histidine is significantly closer (by 1.119 Å) to the substrate oxygen atom in the case of the MD/VB model in comparison to the GROMOS force field case. The length of the water OH bond involved in the hydrogen bond with histidine is larger by 0.021 Å in the case of the VB parametrization, whereas the length of the OH bond not involved in the hydrogen bond is very similar in both
2542 J. Phys. Chem., Vol. 100, No. 7, 1996
Bała et al.
TABLE 2: Comparison of the Averages from 100 ps Trajectories Obtained with the VB Potential Energy Surface (MD/VB), and with the GROMOS Force Field (MD) for the Histidine-Water Complex, and Averages from a 48 ps Trajectory of the Phospholipase A2-Substrate Complexa Im(ND1)-Ca2+ Im(NE2)-Asp(OD2) Im(ND1)-water(O) quantum bond (water) classical bond (water) water(O)-(C2P) water(HWC)-(O22) C2P-O22 (substrate) a
MD
MD/VB
MD70
5.8952 ( 0.0006 2.6061 ( 0.0021 3.1265 ( 0.0043 0.9805 ( 0.0005 0.9987 ( 0.0003 3.8069 ( 0.0088 3.3502 ( 0.0068 1.2367 ( 0.0004
5.3842 ( 0.0039 2.5895 ( 0.0033 2.7685 ( 0.0018 1.0015 ( 0.0003 0.9828 ( 0.0001 2.6876 ( 0.0001 3.2412 ( 0.0029 1.2306 ( 0.0003
5.6 3.1 3.9
All distances are in Å.
cases. In the MD/VB model, the active site is more compact and the calcium Ca2+ ion is about 0.5 Å closer to the histidine than when the GROMOS force field is used. In consequence a higher electrostatic field in the enzyme active site influences the proton position. 6. QCMD Simulations of the Proton Transfer. The proton in the hydrogen bond between the histidine and the water molecule was substituted with a one-dimensional Gaussian wavepacket located in the global energy minimum. The initial classical positions and momenta were taken from the configuration of the last part of the MD/VB simulation. The potential energy surface for the motion of the proton’s wavepacket is calculated on a discrete grid using the VB method in the same way as in the MD/VB simulations. The potential on the grid and the nonstationary extended Hellmann-Feynman forces are recalculated at each time step at new atomic positions. The proton wave function is propagated by integrating the timedependent Schroedinger equation using a Chebychev polynomial expansion method.6,49 Because the potential energy surface used in the QCMD/VB and MD/VB calculations is the same, the transition between the models is, according to the correspondence rule, smooth and no thermalization is required. The QCMD/VB simulation was run for 12 ps. During the first 1.0 ps of the QCMD/VB trajectory, the initial Gaussian wavepacket spreads and, due to its tails, the probability of finding the proton closer to the imidazole increases. Proton transfer occurs shortly thereafter. The probability distribution between different wells of the potential then becomes stabilized (Figure 14). The dividing surface for the probability integration is placed in the middle of the N-O bond at each time step. The observed proton transfer process is relatively fast. The probability of finding the proton close to the imidazole grows from 0.1 to 0.7 in a period of time in the range of 0.1-0.5 ps (Figure 14). Simultaneously, the probability of finding the proton close to the water oxygen decreases. The mean proton position 〈ψ(x)|x|ψ(x)〉 can be used for further analysis. The mean distance between the hydrogen and water oxygen atoms changes in this process from 1.05 to 1.50 Å. At the same time the distance between the proton and the imidazole ring decreases from 1.80 to 1.34 Å and a covalent bond is formed (Figure 14). The atomic charges of the imidazole-water complex calculated within the VB method exhibit a steplike change during the proton transfer process (Figure 15). In the water molecule, the charge of the oxygen atom changes from -0.745 e to -0.935 e. At the same time, the charge at the classical proton decreases and the OH- ion is formed with the total charge of -0.85 e. The charge on the imidazole atoms also change during the proton transfer event by 0.1 e, with a larger change at the nitrogen atom acting as the proton acceptor, which becomes
Figure 14. Probability of finding the proton close to the histidine nitrogen atom, i.e., probability of the proton transfer (a) and the mean bond length (O-H) in the water molecule (b). The results obtained within the QCMD/VB and the MD/VB models are compared (b) showing qualitatively different behavior.
Figure 15. VB charges of the atoms in the water molecule, histidine ring and a part of the substrate during the QCMD/VB dynamics.
more negatively charged by 0.193 e. The small change of the atomic charge on the water hydrogens is associated with the proton delocalization before the quantum hopping. Important changes in the potential energy surface for the proton dynamics can also be caused by the stretching of the
PT Processes in Molecular Complexes and Enzymes
Figure 16. Histidine nitrogen-water oxygen distance obtained within the QCMD/VB and the MD/VB models, respectively.
nitrogen-oxygen distance. Comparing the QCMD/VB and MD/VB trajectories (Figure 16), a slightly longer bond length is found in the quantum classical model (2.832 vs 2.751 Å). This histidine-water distance averaged over the proton transition period (2.847 Å) does not differ significantly from the mean distance. Therefore, the shortening of the nitrogen oxygen distance is not the crucial element promoting the proton transfer. Note that the MD/VB trajectories exhibit slightly higher stretching deformations of this bond than the QCMD/VB trajectories, but the reaction does not occur in the former case (Figure 16). The total energy of the proton wave function increases by 20 kcal/mol during the proton transfer event (see Figure 17). This increase results from the shift of the potential energy surface for the proton transfer caused by motions of the protein atoms. In addition, the shape of the potential is modulated by local changes of the R mode. In most cases, a significant increase in the proton transfer probability is associated with a preceding increase of the proton’s energy. Comparison of the energy of the proton wave function and the VB potential energy of the active site shows that the additional energy required for the proton hopping is transferred to the enzyme’s active site from the classical protein environment. Energy flow inside the protein does not affect conservation of the total energy of the molecules monitored during the simulation. The most important long-range interactions modifying the dynamics of the atoms in the active site is associated with the electrostatic interactions. In the VB method the electrostatic field generated by the charged atoms modifies the potential energy surface of the active site and therefore influences the whole process. This mechanism is also responsible for the energy transfer from the classical enzyme environment to the atoms in the active site. The magnitude of the electric field at the atoms was monitored during the QCMD/VB simulations (Figure 18). The most noticeable changes of the electric field
J. Phys. Chem., Vol. 100, No. 7, 1996 2543
Figure 17. VB energy of the proton wave function. The probability of the proton transfer is shown for comparison (bottom line).
Figure 18. Value of the electrostatic potential generated by the classical protein environment at selected atoms of the histidine ring, water molecule, and the substrate atoms during the QCMD/VB dynamics. The vertical lines show the time interval of the proton transfer.
take place at the water molecule, especially the hydrogen atoms. Changes of the electric field are also observed at the substrate oxygen O22, bonded to the carbonyl carbon. An increase of the electric field is observed at the O22 directly before the proton transfer. This suggests a significant role of this part of the substrate in the first step of the enzymatic reaction. The proton transfer is associated with a large change of the dipole moment
2544 J. Phys. Chem., Vol. 100, No. 7, 1996 of the hydrogen bond which significantly influences the electrostatic interactions with the rest of the system. However, as the simulations have shown, this process takes place during the proton transfer and does not influence interactions before this event. This suggests a mechanism similar to those observed in polar solvents, where proton transfer processes occur to preprepared product states.67-69 Creation of the pre-prepared product state induces the tunneling of the proton to the state localized in the metastable potential well located in the vicinity of the imidazole nitrogen atom. This spreads the proton density along the hydrogen bond, in a way not possible in the classical MD/VB model. After the proton transfer takes place, the OH- ion is formed and its motion toward the substrate carbonyl carbon C2P is observed. Due to the one-dimensional representation of the proton’s wave function and restrictions imposed on the wavepacket being collinear with its donor and the acceptor, the motion of the newly-formed OH- ion does not allow for a fully relaxed nucleophilic substitution. Work which would allow for the unrestrained motion of the OH- ion is in progress. 5. Conclusions A brief overview of quantum and quantum-classical molecular dynamics methods is presented. Two of them, the QCMD model,50,49 developed in our laboratories, and the Q/C TDSCF method41 with a well-established theoretical background, were applied to study proton transfer mechanism in the excited state of Xe-HI, and in the ground state of [H3N-H-NH3]+. The computer simulation results were compared with the accurate QD simulations. We have shown that the QCMD method in which coupling between the quantum particle(s) and the classical atoms is accomplished Via extended Hellmann-Feynman forces, as well as the time dependence of the potential energy function, reproduces sufficiently well the QD results. QCMD in comparison to the Q/C TDSCF conserves the energy somewhat better and gives better tunneling frequencies when compared to the reference QD results. The QCMD/VB model, in which the potential energy function is generated by the VB method, was applied to study proton transfer in phospholipase A2. Comparison of the QCMD/VB method with a conventional MD/VB approach shows the significant contribution of quantum effects, such as proton tunneling, to the enzymatic reaction. The QCMD model provides a new, efficient and practical tool to study enzymatic reactions with proper treatment of the quantum dynamical effects. Acknowledgment. This work was supported in part by the National Science Foundation, the Polish State Committee for Scientific Research, the Curie-Sklodowska Foundation, and the NSF Supercomputer Centers. Partial support from N. Copernicus University is also acknowledged. The calculations were partially performed in the Interdisciplinary Centre for Mathematical and Computational Modelling at Warsaw University. References and Notes (1) Dirac, P. A. M. Proc. Cambridge Philos. Soc. 1930, 26, 376. (2) Kramer, P.; Saraceno, S. Geometry of the Time-Dependent Variational Principle in Quantum Mechanics; Springer: New York, 1981. (3) Schroedinger, E. Ann. Phys. 1926, 79, 361, 489. (4) Schroedinger, E. Ann. Phys. 1926, 81, 109. (5) Leforestier, C.; Bisseling, R. H.; Cerjan, C.; Feit, M. D.; Freisner, R.; Guldberg, A.; Hammerich, A.; Jolicard, G.; Karrlein, W.; Meyer, H.D.; Lipkin, N.; Roncero, O.; Kosloff, R. J. Comput. Phys. 1991, 94, 59. (6) Truong, T. N.; Tanner, J. J.; Bala, P.; McCammon, J. A.; Kouri, D. J.; Lesyng, B.; Hoffman, D. J. Chem. Phys. 1992, 96, 2077. (7) Kosloff, R. Annu. ReV. of Phys. Chem. 1994, 45, 145.
Bała et al. (8) Levine, I. N. Quantum Chemistry, 2nd ed.; Allyn and Bacon: Boston, 1978. (9) Miller, W. H. In The Role of Computational Models and Theories in Biotechnology; Bertran, J., Ed.; NATO ASI Series; Kluwer Academic Publishers: Dordrecht, 1992; pp 193-236. (10) Makri, N. Chem. Phys. Lett. 1992, 193, 435. (11) Fernandez, F. M. Phys. ReV. A 1994, 50, 2953. (12) Jaroszewski, L.; Lesyng, B.; Tanner, J. J.; McCammon, J. A. Chem. Phys. Lett. 1990, 175, 282. (13) Jaroszewski, L.; Lesyng, B.; McCammon, J. A. J. Mol. Struct. (THEOCHEM) 1993, 283, 57. (14) Liu, X.; Murrel, J. N. J. Chem. Soc., Faraday Trans. 1991, 87, 435. (15) Aguado, A.; Paniagua, M. J. Chem. Phys. 1992, 96, 1265. (16) Case, D. A.; Cook, M.; Karplus, M. J. Chem. Phys. 1980, 73, 3294. (17) Hong, G.; Karplus, M. J. Chem. Phys. 1991, 94, 3679. (18) Aida, M.; Corongiu, G.; Clementi, E. Int. J. Quant. Chem. 1992, 42, 1353. (19) Vasilyev, V. V. J. Mol. Struct. (THEOCHEM) 1994, 304, 129. (20) Field, M. J.; Bash, P. A.; Karplus, M. J. Am. Chem. Soc. 1987, 109, 8092. (21) Luzhkov, V.; Warshel, A. J. Comput. Chem. 1992, 13, 199. (22) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (23) Field, M. J. J. Chem. Phys. 1992, 96, 4583. (24) Lunell, S.; Eriksson, L. A.; Fangstrom, T.; Maruani, J. Chem. Phys. 1993, 171, 119. (25) Liu, H.-Y.; Shi, Y.-Y. Chinese Sci. Bull. 1993, 38, 2023. (26) Skokov, S.; Carmer, C. S.; Weiner, B.; Frenklach, M. Phys. ReV. B, 1994, 49, 5662. (27) Carmer, C. S.; Weiner, B.; Frenklach, M. J. Chem. Phys. 1993, 99, 1356. (28) Maluendes, S. A.; Dupuis, M. Int. J. Quantum Chem. 1992, 42, 1327. (29) Adams, G. B.; Page, J. B.; Sankey, O. F.; Sinha, K.; Manendez, J.; Huffman, D. R. Phys. ReV. B 1991, 44, 4052. (30) Oguchi, T.; Sasaki, T. Prog. Theor. Phys., Suppl. 1991, 103, 93. (31) Blaudeck, P.; Frauenheim, T.; Porczag, D.; Seifert, G. Condens. Matter 1992, 4, 6389. (32) Dongqing, W.; Salahub, D. R. Chem. Phys. Lett. 1994, 224, 291. (33) Dongqing, W.; Salahub, D. R. J. Chem. Phys. 1994, 101, 7633. (34) Lynch, D. L.; Troullier, N.; Kress, J. D.; Collins, L. A. J. Chem. Phys. 1994, 101, 7048. (35) Alfonso, D. R.; Yang, S. H.; Drabold, D. A. Phys. ReV. B 1994, 50, 15369. (36) Car, R.; Parinello, M. Phys. ReV. Lett. 1985, 55, 2471. (37) Remler, D. K.; Madden, P. A. Mol. Phys. 1990, 70, 921. (38) Clarke, L. J.; Stich, I.; Payne, M. C. Comput. Phys. Commun. 1992, 72, 14. (39) Theilhaber, J. Phys. Fluids B (Plasma Phys.) 1992, 4, (Part 2), 2044. (40) Theilhaber, J. Int. J. Quantum Chem., Quantum Chem. Symp. 1994, 28, 611. (41) Alimi, R.; Gerber, R. B.; Hammerich, A. D.; Kosloff, R. J. Chem. Phys. 1990, 93, 6484. (42) Deumens, E.; Diz, A.; Longo, R.; Ohrn, Y. ReV. Modern Phys. 1994, 66, 917. (43) Longo, R.; Deumens, E.; Ohrn, Y. J. Chem. Phys. 1993, 99, 4554. (44) Berendsen, H. J. C.; Marvi, J. J. Phys. Chem. 1993, 97, 13464. (45) Mavri, J.; Lesnick, J.; Berendsen, H. J. C. Mol. Phys. 1994, 82, 1249. (46) Marvi, J.; Berendsen, H. J. C.; Van Gunsteren, W. F. J. Phys. Chem. 1993, 97, 13469. (47) Mavri, J.; Berendsen, H. J. C. J. Mol. Struct. 1994, 322, 1. (48) Bala, P.; Lesyng, B.; Truong, T. N.; McCammon, J. A. In The Role of Computational Models and Theories in Biotechnology; Bertran, J., Ed.; NATO ASI Series; Kluwer Academic Publishers: Dordrecht, 1992; pp 299326. (49) Bala, P.; Lesyng, B.; McCammon, J. A. Chem. Phys. 1994, 180, 271. (50) Bala, P.; Lesyng, B.; McCammon, J. A. Chem. Phys. Lett. 1994, 219, 259. (51) Bala, P.; Grochowski, P.; Lesyng, B.; McCammon, J. A. Comput. Chem., in press. (52) Florian, J.; Scheiner, S. J. Comput. Chem. 1994, 15, 553. (53) (a) Grochowski, P.; Bala, P.; Lesyng, B.; McCammon, J. A. Communication presented during the CECAM Workshop “Developing Hybrid Quantum and Classical Mechanical Methods for the Simulation of Biopolymers in Solution”, Lyon, France 9-11 May, 1995. (b) Grochowski, P.; Bala, P.; Lesyng, B.; McCammon, J. A. Submitted for publication in Int. J. Quant. Chem. (54) Bala, P.; Grochowski, P.; Lesyng, B.; McCammon, J. A. In Quantum Mechanical Simulations Methods for Studying Biological Systems; Field, M., Ed.; Les Houches Physics Series, 1995. (55) Alimi, R.; Gerber, R. B.; Hammerich, A. D.; Kosloff, R. J. Chem. Phys. 1990, 93, 6484.
PT Processes in Molecular Complexes and Enzymes (56) Scott, D. L.; White, S. P.; Otwinski, Z.; Gelb, M. H.; Sigler, P. B. Science 1990, 250, 1541. (57) White, S. P.; Scott, D. L.; Otwinski, Z.; Gelb, M. H.; Sigler, P. B. Science 1990, 250, 1560. (58) Scott, D. L.; Otwinski, Z.; Gelb, M. H.; Sigler, P. B. Science 1990, 250, 1563. (59) Bennion, C.; Connolly, S.; Gensmantel, N. P. J. Med. Chem. 1992, 35, 2939. (60) Sessions, R. B.; Dauber-Osguthorpe, P.; Campbell, M. M. Proteins 1992, 14, 45. (61) Rucker, J.; Cha, Y.; Jonsson, T. Biochemistry 1992, 31, 11489. (62) Northrop, D. B.; Duggleby, R. G. Bioorg. Chem. 1990, 18, 435. (63) Cha, Y.; Murray, C. J.; Klinman, J. P. Science 1989, 243, 1325.
J. Phys. Chem., Vol. 100, No. 7, 1996 2545 (64) Kim, H. J.; Hynes, J. T. J. Phys. Chem. 1990, 94, 2736. (65) Van Gunsteren, W. F. GROMOS (Groningen Molecular Simulation Computer Program Package); Biomos, Laboratory of Physical Chemistry, University of Groningen, 1987. (66) Biosym Technologies, San Diego. (67) Bruehl, M.; Hynes, J. T. Chem. Phys. 1993, 175, 205. (68) Borgis, D.; Hynes, J. T. Chem. Phys. 1993, 170, 315. (69) Borgis, D.; Hynes, J. T. J. Chem. Phys. 1991, 94, 361. (70) Jones, S. T.; Alhstrom, P.; Berendsen, H. J. C. Biochim. Biophys. Acta 1993, 1162, 135.
JP952642S
