help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Bala, P.
Right arrow Articles by McCammon, J. A.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Bala, P.
Right arrow Articles by McCammon, J. A.

Biophys J, September 2000, p. 1253-1262, Vol. 79, No. 3

Quantum-Dynamical Picture of a Multistep Enzymatic Process: Reaction Catalyzed by Phospholipase A2

P. Bała,*dagger P. Grochowski,* K. Nowinski,* B. Lesyng,*Dagger and J. A. McCammon§

 *Interdisciplinary Centre for Mathematical and Computational Modelling, Warsaw University, 02-106 Warsaw, Poland;  dagger Institute of Physics, Nikolaus Copernicus University, 87-100 Torun, Poland;  Dagger Department of Biophysics, Warsaw University, 02-089 Warsaw, Poland; and  §Howard Hughes Medical Institute, and Department of Chemistry and Biochemistry and Department of Pharmacology, University of California at San Diego, La Jolla, California 92093-0365 USA


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

A quantum-classical molecular dynamics model (QCMD), applying explicit integration of the time-dependent Schrödinger equation (QD) and Newtonian equations of motion (MD), is presented. The model is capable of describing quantum dynamical processes in complex biomolecular systems. It has been applied in simulations of a multistep catalytic process carried out by phospholipase A2 in its active site. The process includes quantum-dynamical proton transfer from a water molecule to histidine localized in the active site, followed by a nucleophilic attack of the resulting OH- group on a carbonyl carbon atom of a phospholipid substrate, leading to cleavage of an adjacent ester bond. The process has been simulated using a parallel version of the QCMD code. The potential energy function for the active site is computed using an approximate valence bond (AVB) method. The dynamics of the key proton is described either by QD or classical MD. The coupling between the quantum proton and the classical atoms is accomplished via Hellmann-Feynman forces, as well as the time dependence of the potential energy function in the Schrödinger equation (QCMD/AVB model). Analysis of the simulation results with an Advanced Visualization System revealed a correlated rather than a stepwise picture of the enzymatic process. It is shown that an sp2right-arrow sp3 configurational change at the substrate carbonyl carbon is mostly responsible for triggering the activation process.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

Quantum-dynamical processes, including proton transfer reactions, are of considerable importance in chemistry, physics, and biology (Gold and Caldin, 1975; Bamford, 1978; Bell, 1980; Warshel, 1991; Bountis, 1992), particularly for enzyme catalysis (Warshel and Levitt, 1976; Warshel, 1991; Klinman, 1989; Cha et al., 1989; Rucker et al., 1992; Bahnson and Klinman, 1995; Rucker and Klinman, 1999; Alhambra et al., 1999; Karsten et al., 1999). Current experimental studies show that classical transition state theory and its extensions accounting for tunneling corrections do not adequately describe an enzymatic reaction with proton transfer (Basran et al., 1999). Classical molecular dynamics simulations using either ab initio or semiempirical numerical quantum-mechanical potential energy functions have been performed for enzymatic reactions (e.g., Field et al., 1990; Warshel, 1991; Lee and Warshel, 1992; Alagona et al., 1996).

One of the most important problems in the study of enzymatic processes is our ability to describe their complex dynamics using models that are based on first physical principles. In particular, precise description of catalytic processes requires use of the time-dependent Schrödinger equation, which is difficult when modeling enzymatic reactions and dealing with large changes in the electronic charge distribution as well as proton or electron tunneling processes. Our previous quantum-classical molecular dynamics (QCMD) simulations for phospholipase A2, an enzyme that hydrolyzes phospholipids, were able to describe a proton tunneling process that occurs in the enzyme active site (Bala et al., 1995, 1996a,b), but the potential energy surface was not sufficiently flexible to describe the subsequent nucleophilic attack of the OH- group and cleavage of an ester bond. In the present study, the catalytic process is simulated using QCMD with an approximate valence bond method (AVB) (Grochowski et al., 1996) as a fast generator of the potential energy function, so that all steps of the catalytic process occur. The method is referred to as QCMD/AVB. The purpose of this study is to simulate the enzymatic process and to see if the results are consistent with the textbook description. One should note that a priori the method does not assume any specific reaction path. This study should also demonstrate that complex time-dependent biomolecular processes can effectively be studied using the proposed first-principles approach.


    SYSTEM
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

The secondary structure of phospholipase A2 (PLA2) consists of three connected alpha -helices. The active site is located between two parallel alpha -helices and includes a histidine (His48) and a water molecule. A Ca2+ cation coordinated by polar functional groups of the protein and water molecules also plays a role in hydrolysis. The crystal structure shows that His48 is hydrogen bonded to an aspartic acid residue (Asp93) and to structural water (Fig. 1).



View larger version (43K):
[in this window]
[in a new window]
 
FIGURE 1   Active site of phospholipase A2 and the substrate with the C2-O2 ester bond to be cleaved. The water molecule (W) playing a crucial role in the catalytic process is shown in the middle.

Diffusion of a substrate into the active site is the rate-limiting step of the whole enzymatic process. This aspect has been studied by others (Zhou and Schulten, 1996; H. Berendsen, personal communication). Based on crystallographic and biochemical data, a mechanism of hydrolysis was proposed (Scott et al., 1990b; Sessions et al., 1992; Waszkowycz et al., 1990). In the first step, a proton is transferred from the structural water molecule to the histidine residue, creating the OH- group. This allows for the nucleophilic attack of the latter group at the carbonyl carbon of the substrate. The intermediate product of the attack is the tetrahedral substrate-hydroxy anion adduct stabilized by the calcium ion (Sessions et al., 1992). In the next step of the reaction, the substrate ester bond is cleaved. Proton transfer from histidine to the product completes the catalytic cycle.

The proposed reaction scheme (Fig. 2) suggests a crucial role of the proton transfer process in the catalytic cycle, which is typical for a wide class of enzymatic reactions. During catalysis, bulk solvent has no access to the active site. The calcium ion is essential both for the binding of substrate and for catalysis (Sessions et al., 1992). It lowers the activation barrier of the transition state and controls the substrate binding (White et al., 1990; Scott et al., 1990a).



View larger version (11K):
[in this window]
[in a new window]
 
FIGURE 2   Conventional steps of the enzymatic reaction catalyzed by phospholipase A2.


    QCMD MODEL
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

A major challenge in the description of enzymatic reactions is the proper description of the complex bond-breaking process. This cannot be done properly by conventional classical molecular dynamics models. Description of the dissociation process requires nonstandard potential energy functions. Moreover, because of the low mass of the proton, quantum-mechanical tunneling may occur, as has been well documented in some enzymatic systems (Cha et al., 1989; Jonsson et al., 1994).

To meet these challenges, several methods have been developed. Most of them are based on partitioning of the system into quantum and classical subsystems, the evolutions of which are described by quantum and classical dynamics, respectively. The most successful implementations of these approaches are the time-dependent self-consistent method (Gerber et al., 1982), the surface hopping approach (Hammes-Schiffer, 1996), the density matrix evolution method (Marvi et al., 1993), and QCMD (see, e.g., Bala et al., 1992, 1994a, 1996a,b).

Approximation properties and limits of the QCMD model were analyzed by Bornemann et al. (1996). QCMD simulations results for model systems with proton transfer were compared with other approaches, particularly time-dependent self-consistent field (TDSCF) and full quantum-dynamical (QD) simulations, giving in the latter case very similar results (Bala et al., 1996a,b). The QCMD code was implemented on parallel computer architectures and successfully applied in studies of large biomolecular systems (Bala et al., 1997, 1998a,b).

In the current study, the enzyme molecule is 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 His48, a water molecule, and a segment of the substrate (-CH2-CO-O-CH2-) with the ester bond to be hydrolyzed (Fig. 1).

Let x and Xalpha denote position vectors of the quantum and classical degrees of freedom, respectively. The potential energy function is expressed as
V(<B><UP>x, X</UP></B><SUB>&agr;</SUB>)=V<SUP><UP>c</UP></SUP>(<B><UP>X</UP></B><SUB>&agr;</SUB>)+V<SUP><UP>q</UP></SUP>(<B><UP>x, X</UP></B><SUB>&agr;</SUB>), (1)
where Vc(Xalpha ) is the conventional, analytical potential energy term describing mutual interactions of atoms modeled as the classical particles, and Vq(xXalpha ) is the potential energy function of the quantum particle(s), also containing their interactions with the surrounding atoms.

The dynamics of the quantum particles is described by the time-dependent Schrödinger equation, and the rest of the system is governed by the Newtonian equations of motion with classical (Falpha c) and averaged, quantum Hellmann-Feynman (Falpha q) forces (Bala et al., 1994b).

In the simplest representation the QCMD equations are the following:
i&plank;<FR><NU>∂&PSgr;(<B><UP>x</UP></B>, t; <B><UP>X</UP></B><SUB>&agr;</SUB>)</NU><DE>∂t</DE></FR>=<FENCE>−<FR><NU>&plank;<SUP>2</SUP></NU><DE>2m</DE></FR> &Dgr;<SUB><UP>x</UP></SUB>+V<SUP><UP>q</UP></SUP>(<B><UP>x, X</UP></B><SUB>&agr;</SUB>)</FENCE> &PSgr;(<B><UP>x</UP></B>, t; <B><UP>X</UP></B><SUB>&agr;</SUB>), (2)

M<SUB>&agr;</SUB><A><AC><B><UP>X</UP></B></AC><AC>¨</AC></A><SUB>&agr;</SUB>=<B><UP>F</UP></B><SUP><UP>c</UP></SUP><SUB>&agr;</SUB>+<B><UP>F</UP></B><SUP><UP>q</UP></SUP><SUB>&agr;</SUB>, (3)
where
<B><UP>F</UP></B><SUP><UP>c</UP></SUP><SUB>&agr;</SUB>=<UP>−</UP><FR><NU>∂</NU><DE>∂<B><UP>X</UP></B><SUB>&agr;</SUB></DE></FR>V<SUP><UP>c</UP></SUP>(<B><UP>X</UP></B><SUB>&agr;</SUB>), (4)

<B><UP>F</UP></B><SUP><UP>q</UP></SUP><SUB>&agr;</SUB>=<UP>−</UP><FR><NU>∂</NU><DE>∂<B><UP>X</UP></B><SUB>&agr;</SUB></DE></FR>⟨&PSgr;(<B><UP>x</UP></B>, t; <B><UP>X</UP></B><SUB>&agr;</SUB>)‖V<SUP><UP>q</UP></SUP>(<B><UP>x, X</UP></B><SUB>&agr;</SUB>)‖&PSgr;(<B><UP>x</UP></B>, t; <B><UP>X</UP></B><SUB>&agr;</SUB>)⟩. (5)
The coupling terms determine the energy transfer between the quantum and classical domains, which influences the time evolution of the whole system (Bala et al., 1994a,b, 1996a,b). In this study the Hellmann-Feynman forces are evaluated using a spectral decomposition of the wavefunction into instantaneous eigenstates (Bala et al., 1994b). Molecular dynamics with the Hellmann-Feynman forces works very well in the situation where the quantum particles are localized in a compact area, which is the case for the studied proton transfer process in the enzyme active site. The forces between quantum particle(s) and classical atoms located sufficiently far from the quantum subsystem can be calculated in the classical limit. In this study, the classical forces between the proton charge located in the mean proton position and remote classical atoms have been used instead of the Hellmann-Feynman forces.


    AVB METHOD
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

Evaluation of the potential for the quantum domain Vq (x, Xalpha (t)) is a bottleneck of all theoretical calculations of enzymatic reactions. Determination of high-quality potential energy surfaces (Ohrn et al., 1996) requires usage of advanced quantum mechanical methods. However, reliable ab initio or density functional theory (DFT) methods are very time consuming and cannot be directly applied in hybrid quantum-classical models. However, reliable ab initio or DFT methods are very time consuming and are hard to apply to whole enzymatic processes (see, e.g., Parrinello, 1997).

To overcome this difficulty we have formulated and parameterized an approximate valence bond (AVB) method (Grochowski et al., 1996) to describe the potential energy functions for proton transfer and hydroxy anion nucleophilic reactions in enzymes. AVB is similar to the EVB approach (Warshel, 1991; Lee and Warshel, 1992) but uses ab initio rather than empirical parameterization. In particular, this parameterization has been based on DFT calculations with the Heiden-Lundqvist/Janak-Morruzi-Williams functional with the DNP basis. The most crucial calculations, including the energy profiles for the proton transfer processes, have been repeated with the conventional ab initio method at the MP2-6-31(d) level, using the Gaussian program. The two methods gave very similar energy profiles. Dissociation energies of isolated molecules were verified against available experimental data in the vapor state and scaled to the precise experimental observables (proton affinities: 390 kcal/mol for OH-, 229 kcal/mol for imidazole, and 381 kcal/mol for CH3O-). Influence of the molecular environment on the energy profiles was accounted for via the electrostatic field generated by the atomic charges. The AVB method reproduces bond dissociation processes as well as partial atomic ESP charges, which allows for proper description of the electrostatic potential in the area of the quantum site (Tables 1-3).


                              
View this table:
[in this window]
[in a new window]
 
TABLE 1   Parameters of the modified Lennard-Jones potentials according to Shiratori and Nakagawa (1991), describing interactions of the water molecules with the Ca2+ ion


                              
View this table:
[in this window]
[in a new window]
 
TABLE 2   The AVB atomic charges for the initial geometry


                              
View this table:
[in this window]
[in a new window]
 
TABLE 3   The atomic charges in selected groups of the substrate molecule

The AVB method has been applied in calculations of the potential energy surface for selected atoms as well as the proton of the water molecule in the phospholipase A2 active site. The quantum region of the PLA2-substrate complex is described by 14 different VB structures (see Fig. 3).



View larger version (17K):
[in this window]
[in a new window]
 
FIGURE 3   Structures used in the AVB model of the enzymatic reaction catalyzed by phospholipase A2.

The instantaneous, microscopic Born-Oppenheimer potential energy surface for the proton motion in the isolated active site is bistable, with a deep minimum located 0.98 Å from the water oxygen. The second minimum, near the imidazole nitrogen, is ~40 kcal/mol above the global one (Grochowski et al., 1996). The presence of the structural Ca2+ ion in the active site reduces this barrier by up to 20 kcal/mol. The structural as well as electrostatic field fluctuations further reduce this barrier (see, e.g., Bala et al., 1995; Grochowski et al., 1996), and therefore the mean barrier for the proton transfer is much lower. Changes in the shape of the potential energy surface play a crucial role in the proton quantum dynamics and lead to the proton transfer and, later on, to the successive reaction steps (cf. Warshel, 1984; Hwang et al., 1991).


    MD/AVB AND QCMD/AVB MODELS
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

The AVB method is used to generate forces acting on the atoms in the active site in both classical (MD/AVB) and quantum-classical molecular dynamics (QCMD/AVB) simulations. The AVB potential energy surface also accounts for changes in hybridization of atoms in the active site, in particular water oxygen, His48 nitrogen, and the substrate carbon atom (see Fig. 1).

The conventional interaction parameters used in the simulations are based on the GROMOS (Van Gunsteren, 1987) parameterization with implicit solvent. Several modifications of the electrostatic charges of atoms in the neighborhood of the Ca2+ cation were introduced as in our previous studies (Bala et al., 1994b, 1995). Atomic charges on the water molecule, the imidazole ring of His48, and on the substrate directly involved in hydrolysis are obtained within the AVB method and recalculated at the each time step. The charges on the substrate PO4- group and on the terminal NH3 group were obtained from DFT calculations, using DMol (Molecular Simulations) and are kept fixed during the dynamics.

The total potential energy surface used in the MD/AVB and QCMD/AVB models is given by Eq. 2 with
V<SUP><UP>q</UP></SUP>(<B><UP>x, X</UP></B><SUB>&agr;</SUB>)=V<SUP><UP>AVB</UP></SUP>(<B><UP>x, X</UP></B><SUB>&agr;</SUB>). (6)
The Vc potential is calculated using the GROMOS empirical potential energy functions and describes interactions within the classical region as well as the bonding and Lennard-Jones interactions between the classical and quantum regions. The AVB component is the total Born-Oppenheimer energy resulting from the AVB calculations for the active site region, including the electron polarization effects inside this site as well as the electrostatic interactions with the classical region (Grochowski et al., 1996; Bala et al., 1996a,b, 1998a,b). x is the position vector of the proton moving from water to imidazole (Fig. 1). As mentioned before, the charges of the atoms in the quantum region associated with the enzymatic reaction and other electrostatic properties are evaluated at each time step according to the AVB model.


    THE ENZYME SYSTEM: ITS THERMALIZATION AND EQUILIBRATION
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

In our studies, crystal structures of both native (Scott et al., 1990b) and inhibited (White et al., 1990; Scott et al., 1990a; Bennion et al., 1992) enzymes from a cobra venom were used. The inhibitor molecule was changed into the substrate by replacing the phosphonate group of the transition state analog by a carbon atom and one water molecule. Before the MD/AVB simulations, the structure of the system was relaxed. One thousand steps of steepest-descent minimization were executed, followed by 10 0.1-ps dynamics runs at a temperature of 5 K performed in the microcanonical ensemble (NVT), with a time step of 0.5 fs. The final coordinates were chosen as a starting point for further energy minimization. One thousand one-step runs at a temperature of 5 K were performed, and the velocities were reassigned at each time step to remove excess kinetic energy 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) was performed at 300 K. The last stage of the thermalization procedure consisted of 55-ps dynamics in the NVT ensemble with a relaxation time of 1 ps. During the thermalization and equilibration process a restraint on the water oxygen-hydrogen distance was applied to prevent proton transfer to the imidazole ring. All simulations were performed with the MD/AVB model.


    PROTON TRANSFER AND HYDROLYTIC PROCESSES IN THE MD/AVB AND QCMD/AVB SIMULATIONS
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

The productive MD/AVB trajectory started at 150 ps. 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. At the 224-ps point of the simulation, the water molecule was stable and the OH bond length was 1.02 Å, a typical value for water forming a hydrogen bond. At 224.5 ps a change in the substrate carbon C2 hybridization from planar to tetrahedral is observed (Fig. 4). Shortly thereafter, at 224.82 ps, the proton transfer from the water molecule to the imidazole ring of His48 takes place.



View larger version (19K):
[in this window]
[in a new window]
 
FIGURE 4   Changes in the carbonyl carbon dihedral angle C2-O2-O22-C22 in the substrate molecule. MD/AVB simulation.

The changes in the active site induced by the enzymatic environment are confirmed by the time dependence of the formal charges of the atoms in the active site. The charges were calculated at each time step within the AVB model. The C2 carbon charge increases by 0.2e at 224.82 ps, before the proton transfer takes place, and at the same time the hybridization change from sp2 to sp3 at the carbon atom is observed. The charge transfer and the change in hybridization are induced by changes in the electrostatic field generated by the protein environment and electron polarization effects inside the active site induced by the proton transfer.

At 224.85 ps the mean distance between the hydrogen and water oxygen increases to 1.7 Å. At the same time, the distance between the proton and imidazole ring nitrogen decreases to 1.0 Å, and a covalent bond is formed (Fig. 5 a). The proton transfer is fast and is completed within 0.02 ps. However, the stabilization of the bond and dissipation of the excess energy of the proton take a longer time (0.2 ps).



View larger version (26K):
[in this window]
[in a new window]
 
FIGURE 5   The water oxygen (OW)-hydrogen (gray trace) and imidazole nitrogen (ND1)-hydrogen (black trace) distances obtained in the MD/AVB (a) and QCMD/AVB (b) simulations.

The proton transfer process is associated with a slight increase in the water oxygen-imidazole nitrogen distance. This observation is contrary to the common conviction that shortening of the hydrogen bond, which lowers the energy barrier for the proton motion, is one of the most important promoters of the proton transfer. Such correlations between the proton transfer process and length of the hydrogen bond have not been observed in these simulations.

While the proton transfer takes place, the OH- ion is formed, and instantaneously it begins moving toward the substrate carbonyl carbon (Fig. 6 a). Decreases in the OH--C2 substrate distance, first to 1.9 Å and then after 0.2 ps to 1.6 Å, are observed. The cooperative motion of the proton and the OH- group does not support the textbook sequential mechanism of the enzymatic reaction presented in Fig. 2. Although the first, most significant part of the nucleophilic attack is observed at the same time as the proton transfer, the stabilization of the formed bond takes ~0.6 ps. These processes are associated with the increase in the mean C2-O2 distance from 1.5 Å to 1.75 Å (Fig. 6), and the covalent bond C2-O2 is cleaved. One can say that the enzymatic process occurs in a concerted manner. Such processes, among others, were discussed by Jencks (1997).



View larger version (21K):
[in this window]
[in a new window]
 
FIGURE 6   Change in the water oxygen (OW)-carbon (C2) distance (gray trace) and the length of the substrate ester bond (C2-O2) (black trace), obtained in the MD/AVB (a) and QCMD/AVB (b) simulations.

The simulation shows that successful proton transfer takes place only when it is accompanied by the nucleophilic attack and the change of hybridization at the carbon atom, which is visible in the changes of the dihedral angles involving the C2 substrate atom. Because the proton transfer and the nucleophilic attack take place after the C2 hybridization starts to change, one sees that changes at the substrate carbon atom are crucial for the reaction. One should stress that an additional 50-ps MD/AVB simulation with a restraining potential forcing the planar conformation of the substrate C2 atom did not lead to any proton transfer, and during this simulation no catalytic process was observed.

More accurate simulations were also performed using the QCMD/AVB model, allowing the proton to be a fully quantum-dynamical particle. In this case, the proton in the hydrogen bond between the histidine and the water molecule was replaced with a three-dimensional Gaussian wavepacket located in the global energy minimum. The initial classical positions and momenta were taken from the MD/AVB at 224.4 ps, before the classical proton transfer process occurred. The momentum of the wavepacket was set to zero. The potential energy surface for the motion of the proton's wavepacket was calculated on a discrete 64 × 32 × 32 point grid with the AVB method. The potential on the grid and the Hellmann-Feynman forces were recalculated at each time step at new atomic positions. The proton wavefunction was propagated on the discrete grid by numerical integration of the time-dependent Schrödinger equation, using a Chebychev polynomial expansion method (Bala et al., 1994a; Tel-Ezer and Kosloff, 1984; Truong et al., 1992).

Because the potential energy surfaces used in the QCMD/AVB and MD/AVB calculations are the same one, the transition between the models is smooth and no thermalization is required according to the correspondence rule. The QCMD/AVB simulation was run for 0.7 ps. The time step, temperature, and coupling to the thermal bath were the same as in the preceding MD/AVB simulations.

The mean distance between the hydrogen and water oxygen atoms changes during the first 0.2 ps from 1.05 Å to 1.20 Å (Fig. 5). This effect is not observed in the reference classical trajectory and is a signature of the quantum nature of the proton in the hydrogen bond (Borisenko et al., 1995). In the same time period, a shortening of the water oxygen-substrate carbon distance from 2.2 Å to 1.95 Å is observed (Fig. 6 b). Shortly thereafter the proton transfer occurs. The mean distance between the hydrogen and water oxygen increases to 1.9 Å. At the same time the distance between the proton and nitrogen of the imidazole ring decreases to 1.0 Å, and a covalent bond is formed (Fig. 5 b). The proton transfer process occurs in a smoother way than the classical one and is completed within 0.2 ps. During this time the proton's wavefunction is strongly delocalized. After this period the delocalization and splitting of the proton's density between the two potential minima disappear. The wave function localizes near the imidazole ring. The observed time of the proton transfer event is similar to the value obtained in the MD/AVB simulations; however, in the QCMD/AVB model the process is smooth and cannot be decomposed, as in the previous case, into fast proton transfer and slow vibronic relaxation.

The quantum-dynamical mechanism of the proton transfer reaction can be characterized by analysis of the changes in the potential energy surface together with changes in the proton's probability density. The potential energy surface for the proton motion changes significantly during the QCMD/AVB simulations. The energy barrier for the proton transfer is initially 20 kcal/mol and decreases during the process. The barrier almost vanishes at time 224.48 ps, and shortly afterward (0.01 ps) a rapid increase in the proton transfer probability is observed (Fig. 7). Before the proton transfer, the minimum located near the water molecule is 10 kcal/mol lower than the second one. At time 224.45 ps, just after the proton transfer process is started, the energy minimum close to nitrogen is lowered by 25 kcal/mol and becomes lower than the other one (Fig. 7). At all times, the proton's wavefunction is strongly localized in the potential minimum near the water molecule, and the probability of finding a proton near nitrogen is lower than 0.1. At 224.45 ps, the minimum near water oxygen becomes higher than the other one by 10-15 kcal/mol, and after 0.15 ps it cannot be resolved (Fig. 7). At 224.6 ps the potential energy surface has only one minimum, close to the imidazole nitrogen. One should note that most of the proton transfer process takes place during this period, and ~60% of the proton's density is transferred to the nitrogen atom. Full localization of the proton's wavepacket near imidazole is achieved after 0.4 ps from the time at which the minimum close to imidazole becomes deeper. At the same time additional stabilization of the potential minimum is observed. Because the energy minimum level is at this time clearly correlated with the proton's total density in the area near the imidazole nitrogen, this stabilization can be treated as an accommodation of the active site region to the situation after the proton transfer and C-O bond cleavage.



View larger version (18K):
[in this window]
[in a new window]
 
FIGURE 7   The height of the energy barrier for the proton transfer (EB) and the energy difference between the potential minima, one close to the water molecule and another close to the imidazole nitrogen (Delta Emin). The integrated probability of the proton being found in the area close to imidazole is also presented.

Effective barriers for particle transfer can differ noticeably, depending on whether the particle is treated classically or quantum-dynamically. This has clearly been shown by Parrinello and co-workers for the hydrated excess proton in water. It appears that a barrier with a height of ~2 kcal/mol is washed out to zero after the classical proton is replaced by the quantum-dynamical one (Tuckerman et al., 1997).

We observed a similar effect in our study. By applying the same averaging procedure as described in the paper above, we computed the effective barrier for the proton transfer, first treating proton classically (AD/AVB), and then treating the proton quantum-dynamically (QCMD/QVB). The averaging window included the proton transfer process and cleavage of the ester bond. In the classical case we got a barrier height of 2.5 kcal/mol, and in the quantum one, 0.5 kcal/mol. This phenomenon is due to the fact that the quantum particle is always localized above the zero-vibrational level, and in addition the delocalized probability density probes areas with lower energy through the instantaneous barrier "smoothing." For such cases the simple transition state theory does not apply.

A better understanding of the catalytic process can be achieved by visualizing in three dimensions the structural changes of the enzyme active site containing the substrate molecule, along with the changes in the electrostatic potential and changes in the total potential energy function for the proton transfer that drives the evolution of the proton probability distribution. An AVS/Express environment (Advanced Visual Systems, Waltham, MA) was used to present and analyze selected snapshots of the QCMD/AVB simulation. Fig. 8, a-d, presents the proton transfer process from the water molecule to the histidine ring. Over time the proton charge distribution loses its initial convex shape. The transition state is characterized by the split and roughly symmetrical distribution being shared by the donor and the acceptor (see Fig. 8 c). The equipotential energy surfaces surrounding the proton charge distribution evolve in time in a smooth way. Its bistable shape changes, shifting the global energy minimum from the water oxygen atom to the imidazole. This shift results mostly from the electrostatic interactions, which are apparent when the electrostatic potential field is analyzed. The proton transfer process initiates the nucleophilic attack of the OH- group on the substrate carbonyl carbon. Fig. 8, c-e, presents the attack. The OH- nucleophilic attack starts before the proton transfer is fully completed. In turn, as already mentioned, the OH- attack is strongly coupled to hydrolysis of the C2-O2 ester bond. Fig. 8, d and f, presents configurations characteristic for the beginning and the end of the hydrolytic process. One sees the sp2right-arrow sp3 hybridization change at the C2 carbon (Fig. 8 d), which initiates the hydrolytic process. The distances OH--C2 and C2-O2 are equal, respectively, to 1.81 Å and 1.66 Å at the beginning (Fig. 8 e) and to 1.57 Å and 1.74 Å at the end of this process.



View larger version (62K):
[in this window]
[in a new window]
 
FIGURE 8   Selected snapshots of the QCMD/AVB simulation. Figures present the proton transfer process from the water molecule to the histidine ring (a-d), the nucleophilic attack of the OH- group (c-e), and hydrolysis of the substrate C2-O2 bond (e-f). The potential energy surface for the proton motion is represented by a contour plot, and color at the plane denotes total amplitude of the electrostatic potential generated by the protein.


    CONCLUDING REMARKS
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

In this study, for the first time, we were able to use time-dependent quantum theory and a quantum approximation for the interatomic forces to simulate the whole catalytic process of an enzyme. In particular, the quantum-dynamical dissociation of the water molecule in the enzyme active site and the nucleophilic attack of the ensuing OH- group at the ester bond of the substrate were simulated. The results demonstrate the capability of the MD/AVB and QCMD/AVB models for the study of enzymatic reactions (Fig. 9). The QCMD approach is more accurate and automatically accounts for the zero-level vibrational correction, a characteristic of the time-independent quantum theories that always keeps the wave packet above the potential energy minimum and increases the transfer probability (see, e.g., Tuckerman et al., 1997, Benoit et al., 1998). Changes in the potential energy surface driving the reaction process are determined by the changes in the atomic positions in the active site and by the electrostatic potential generated in the active site region by the protein environment.



View larger version (20K):
[in this window]
[in a new window]
 
FIGURE 9   The effective energy barriers for the classical (MD/AVB) and quantum-dynamical (QCMD/AVB) proton transfer processes. A dimensionless reaction coordinate q is defined through proton-oxygen and proton-nitrogen distances as q = (rOH - rNH)/(rOH + rNH).

The results show that shortening of the hydrogen bond is not the only process that can promote the proton transfer. It is clearly seen that the protein scaffold, especially fluctuations in the electrostatic field generated by the charged residues, plays a significant role in the catalytic cycle. The most important are structural and electronic changes occurring in the substrate molecule, especially at the carbonyl carbon atom C2.

A second important implication of the present results is that the nucleophilic attack and the proton transfer from water to His48 are almost simultaneous. Moreover, successful proton transfer cannot take place when the nucleophilic attack is prohibited by constraining the carbonyl carbon at the sp2 hybridization. Both processes require proper preparation of the active site region. In particular, the successful reaction cycle must be preceded by the charge shift at the substrate atoms, especially carbon C2. This process must take place before the proton transfer to stabilize the reaction products.

    ACKNOWLEDGMENTS

This work was supported by the Polish State Committee for Scientific Research (8T11F 016 16), Nikolaus Copernicus University, the San Diego Supercomputer Center, and the U.S. National Science Foundation.

    FOOTNOTES

Received for publication 24 February 2000 and in final form 5 June 2000.

Address reprint requests to Dr. J. Andrew McCammon, Department of Chemistry and Biochemistry, University of California-San Diego, Urey Hall, Room 4238, 9500 Gilman Drive, La Jolla, CA 92093-0365. Tel.: 858-534-2905; Fax: 858-534-7042; E-mail: jmccammon{at}ucsd.edu.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
SYSTEM
QCMD MODEL
AVB METHOD
MD/AVB AND QCMD/AVB MODELS
THE ENZYME SYSTEM: ITS...
PROTON TRANSFER AND HYDROLYTIC...
CONCLUDING REMARKS
REFERENCES

Biophys J, September 2000, p. 1253-1262, Vol. 79, No. 3
© 2000 by the Biophysical Society   0006-3495/00/09/1253/10  $2.00



This article has been cited by other articles:


Home page
Protein Sci.Home page
J. Trylska, P. Grochowski, and J. A. McCammon
The role of hydrogen bonding in the enzymatic reaction catalyzed by HIV-1 protease
Protein Sci., February 1, 2004; 13(2): 513 - 528.
[Abstract] [Full Text] [PDF]


Home page
Biophys. JHome page
J. Trylska, P. Bala, M. Geller, and P. Grochowski
Molecular Dynamics Simulations of the First Steps of the Reaction Catalyzed by HIV-1 Protease
Biophys. J., August 1, 2002; 83(2): 794 - 807.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
M. A. Lill and V. Helms
Proton shuttle in green fluorescent protein studied by dynamic simulations
PNAS, March 5, 2002; 99(5): 2778 - 2781.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Bala, P.
Right arrow Articles by McCammon, J. A.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Bala, P.
Right arrow Articles by McCammon, J. A.