11266
J. Phys. Chem. 1995,99, 11266- 11275
Dynamic Force Field Models: Molecular Dynamics Simulations of Human Carbonic Anhydrase I1 Using a Quantum MechanicalMolecular Mechanical Coupled Potential David S. Hartsough and Kenneth M. Merz, Jr.* Department of Chemistry, I52 Davey Laboratory, Pennsylvania State University, University Park, Pennsylvania I6802 Received: February 7, 1995; In Final Form: April 13, 1 9 9 9
We report extended molecular dynamics simulations of the enzyme human carbonic anhydrase I1 using a quantum mechanicaVmolecular mechanical coupled potential. This method allows us to treat the key active site residues of this metalloenzyme using semiempirical quantum mechanics, while directly incorporating the effects of the surrounding enzymatic environment into the electronic structure calculations. Furthermore, this approach provides complete geometric freedom within the zinc coordination sphere and also enables examination of the relative importance of geometric and environmental effects upon atomic charge.
Introduction Computer simulation of proteins and other macromolecules has proven to be an invaluable aid in understanding their structural and dynamic proper tie^.'-^ Because of the importance of understanding these systems, much effort has been directed toward refining the computational methods used to represent them. However, considerations of computational expense severely limit the degree of rigor which can be applied to any system. As a consequence of such considerations, molecular mechanical force fields are typically utilized in the study of macromolecular systems in s ~ l u t i o n . ~However, -~ using such methods implicitly eliminates the possibility of studying inherently quantum mechanical events such as bond breaking or charge reorganization. Although methods have been proposed which incorporate terms for such events in a MM force field,7.8 these treatments will by their very nature be only an approximation to the quantum mechanical reality and often ignore the perturbing effect of solvendenzyme on the electronic structure of the reactive substrate. Other approaches use the so-called supermolecule approach and study enzymatic reactions in the gas phase.9-' While this approach is computationally attractive and is capable of giving mechanistic insights, it is not capable of giving an accurate accounting of environmental effects that arise within an inhomogeneous environment like an enzyme active site.'3-'5 A method which potentially affords an escape from this difficulty is a quantum mechanicaYmolecular mechanical (QM/ MM) coupled potential.'6-20 In such a method, the computational model is partitioned into QM and MM regions (Figure 1). Each portion is then treated using the appropriate computational method, and communication between the two regions is facilitated via van der Waal's and electrostatic interactions. Using such a method, it is thus possible to carry out a quantum mechanical examination of bond breaking/forming events in an explicit solvendenzymatic environment. The QM method to be employed may vary from empirical valence semiempirical molecular orbital,I6.I7density f u n ~ t i o n a l ,and ~~.~~ even Hartree-Fock methods. 1924 For reasons of computational expense, NDDO-based semiempirical quantum mechanical methods have seen the most widespread use to date.'7.'8.25-3' However, it is important to note that there is nothing hindering the use of more sophisticated quantum mechanical treatments, @
Abstract published in Advance ACS Abstracts. June 15, 1995.
0022-3654/95/2099-11266$09.00/0
1%K
Figure 1. Division of the molecular system within a coupled potential calculation. The quantum mechanical (QM) region is surrounded by a molecular mechanical (MM) region and then a boundary region (BR).
other than the obvious demands placed upon CPU time. The utility of these methods for the study of quantum mechanical processes such as charge reorganization and chemical reactions in aqueous solution has been amply demonstrated by the work of K a r p l u ~ Ga0,27.30 ,~~ and W a r ~ h e l . ~ The representation of metal ions in force fields is difficult because of the large charge associated with many ions (e.g., Zn2+,Fe3+,etc.) which results in substantial polarization effects, the ability of metal ions to undergo coordination number changes, as well as the ability to undergo facile rearrangements (e.g.,Berry pseud~rotation).~~ Several limiting strategies have been employed to study metal ions in metal lo enzyme^.^^.^^ One set of approaches, the so-called "bonded" approach involves placing explicit bonds between a metal ion and its surrounding ligands,32 while the nonbonded approach involves the use of electrostatic and van der Waals contacts to retain the shape of the coordination sphere around the metal ion.33.34The bonded approach has as its major limitation the inability to allow for coordination changes around the metal ion. On the other hand, the nonbonded approach is very sensitive to the electrostatic model chosen33and can suffer from the inability to retain low coordination number^.^' Coupled potential methods can offer a solution to these problems by treating the metal ion and its ligands quantum mechanically while treating the surrounding environment with a force field. Using this "dynamic" approach, one can effectively deal with the charge representation since polarization of the metal ion and its ligands is effectively included. Furthermore, the dynamics of the coordination sphere is included as well as the ability to undergo coordination
0 1995 American Chemical Society
QM-MD Simulations of Human Carbonic Anhydrase I1
SCHEME 1 Rgn-OH
J. Phys. Chem., Vol. 99, No. 28, 1995 11267 SCHEME 2
coz
RgnOH +COP
Rgn-OH2
HCO,
L
T-
:NS# %s
E binding site €4
0
3CP
changes. Determining whether this is indeed the case is one of the interests of this work. Our long-term interest in using QMlMM coupled potentials is as a tool in our effort to understand enzymatic reactivity. In particular, we are interested in understanding the molecular level details in the catalytic cycle of the enzyme human carbonic anhydrase I1 (HCAII).35-37 HCAII, one of seven isozymes of the zinc metalloprotein human carbonic anhydrase (HCAL3*is a 260-residue protein with a mass of -29 kDa. A single zinc atom is located in the enzyme active site and is necessary for catalytic a ~ t i v i t y . ~ ~The - ~ ' active site itself lies at the bottom of a deep cavity (15 A.deep) in the protein, which is readily accessible to solvent.39 The active site cavity is divided into hydrophobic and hydrophilic regions, with a network of hydrogen-bonded water molecules connecting the active site region and the surrounding solvent e n ~ i r o n m e n t .The ~ ~ catalytically necessary zinc atom lies at the bottom of the active site cleft and is tetrahedrally coordinated by three histidine residues (His-94, -96, and -1 19) and a fourth ligand,39whose identity is pH-dependent. At high pH ('8) the fourth ligand is a hydroxide ion, while at acidic pH the fourth coordination site is occupied by a water m ~ l e c u l e . ~ ~ - ~ ~ The catalytic mechanism of HCAII has been studied in detail, yet much still remains ~ n c l e a r . ~Catalysis ~ - ~ ~ has been found to depend upon a group with a pKa of around 7, with the fourth zinc ligand (hydroxide/water) appearing to fulfill this requirement. This consideration, in conjunction with the observed ping-pong kinetics, gave rise to the mechanism shown in Scheme 1.37 The proton transfer step converting D into E has been implicated as the rate-limiting step at high concentration of extemal buffer, while E to A is thought to be rate limiting a low buffer concentration^.^^-^^ The conversion of D to A is kinetically distinct from the sequence of steps converting A into D, via B and C. Although experiments have not completely elucidated the detailed structural changes in the mechanism for catalysis, there is considerable evidence that certain residues are catalytically important (see Scheme 2). These include His64, Glu-106, Thr-199, and several water molecules near the active site.35-37 Thr-199 is positioned with Thr-200 on the opposite side of the active site cavity from the zinc atom. These threonine residues, His-64 (located at the entrance of the active site cavity), and Glu-106 combine with other polar residues to constitute the hydrophilic half of the cavity. Thr-199 is a key residue which is centered between the hydrophilic and hydrophobic halves of the cavity. It is locked into this position as part of a key hydrogen-bonding network. The hydrogen of the hydroxolaquo zinc ligand is donated for hydrogen bonding to
Ales
H2NYo Oh44
the y-oxygen of Thr-199. The proton of the hydroxyl g o u p on Thr-199 is then in tum donated to an +oxygen of Glu-106, forming a second hydrogen bond interaction (see Scheme 2). The proximity of Thr-199 to the active site zinc atom and the rigidity of this hydrogen bond network are considered to be crucial for catalysis and inhibitor b i n d i r ~ g . ~This ~ . ~group ~ of hydrogen bonds may also play an important role in CO? binding and catalysis. It has been suggested that this hydrogen bond network serves to properly orient the lone pair electrons of the hydroxide ligand allowing for rapid addition to CO?, as the CO2 molecule approaches from the hydrophobic cavity (A binding site in Scheme 2).40-42 Bordering on this hydrogen bond network is a group of eight water molecules which extends toward bulk water. It has been proposed that these water molecules serve to shuttle a proton out of the active site and into bulk solution via H i ~ - 6 4 . ~ ~ Although our stated long-term goal is an examination of the entire HCAII catalytic cycle using coupled QM/MM methods, we must first establish that this methodology is able to accurately represent the active site environment and reproduce key interactions in this region. Besides examining gross structural changes in the active site region, we are also interested in examining the structure and dynamic properties of the aforementioned hydrogen bond network between ZnOH(H2), Thr-199, and Glu106. In addition to providing an opportunity to examine the capabilities of our coupled QM/MM method, the calculations reported herein will also allow us to evaluate the importance of quantum mechanical events, such as charge reorganization, which have been ignored in many prior simulations of metalloenzyme~.~~ Computational Methods QM/MM Coupled Potential. The QM/MM method employed is essentially identical to the method employed by Field, Bash, and Karplus.I6 Because they have already outlined the theoretical basis and technical details of this method, we will only comment upon some significant differences between our model and theirs. The coupled potential we have implemented couples together the MD program AMBER 4.043 and the
Hartsough and Merz
11268 J. Phys. Chem., Vol. 99, No. 28, I995 semiempirical quantum mechanical program MOPAC 5.0. We have used standard AMBER force field parameters throughexcept for the active site region, where we have used the MM parameters of Hoops e f ai. developed especially for HCAII.32 We have also found it necessary to add a special hydrogen-bonding 10- 12 interaction between the protons of the fourth zinc ligand (water or hydroxide) and the hydroxyl oxygen of Thr-199 (parameters for this interaction C = 4019 and D = 21004.5). Similar to the work of Karplus and c o - w ~ r k e r s ’and ~ . ~Singh ~ and K ~ l l m a n , we ’ ~ have used link atoms to cap exposed valencies due to bonds which cross the QM-MM boundary. In this method the QM region of the system is treated as a closed-shell molecule with no exposed valencies. In our system the imidazole rings are capped by hydrogen atoms at the C-y carbon atom. These carbon atoms (three total) are then bonded to their respective C-/3 carbon atoms through molecular mechanical bonds (using a standard AMBER carbon-carbon parameter set) between the two carbon atoms. Finally, the nonbonded interactions between the capping atoms and the remainder of the protein molecule are not evaluated. Computational Protocol. Starting coordinates for HCAII and crystallographically located water molecules were taken from the experimental crystal structure of Liljas and cow o r k e r ~ .TIP3P4 ~~ water molecules were added, where possible, around the active site zinc atom to a distance of 15.0 8, using the EDIT module of AMBER 4.0.43 All residues were represented using the AMBER united atom model,4 except for the active site residues His-94, -96, and -119, which were represented as all atom residues. Minimization of all atoms was carried out for both forms of HCAII using an all MM model and a standard version of AMBER 4.0. Following removal of bad contacts by MM minimization, a second minimization was carried out using a modified version of AMBER 4.0, in which the residues His-94, -96, and -119 and the active site zinc atom and its fourth ligand (either water or hydroxide) were treated as QM atoms using the PM3 Hamiltonian (for a total of 30 QM atoms for the zinc-hydroxide form and 31 for the zincwater f ~ r m ) . ~ ~ The - ~ ’ junction between the QM and MM regions was made between C-P and C-y of the His residues. In order to preserve integral charge in the MM region, the partial charges of the P-carbons of the QM His residues and the hydrogens attached to these carbons were changed to -0.080 and 0.048, respectively. Following the second minimization, MD simulations were carried out on both forms (water and hydroxide) of HCAII. A 15 8, sphere was defined around the active site zinc atom, and only residues within this sphere as well as the cap water molecules were allowed to move during the MD simulations. Each MD simulation consisted of 100 ps. The temperature was raised from 0 to 300 K during the first 2 ps of simulation, and the temperature was maintained at 300 K by coupling to a constant temperature heat bath.48 A 10 8, nonbond cutoff distance was used, and the nonbond pair list was updated every 25 time steps. The SHAKE algorithm was used to constrain all bonds between pairs of MM atoms.49 A 1 fs time step was employed during the MD simulations, and coordinates were saved for analysis every 50 time steps.
Results and Discussion
Rms Deviations. In no previous reports have protein dynamics been investigated using a semiempirical MO/MM coupled potential. Therefore, prior to any analysis of protein dynamics, it must first be established that the gross structural properties of the simulated protein are in at least rough agreement with the experimentally determined three-dimensional
1.0 1
.-
0.8
0
Y
.c
>
36
I
I
0.6
Oa4
0
20
40
60
80
100
Time (ps) Figure 2. Rms deviations of all atoms (top curve) and the backbone atoms (bottom curve).
structure. The easiest method by which to gauge the level of agreement between simulated and experimental structures is by examination of the root-mean-square (rms) deviation between the atomic positions in the experimental and simulated structures. Given the fact that large portions of the protein are frozen in place in our computational protocol (see Computational Procedures), it would be quite difficult to envision significant alterations in protein structure arising during the course of our simulations. It is nevertheless reassuring to observe that the rms deviation between the experimental and simulated structures for all of the mobile atoms in the simulation of the zinchydroxide form of HCAII is in the range 0.8-0.9 A, while the rms for only the backbone atoms is 0.5-0.6 8, (see Figure 2). Similar results were obtained from simulation of the zinc-water form of the enzyme (results not shown). These results are wholly consistent with the results one might expect to obtain by using an all MM force field (using a fixed “Belly” region), indicating that our combined QM/MM method has not introduced any severe perturbations into the computational model. Hydrogen Bond Interactions. The network of hydrogen bond interactions (ZnOH(H2)- - -Thr-199- - -Glu-106) which is proposed to orient the oxygen lone pairs favorably for catalysis is considered crucial for enzymatic reactivity (see Scheme 2).40.41 Accordingly, we feel it is of the utmost importance that our computational model retain this particular facet of the active site geometry. Upon minimization, both the zinc-water and zinc-hydroxide forms of the enzyme retain this group of hydrogen-bonded interactions. However, of greater import is the ability of this method to retain these interactions during the course of molecular dynamics simulations. Preliminary investigations using short time scale MD simulations (520 ps) showed that the hydrogen bond interactions in the zinchydroxide form of the enzyme were going to be problematic. Specifically,the stability of the ZnOH to Thr-199 hydrogen bond appeared to be significantly underestimated, as this hydrogen bond was observed to fragment with ease. Gao and co-workers have previously noted the difficulties arising from hydrogen bond interactions which traverse the QM-MM interface in which the hydrogen bond donor is a QM residue (as in the present instance).29 Accordingly, we decided upon adding a supplemental hydrogen bond potential between the protons of the fourth zinc ligand and the 0 - y oxygen of Thr-199. Keeping with the spirit of the AMBER force a 10-12 potential energy form was utilized for this supplemental interaction. The parameters for the 10-12 potential were determined by examination of a small model system ((NH3)3Zn-OH- - -H?O) incorporating the essential features of the enzymatic system.50 The 10-12 parameters were varied in an effort to obtain a potential energy surface for the hydrogen bond interaction
QM-MD Simulations of Human Carbonic Anhydrase I1
a
-
J. Phys. Chem., Vol. 99, No. 28, 1995 11269
6
'5 e
4
1
'
0
'
I
20
.
,
.
40
,
60
.
!
80
.
,
I
100
Time (ps)
5
! 0
'
I
20
'
I
40
1
'
60
.
I
'
I
80
100
BO
I00
Time (ps)
1 1 3
I
2
0
20
40
60
BO
100
Time (ps)
0
20
40
60
Time (ps)
Figure 3. (a) Fluctuations of the ZnOH- - -Thr-199 0 - y hydrogen bond. (b) Fluctuations of the Thr-199 HO-y- - -Glu- 106 0 - 6 1 hydrogen bond in the zinc-hydroxide form of HCAII. See Scheme 2 for atom labeling.
Figure 4. (a) Fluctuations of the ZnOH2- - -Thr-199 0 - y hydrogen hydrogen bond. (b) Fluctuations of the Thr-199 HO-y- - -Glu-106 0-€1 bond in the zinc-water form of HCAII. See Scheme 2 for atom labeling.
similar to that obtained for this same system from either ab initio calculations (using the 6-31G* basis set") or all MM calculations using the charge model of Hoops et aL3* In this effort, we attempted to strike a balance between retaining this crucial hydrogen bond interaction in the enzymatic system (suggesting the use of a deep, clearly defined 10-12 potential energy well) and avoiding disruption of other facets of the active site geometry, such as the Thr-199 to Glu-106 hydrogen bond (suggesting the use of a shallow ill-defined 10- 12 potential). Upon addition of this supplemental hydrogen bonding term, the active site geometry is fairly well preserved during MD simulation of the zinc-hydroxide form of the enzyme. As shown in Figure 3a, the hydrogen bond interaction between ZnOH and Thr-199 fluctuated between ca. 1.8 and 2.5 A, with occasional excursions to distances greater than 3.0 8, (average distance 2.10 A). Fortunately, while causing retention of the ZnOH to Thr-199 hydrogen bond interaction, the additional 1012 interaction term does not severely disrupt the Thr-199 to Glu-106 interaction. As seen in Figure 3b, this hydrogen bond possesses an average value of 2.20 A, with occasional elongation in excess of 3.0 A. However, in many of the instances in which this elongation occurs, a water molecule is observed to have inserted itself between Thr-199 and Glu-106, acting as a hydrogen bond bridge. Thus, while the hydrogen bond interaction between Thr-199 and Glu-106 may not be consistently present in the strictest sense, the geometric results of this interaction are well preserved throughout the course of the simulation. In the zinc-water form of the enzyme, both of the hydrogen bond interactions are significantly more stable than in the zinchydroxide case. In this instance, both hydrogen bond interac-
tions show little or no fluctuation throughout the entire MD simulation (see Figure 4a,b). Instead, the hydrogen-bonding residues are clearly locked into favorable geometries and show no inclination to vary from these positions. Average values for the Zn-OH? to Thr-199 (Figure 4a) and Thr-199 to Glu106 (Figure 4b) hydrogen bonds are 1.85 and 1.95 A, respectively. The greater stability of the hydrogen bond network in the aquo form may be attributed to the significantly lower pKa of the water ligand relative to the hydroxide. In fact, the system of hydrogen bonds was even well preserved in short test runs for the zinc-water form in the absence of the additional hydrogen bond interaction. However, in order to be consistent, this term was included in all calculations reported herein. Bonded Interactions. As a further check on the validity of the QM/MM method, we have also examined the stability of the bonded interactions. Unlike an all MM force field calculat i ~ n the , ~ present ~ simulations include no constraints on the bonds between pairs of QM atoms. Because of this freedom, the potential exists for alterations of bonding arrangements around the active site zinc atom. Although such processes would, of course, be quite energetically unfavorable, we must establish that the QMMM c o u p l e d potential provides a faithful rendering of bonded interactions about the zinc atom. An examination of instantaneous structures along the MD trajectory reveals that the bonding arrangements in the active site are well preserved throughout the MD simulations (for example, see the stereopictures given in Figure 5a.c). Further inspection reveals that the coupled potential method is also able to predict some of the important features of metal-ligand bonding in the active site. As shown in Figure 5, there are distinct differences in zinc-oxygen interactions in the zinc-
Hartsough and Merz
11270 J. Phys. Chem., Vol. 99, No. 28, 1995 a
‘a_ *U g *
2.5
-
2.3
-
2.1
5 i .9
1.7 0
20
40
60
80
I00
Time (ps) C
d
1.9
1
I
4
I
U
I .7 0
20
40
60
80
100
Time (ps)
Figure 5. Stereoplots of the superimposition of the experimental (labeled) and calculated (unlabeled) structures of (a) the Zn-OH and (c) ZnOH? forms of HCAII. The calculated structures are taken from a single representative snapshot along the 100 ps MD trajectory. Fluctuations of the (b) Zn-OH and (d) Zn-OH? bonds. See Scheme 2 for atom labeling.
water (see Figure 5c,d) and zinc-hydroxide (see Figure 5a,b) forms of the enzyme. The zinc-oxygen bond length is substantially shorter in the zinc-hydroxide form than in the zinc-water form (average length 1.9 vs 2.2 A). The high resolution (1.54 A) structure of HCAII gives the Zn-OH
distance as 2.05 A,51 while a lower resolution (2.1 A) structure gives a distance of 2.4 A*: both of which are longer than the average value observed herein. For the zinc-water form we have a 2.3 8, resolution structure obtained at pH 5.7 which gives a Zn-OH? distance of 2.7 A.53 However, this elongated
QM-MD Simulations of Human Carbonic Anhydrase I1
a
-
:
2.5 1
2.3
-I
1
I
I
2.1
*
a
‘I
I
I .7
0
20
40
60
EO
loo
Time (ps)
Figure 6. Fluctuations of the (a) Zn-N-BI(His-119) bond and the (b) Zn-N-€2 (His-96) bond in the zinc-hydroxide form of HCAII. See Scheme 2 for atom labeling.
distance could be due to the presence of an equilibrium mixture of water, chloride, and azide ion bound in the fourth coordination site.53 Nonetheless, our distance is again shorter than observed experimentally. However, given the resolution of the crystal structures it is difficult to make a definitive conclusion regarding the magnitude of the calculated errors. From our simulations, we also find that dynamic fluctuations of the zinc-oxygen bond’ appear to be significantly damped out in the zinc-hydroxide form relative to the zinc-water. Both of these observations are in good accord with expectations and e ~ p e r i m e n t . ~ In ’.~~ the zinc-water form the electron pair used for formation of the metal-oxygen bond is much more tightly held by the oxygen than it is in the anionic hydroxide ligand. Accordingly, a weaker bond is produced, as reflected by the longer interaction distance and weaker force constant for bond stretching. Variability is also observed in the linkages connecting the three histidine ligands to the active site zinc atom. In a conventional force field representation, all three histidine ligands would be constrained to a single zinc-nitrogen distance. Conversely, in the coupled potential, the zinc-nitrogen distance is allowed to vary in response to the molecular environment. As a result of this freedom, significant differences are observed in the zinc-nitrogen bond lengths (see Figure 6a,b). In particular, the nitrogen-metal bond between His-96 (Figure 6b) and zinc was observed to be elongated, on average, relative to those between zinc and His-94 (see Figure 6a) and His-119 (average bond lengths 2.10,2.00, and 2.00 A, respectively). The highest resolution structure of HCAII (1.54 A) indicates that the Zn-N bonds for His-94 and His- 119 are 2.1 8, , while this distance is 2.12 8, for H i ~ - 9 6 . ~Hence, ’ the best experimental information shows that the Zn-His-96 bond is elongated, but
J. Phys. Chem., Vol. 99, No. 28, 1995 11271 to a lesser extent than that observed herein. We note, however, that lower resolution structures5*show significant more variation in the Zn-N bond distances that seen in our simulations or the the high-resolution5‘ crystal structure. Examination of instantaneous structures along the MD trajectory reveals that a hydrogen bond is formed during the course of the simulations between the proton attached to the N-€2 nitrogen of His-96 and Glu-106, possibly pulling the His ligand away from the zinc atom. Arguing against this interpretation, however, is the existence of a hydrogen bond between Glu- 116 and His- 119 and the lack of a similar zinc-nitrogen bond elongation in this instance. Another possible explanation is the increased freedom of motion afforded His-96. This zinc ligating residue is much more exposed to the active site cleft and accordingly possesses a greater ability to move about relative to the other more tightly packed histidine residues. In any case, in addition to possessing an elongated zinc-nitrogen bond, the interaction between His96 and zinc also shows much greater dynamic fluctuation than those between zinc and the other two histidine residues. Atomic Charges. Having established that the QM/MM coupled potential produces no significant distortions of the active site geometry, we may now turn our attention to an area in which the coupled potential will excel in comparison to MM methodologies. Specifically, we wish to examine the alterations in electronic structure which occur during the course of the MD simulations. These changes may be most easily evaluated by examination of electrostatic potential (ESP) derived charges obtained along the MD trajectory. Gao has previously examined the polarization effects during the course of a QM/MM ~imulation.?~ However, that work was limited to small organic solutes in aqueous solution and did not consider systems such as metalloenzymes, which might be expected to show greater charge variation. In previous simulations of zinc metalloenzymes, one of the most important decisions made was how to represent the electrostatic properties of the metal atom and its ligand^.^*,^^ In some the earliest simulations of metalloenzymes, the entire problem was obviated by simply neglecting the metal atom charge or by simulating the apo-en~yme.~? A step beyond this approach was to assign the metal atom a charge equal to the formal charge (in the case of HCAKI a value of f l or f 2 ) . 5 4 However, this simplistic approach can produce severe geometric distortions due to the unrealistic electrostatic interaction^.^' In a more recent effort, Hoops er al. assigned atomic charges to the active site zinc atom of HCAII and its ligands based upon ESP calculations on small model systems3? Using this method, substantially smaller charges were obtained for the zinc atom than would have been expected based on the formal charge (0.688 and 0.866 vs f l and +2). The remaining positive charge was dispersed over the ligands due to charge transfer interactions. This electrostatic model, in conjunction with explicit bond, angle, and dihedral terms between the metal ion and its ligands (Le., the “bonded” model), produced a reasonable representation of the enzyme active site, as judged by deviations from the experimental geometry after 12 ps of MD simulation. However, while an ESP-based “bonded model” has proven quite successful, it still fails to convey the “dynamic” geometric and charge fluctuations inherent in the dynamics of a metalloen~yme.~~ Using the QM/MM coupled potential, we are able to examine the flow of charge between the metal atom and its ligands. Shown in Figure 7a,b are the values of the ESP-derived atomic charge for the zinc atom in the zinc-hydroxide and zinc-water forms of HCAII, respectively, during the last 70 ps of the 100 ps MD simulation. As might have been expected, the ESP-
11272 J. Phys. Chem., Vol. 99, No. 28, 1995
a
0.6 1
-0.2
I
0.8 I
50
70 Time (ps)
1
0.4
a
30
b
Hartsough and Merz
-0.2 4 30
90
1
I 70
50
90
Time (ps)
,
1
Figure 8. Fluctuations of the charge on the N-6, atom in His-1 19.
See Scheme 2 for atom labeling.
0.0
30
50
70
90
Time (ps) Figure 7. Fluctuations of the charge on the Zn atom in (a) the zinchydroxide and (b) the zinc-water forms of HCAII.
derived atomic charge of zinc shows rather substantial variation during both MD trajectories. The variation of these charges could be ascribed to the use of an ESP-based procedure, but from previous workss we have found that the error in the calculated charges arising from surface and orientation dependencies is -0.02 charge unit for buried atoms (like zinc here) and -0.005 charge unit for more exposed atoms. Hence, the observed charge variations are significant and arise due to various electronic and geometric effects. The average values obtained (0.25 and 0.60 charge units) here using the PM3 Hamiltonian are somewhat lower than those found by Hoops et al. based on a similar model system using the MNDO H a m i l t ~ n i a n .It~ ~is also immediately clear from these results that a static charge model is inadequate for the modeling of systems with large charge such as metalloenzymes. Of course, as the metal atom charge varies during the course of the MD simulation, the charge on the ligands must also vary in response. Such behavior is in fact observed during the MD trajectories. For example, the charge of the N-61 nitrogen atom of His- 119 (the nitrogen atom bound to the zinc) varies quite dramatically during the course of the MD simulation, covering a range of ca. 0.5 charge unit (see Figure 8). Similar variation may also be observed in the charges of atoms not directly bound to the metal atom. The atomic charge of the C-ci carbon of His-94, for example, also shows substantial variation during the MD simulation (see Figure 9). However, not all atoms show variations of such magnitude in their atomic charge. It is altogether reasonable that only those atoms involved in bonding to the metal center, either directly or through resonance, should show significant charge fluctuations. This is partly due to the weak metal-ligand bonds which can undergo more substantial bond distance fluctuations than can a carbon-carbon single bond for example. Furthermore,
30
so
70
90
Time ,(ps) Figure 9. Fluctuations of the charge on the C-EIatom in His-94. See Scheme 2 for atom labeling.
c
2
0.6
V u
-
*g 0.4-
1
.-
B L W
0.2
-
0.0 4 30
50
70
90
Time (ps) Figure 10. Fluctuations of the charge on the hydrogen atom bound to C-62 atom in His- 119. See Scheme 2 for atom labeling.
one would anticipate little or no charge variation for relatively nonpolarizable atoms which are not in direct (or resonance) contact with the metal center. In accord with these expectations, the protons attached to the heavy atoms of the His residues ligating the metal atom display only minor variations in atomic charge (see for example Figure 10). Thus, while molecular models which incorporate charge flux are of paramount importance for polarizable systems such as metal atoms and their ligands, these results suggest that such effects will be of relatively lesser importance for many amino acid side chains, which should be quite nonpolarizable in comparison. A second consideration in examining dynamic charge variation is the relative importance of electrostatic perturbations of the electron density and charge fluctuations arising from structural changes. This issue can most easily be addressed by
QM-MD Simulations of Human Carbonic Anhydrase I1
J. Phys. Chem., Vol. 99, No. 28, 1995 11273
a
f d
;
-o.2
0.4
V
u
*g a
0.2
.L
0.0
a
I
a VJ
Y
w -0.2
I 30
I 50
70
90
30
50
70
90
Time (ps)
Time (ps)
I
b
0.00
-0.04
V
0.00
! 30
I 50
70
90
Time (ps)
1
-0.16
-0.20
1 50
70
90
30
Time (ps)
Figure 11. (a) Fluctuations of the charge on the Zn atom in the zinchydroxide form of HCAII. The dotted line includes MM interactions, while the solid is in the absence of MM interactions. (b) Difference in the charges between the case including MM interactions at the neglecting them.
Figure 12. (a) Fluctuations of the charge on the Zn-oxygen atom in the zinc-hydroxide form of HCAII. The dotted line includes MM interactions, while the solid is in the absence of MM interactions. (b) Difference in the charges between the cases sincluding MM interactions and neglecting them.
recalculating the atomic charges of the QM atoms along the MD trajectory without the electrostatic influence of the enzymatic environment. The charge variations which are observed in the system not affected by the enzymatic environment can thus be regarded as the level of charge polarization due to geometric changes. The remaining difference in atomic charge between these charges and those obtained in the enzyme environment corresponds to polarization due to the electrostatic surroundings. Such an analysis was carried out for the zinc-hydroxide form of HCAII. In general, the charges obtained for the QM atoms in the absence of the MM environment track quite closely with those obtained in the presence of the enzyme (see Figures 11 and 12). In Figure 1la, the charge fluctuations in the presence of the MM interactions (dotted line) and without the MM interactions (solid line) are given for the zinc atom. In this case the charge of the zinc ion becomes more positive in the presence of the MM environment by -0.1 charge unit on average (see Figure 1lb). We can also plot the two charge sets against one another to determine whether there is a constant shift between the two charge sets. From this analysis we find that there is a reasonable correlation between the two charge sets (correlation coefficient = 0.95), which suggests in this case that polarization effects can be incorporated in an average way through the scaling of atomic charges. In Figure 12a,b we give the same analysis for the zinc bound oxygen, and in this case we find that the charge becomes more negative in the presence of the MM environment. The shift in charge is again on the order of -0.1 charge unit (see Figure 12b). However, the fluctuation in the charge is significantly less than in the case of
the zinc ion (see Figure 1la), and this is presumably due to the differences in the electronegativities of the atoms and the small fluctuations in the Zn-0 bond (see Figure 6a) as opposed to the Zn-N bonds (see Figure 7). Finally, by plotting the two charge sets relative to one another, we find that the correlation is quite low (correlation coefficient = 0.66) in the case of the oxygen atom, suggesting that polarization effects cannot be accurately included via a scaling scheme in this case. It is readily apparent that the level of charge polarization induced by the enzyme environment is significantly less than the charge variation observed during the MD trajectory for the zinc atom (and, for example, the zinc bound nitrogen atoms), while for the oxygen atom in the zinc-hydroxide form of the enzyme polarization appears to play the key role. The dichotomy appears to arise due to the differences in the bond distance fluctuations of the Zn-N bonds (relatively flexible) and the Zn-OH bond (relatively inflexible). Our results suggest that when a bond distance fluctuates considerably polarization plays a smaller relative role in the total charge variation. On the other hand, when the bond fluctuates less polarization becomes the dominant mechanism by w h i c h c h a r g e redistribution is accomplished. In particular, these results suggest that geometric changes (during the course of an MD simulation) play a more significant role in causing charge reorganization in metal systems than electron redistribution due to polarization effects, while in organic molecules the latter effect might be expected to be relatively more important due to reduced fluctuations in fully covalent bonds. However, while polarization has enjoyed significant attention, it should be stressed that charge redistribution due to polarization effects is only part of
11274 J. Phys. Chem., Vol. 99, No. 28, I995
a more complicated picture with regard to charge reorganization in complex molecular systems.56
Conclusions The brief analysis presented herein has demonstrated the utility of the QM/MM coupled potential method for the study of metalloenzymes. This coupling of QM and MM methodologies has been shown to provide a reasonable depiction of active site geometry and dynamics. In addition to allowing simulation of protein dynamics, this method also allows the study of quantum mechanical charge transfer interactions. Unlike purely classical MD simulations, the QM/MM method is able to model the dynamic flow of electron density as a result of protein motions. These changes in electronic structure, due to environmental influences and changes in geometry, can be expected to have a significant impact on energetics and forces in the active site region of metalloenzymes. Accordingly, we feel that this method holds great promise for the study of chemical reactions in an enzymatic environment.
Acknowledgment. We thank the NIH for supporting this research through Grant GM44974. We also thank the Pittsburgh Supercomputer Center and the Cornel1 Theory Center for generous allocations of supercomputer time. References and Notes ( I ) McCammon. J. A,: Harvey, S . C. Dvnamics of Proteins and Nucleic Acids: Cambridge University Press: New York. 1987. (2) Brooks, C. L.. 111: Karplus. M.; Pettit, B. M. Proteins: A Theoretical P erspective of Dynamics, Strucrure, and Thermodynamics: John Wiley and Sons: New York. 1988: Vol. LXXI. (3) Warshel. A. Computer Modelling of Chemical Reactions in Enpmes and Solutions: John Wiley & Sons: New York, 1991: p 236. (4) Weiner, S . J.: Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio. C.: Alagona. G.: Profeta. S.: Weiner. P. A New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J. Am. Chem. SOC. 1984. 106. 765-784. ( 5 ) Weiner. S . J.; Kollman. P. A.: Nguyen. D. T.: Case, D. A. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. J. Commit. Chem. 1986. 7. 230-252. (6) Brooks. B. R.: Bruccoleri, R. E.: Olafson. B. D.: States, D. J.: Swaminathan. S.: Karplus. M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Compur. Chem. 1983. 4. 187-217. (7) Dang. L. X.: Rice. J. E.: Caldwell, J.: Kollman, P. A. Ion Solvation in Polarizable Water: Molecular Dynamics Simulations. J. Am. Chem. Soc. 1991. 113. 2481-2486. (8) Sprik. M.: Klein, M. L. A Polarizable Model for Water Using Distributed Charge Sites. J. Chem. Phvs. 1988, 89. 7556-7560. (9) Jacob. 0.:Cardenas. R.; Tapia. 0. An ab initio Study of Transition Structures and Associated Products in [ZnOHCO?]+, [ZnHCO3H?O]+. [Zn(NH3)3HCO?]+Hypersurfaces. On the Role of Zinc in the Catalytic Mechanism of Carbonic Anhydrase. J. Am. Chem. SOC. 1990. 112. 86928705. ( 10) Jacob, 0.;Tapia, 0. On the Possible Roles of the Zn Cation in the Carbon Dioxide Hydration Reaction: An Ab Initio RHF-SCF MO Study of Transition Structure on the [OH? CO?] [ZNOH?CO?]?+and [Zn(NH3)3OH?CO?]?+ Reactive Hypersurfaces. Inr. J. Quantum Chem. 1992. 42, 1271-1289. ( 1 1) Zheng, Y. J.: Merz. Jr.. K. M. Mechanism of the Human Carbonic Anhydrase I1 Catalyzed Hydration of Carbon Dioxide. J. Am. Chem. Soc. 1992, 114. 10498-10507. (12) Merz. K.M.. Jr.: Hoffmann. R.: Dewar. M. J. S . Mode of Action of Carbonic Anhydrase. J. Am. Chem. SOC.1989, 111. 5636-5649. ( 13) Tapia, 0.:Johannin. G. An Inhomogeneous Self-Consistent Reaction Field Theory of Protein Core Effects. Towards a Quantum Scheme for Describing Enzyme Reactions. J. Chem. Phys. 1981, 75, 3624-3635. (14) Warshel. A.: Levitt. M. Theoretical Studies of Enzymatic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103. 227-249. (15) Tapia. 0.: Lamborelle. C.; Johannin. G. Towards a QuantumChemical Representation of Enzyme Activity. A SCRF PCE CNDOR Study of the LADH Proton Relay System. Chem. Phys. Lprr. 1980. 72. 334-341. (16) Field, M. J.: Bash, P. A.: Karplus. M. A Combined Quantum Mechanical and Molecular Mechanical Potential for Molecular Dynamics Simulations. J. Comput. Chem. 1990. 11. 700-733.
+
Hartsough and Merz (17) Luzhkov. V.; Warshel, A. Microscopic Models for Quantum Mechanical Calculations of Chemical Processes in Solutions: LD/AMPAC and SCAASIAMPAC Calculations of Solvation Energies. J. Compur. Chem. 1992. 13. 199-213. (18) Gao. J. Absolute Free Energies of Solvation from Monte Carlo Simulations Using Combined Quantum and Molecular Mechanical Potentials. J. Phys. Chem. 1992, 96. 537-540. (19) Singh. U. C.: Kollman, P. A. A Combined Ab initio Quantum Mechanical and Molecular Mechanical Method for Carrying out Simulations on Complex Molecular Systems: Applications to the CH3C1 C1Exchange Reaction and Gas Phase Protonation of Polyethers. J. Compur. Chem. 1986. 7 , 718-730. (20) Aqvist. J.: Warshel, A. Simulation of Enzyme Reactions Using Valence Bond Force Fields and Other Hybrid Quantum/ Classical Approaches; Chem. Rev. 1993. 93. 2523. (21) Aqvist. J.: Warshel. A. Computer Simulation of the Initial Proton Transfer Step in Human Carbonic Anhydrase I. J. Mol. Biol. 1992. 224. 7-14. (22) Stanton, R. V.: Hartsough, D. S.; Merz. Jr.. K. M. Calculation of Solvation Free Energies Using a Density FunctionaVMolecular Dynamics Coupled Potential. J. Phys. Chem. 1993. 97. 11868-11870. (23) Stanton. R. V.: Hartsough. D. S.: Merz, Jr., K. M. An Examination of a Density FunctionaVMolecular Mechanical Coupled Potential. J. Compur. Chem. 1995. 16, 113-128. (24) Stanton. R. V.; Merz. Jr.. K. M. Unpublished results. (25) Bash. P. A.: Field. M. J.: Davenport. R. C.: Petsko. G. A,; Ringe, D.: Karplus. M. Computer Simulation and Analysis of the Reaction Pathway of Triosephosphate Isomerase. Biochemisrp 1991, 30. 5826-5832. (26) Gao, J.: Pavelites. J. J. Aqueous Basicity of the Carboxylate Lone Pairs and the C - 0 Barrier in Acetic Acid: A Combined Quantum and Statistical Mechanical Study. J. Am. Chem. SOC.1992, 114. 1912-1914. (27) Gao. J.: Xia, X. A Priori Evaluation of Aqueous Polarization Effects Through Monte-Carlo QM-MM Simulations. Science IWashingron. D.C.) 1992. 258. 63 1-635. (28) Gao. J. Absolute Free Energies of Solvation from Monte Carlo Simulations Using Combined Quantum and Molecular Mechanical Potentials. J. Phys. Chem. 1992. 96. 537-540. (29) Gao. J. Comparison of the Hybrid AMliTIP3P and the OPLS Functions Through Monte Carlo Simulations of Acetic Acid in Water. J. Phys. Chem. 1992, 96. 6432-6439. (30) Gao. J.; Xia. X. A Two-Dimensional Energy Surface for a Type I1 SN2 Reaction in Aqueous Solution. J. Am. Chem. SOC. 1993. 115. 96679675. (31) Gao. J. Solvent Effect on the Potential Surface of the Proton Transfer in [NHj - H - NH?]+. Inr. J. Quantum Chem.: Quantum Chem. Symp. 1993. 27. 491-499. (32) Hoops. S . C.: Anderson, K. W.: Merz. Jr., K. M. Force Field Design for Metalloproteins. J. Am. Chem. SOC.1991, 113. 8262-8270. (33) Vedani, A.: Huhta. D. W. A New Force Field for Modeling Metalloproteins. J. Am. Chem. SOC. 1990, 112. 4759-4767. (34) Vedani. A,: Huhta, D. W.: Jacober, S . P. Metal Coordination. H-bond Network Formation. and Protein-Solvent Interactions in Native and Complexed Human Carbonic Anhydrase I: A Molecular Mechanics Study. J. Am. Chem. SOC. 1989. 111. 4075-4081. (35) Silverman. D.: Vincent. S . H. Proton Transfer in the Catalytic Mechanism of Carbonic Anhydrase. CRC Crir. Rev. Biochem. 1983. 14. 207-255. (36) Silverman. D. N.: Lindskog. S . The Catalytic Mechanism of Carbonic Anhydrase: Implications of a Rate-Limiting Protolysis of Water. Acc. Chem. Res. 1988, 21. 30-36. (37) Lindskog, S . In Zinc Enzymes: Spiro, T. G.. Ed.: Wiley: New York. 1983: pp 77-121. (38) Tashian. R. E. The Carbonic Anhydrases: Widening Perspectives on Their Evolution. Expression and Function. BioEssays 1989, 10. 186192. (39) Eriksson. A. E.: Jones. A. T.; Liljas. A. Refined Structure of Human Carbonic Anhydrase I1 at 2.0 A Resolution. Proreins 1988. 4. 274-282. (40) Xue. Y.; Liljas, A,: Jonsson. B.-H.: Lindskog, S . Structural Analysis of the Zinc Hydroxide-Thr 199-Glu106 Hydrogen Bond Network in Human Carbonic Anhydrase 11. Proreins 1993, 17, 93- 106. (41) Merz. K. M., Jr. Insights into the Function of the Zinc HydroxideThr 199-Glu106 Hydrogen Bonding Network in Carbonic Anhydrase. J. Mol. Bioi. 1990. 214. 799-802. (42) Merz, K. M.. Jr. CO? Binding to Human Carbonic Anhydrase 11. J. Am. Chem. Soc. 1991. 113. 406-411. (43) Pearlman, D. A,: Case. D. A,: Caldwell. J. C.: Seibel. G. L.: Singh, U. C.: Weiner. P.: Kollman. P. A. In University of California. San Francisco: 199 1. (44) Jorgensen. W. L.: Chandrasekhar, J.: Madura, J.: Impey, R. W.: Klein. M. L. Comparison of Simple Potential Functions for the Simulation of Liquid Water. J. Chem. Phys. 1983. 79. 926. (45) Stewart. J. J. P. Optimization of Parameters for Semiempirical Methods I. Method. J. Compur. Chem. 1989. 10. 209-220. (46) Stewart. J. J. P. Optimization of Parameters for Semiempirical Methods 11. Applications. J. Compur. Chem. 1989. 10. 221-264.
+
QM-MD Simulations of Human Carbonic Anhydrase II (47) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods. III Extension of PM3 to Be, Mg, Zn, Ga, Ge, As, Se, Cd, In, Sn, Sb. Te. Hg, Tl,Pb and Bi. J. Comput. Chem. 1991, 12, 320-341. (48) Berendsen, H. J. C.; Potsma, J. P. M.: van Gunsteren, W. F.; DiNola, A. D.; Ha&, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. P h y . 1984, 81, 3684-3690. (49) van Gunsteren, W. F.; Berendsen, H. J. C. Algorithm for Macromolecular Dynamics and Constraint Dynamics. Mol.Phys. 1977.34, 13 11. (50) Hartsough, S. D.: Merz. Jr., K. M. Unpublished results. (51) H h s s o n . K.: Carlsson, M.; Svensson. L. A.; Liljas, A. The Structure of Native and Apo Carbonic Anhydrase 11, and some of its Anion Ligand Complexes. J. Mol.Bid. in press. (52) Alexander, R. S.; Nair, S. K.; Christianson, D. W. Engineering the Hydrophobic Pocket of Carbonic Anhydrase 11. Biochemist, 1991, 30. 11064- 11072.
J. Phys. Chem., Vol. 99, No. 28, 1995 11275 (53) Nair, S. K.: Christianson, D. W. Unexpected pH-Dependent Conformation of His-64, the Proton Shuttle of Carbonic Anhydrase 11. J. Am. Chem. SOC. 1991. 113, 9455-9458. (54) Merz, Jr., K. M.: Kollman, P. A. Free Energy Perturbation Simulation of the Inhibition of Themolysin: Prediction of the Free Energy of Binding of a New Inhibitor. J. Am. Chem. SOC. 1989, 1 11, 5649-5658. (55) Men, K. M., Jr. Analysis of a Large Database of Electrostatic Potential Derived Atomic Point Charges. J. Compur. Chem 1992, 13,749767. (56) Dinur, U. Charge Flux and Electrostatic Forces in Planar Molecules. J. P h y . Chem. 1991, 95, 6201-6211. (57) Dinur, U. “Flexible” Water Molecules in External Electrostatic Potentials. J. Phys. Chem. 1990, 94, 5669-5671. Jp950374S
