Catalytic Mechanism of Angiotensin-Converting Enzyme and Effects of

May 15, 2013 - Elena V. Kugaevskaya , Alexander V. Veselovsky , Maria I. Indeykina , Nina I. Solovyeva , Maria S. Zharkova , Igor A. Popov , Eugene N...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCB

Catalytic Mechanism of Angiotensin-Converting Enzyme and Effects of the Chloride Ion Chunchun Zhang,† Shanshan Wu,‡ and Dingguo Xu*,‡ †

Analytical&Testing Center and ‡MOE Key Laboratory of Green Chemistry and Technology, College of Chemistry, Sichuan University, Chengdu, Sichuan, 610064, P. R. China S Supporting Information *

ABSTRACT: The angiotensin-converting enzyme (ACE) exhibits critical functions in the conversion of angiotensin I to angiotensin II and the degradation of bradykinin and other vasoactive peptides. As a result, the ACE inhibition has become a promising approach in the treatment of hypertension, heart failure, and diabetic nephropathy. Extending our recent molecular dynamics simulation of the testis ACE in complex with a bona fide substrate molecule, hippuryl-histidyl-leucine, we presented here a detailed investigation of the hydrolytic process and possible influences of the chloride ion on the reaction using a combined quantum mechanical and molecule mechanical method. Similar to carboxypeptidase A and thermolysin, the promoted water mechanism is established for the catalysis of ACE. The E384 residue was found to have the dual function of a general base for activating the water nucleophile and a general acid for facilitating the cleavage of amide C−N bond. Consistent with experimental observations, the chloride ion at the second binding position is found to accelerate the reaction rate presumably due to the long-range electrostatic interactions but has little influence on the overall substrate binding characteristics. (CPA).8 In 2003, the first crystallographic structure of the ACE in complex with an inhibitor, lisinopril, was determined at 2.0 Å resolution.9 Later, some other X-ray structures of tACE in complex with inhibitors were obtained.10 Interestingly, apart from the conserved HEXXH zinc binding motif, there are a few structural similarities between ACE and CPA. More importantly, these crystallographic structures of ACE-inhibitor complexes can provide a firmer basis for further rational design. The wide usage of ACE inhibitors in the treatment of vascular diseases underscores its importance. Currently, much attention has been devoted to understand both the inhibition mechanism11,12 and the inhibitor development.13−17 Early stage inhibitors have little domain selectivity, thus having adverse side effects. With the X-ray structures solved for both the C domain9,10,18,19 and N domain,20 some successful domainselective inhibitors have been reported, for example, RXP407 to the N domain7 and RXPA380 to the C domain.18 Nevertheless, the lack of information concerning the native substrate binding pattern and the catalytic mechanism still restricts the rational design of more efficient inhibitors. Theoretical simulations represent a useful way to address these issues from the molecular level. Recently, Papakyriakou et al.21 reported a classical molecular dynamics (MD) study on the enzyme in complex with a long peptide substrate, the gonadotropinreleasing hormone (GnRH), in which they suggested that the

1. INTRODUCTION The angiotensin-converting enzyme (ACE) is a zinc-containing dipeptidase, which was discovered by Skeggs et al. in 1956.1 This enzyme exhibits important biological functions in the conversion of angiotensin I to angiotensin II,2 a potent vasoconstrictor. It also shows functions in the degradation of bradykinin, a potent vasodilator, and other vasoactive peptides.3 These two actions endow this enzyme with critical roles in blood-pressure regulation. As a result, the ACE has become a promising target in the treatment of hypertension, heart failure, and diabetic nephropathy. Two types of ACE have been identified, human somatic angiotensin-converting enzyme (sACE) and testis angiotensinconverting enzyme (tACE). It was suggested that tACE plays an important role in male fertility.4,5 The human sACE contains two catalytic domains, namely, the N domain and the C domain.6 They have ∼55% sequence similarity. Both domains possess the same function for converting angiotensin I to angiotensin II with similar efficiency. The C domain was found to be the dominant angiotensin-converting site in controlling blood pressure and cardiovascular functions based on the observation that the inhibition of the N domain by RXP407 has no effect on these functions.7 Testis ACE is almost identical to the C domain of sACE, except the first 36 residues. Remarkably, the FDA-approved drug, L-captopril, was designed for the orally active ACE inhibitor for treatment of human hypertension without any structural information of the enzyme but based on the assumption of a similar mechanism with another metallopeptidases, for example, carboxypeptidase A © XXXX American Chemical Society

Received: January 29, 2013 Revised: May 12, 2013

A

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Scheme 1. Proposed Promoted-Water Mechanism for the Hydrolysis of Hip-His-Leu Catalyzed by ACE

zinc cofactor is penta-coordinated with both scissile carbonyl oxygen atom of the substrate and active site water as its ligands. Subsequently, combined quantum mechanical and molecular mechanical (QM/MM) MD simulations of the ACE in complex with the hippruyl-histidyl-leucine (Hip-His-Leu), a synthetic analogue of angiotensin I, were reported by us.22 Our simulations showed a different binding pattern in which no direct interaction between the substrate and zinc cofactor can be observed. The zinc ion is tetra-coordinated with H383, H387, and E411 and one putative active site water molecule. This water molecule is hydrogen-bonded with a nearby E384 and stays at a near attack position. As shown in Scheme 1, a glutamic acid residue assisting water addition mechanism was speculated for the hydrolysis of Hip-His-Leu. Such a mechanism is similar to those in the hydrolysis of peptide catalyzed by CPA23,24 and thermolysin (TLN).25 Another interesting feature in ACE is the biological functions of two chloride ions. Two chloride ions have been identified outside the active site, which are proposed to be important in either the stabilization of the overall enzyme structure or facilitation of the catalytic reaction.26−29 This is unique in metallopeptidases. It should be pointed out that the functions of chloride ions are very complicated and still under much debate. Experimental evidences indicated that the activation of the hydrolysis rate by the chloride ions is highly depended on the nature of the substrate molecule.27,29,30 Indeed, about 13fold increasing of the hydrolysis rate for Hip-His-Leu was observed in the presence of the 300 mM NaCl compared with the absence of the chloride ion.30,31 The presence of chloride ion may also inhibit the hydrolysis reaction in some cases, for example, the hydrolysis of Hip-Ala-Pro, which is 80% inhibited with 100 mM NaCl.30 However, none of these observations provides microscopic insights into their actions. On the basis of the X-ray structures and shown in Figure 1, the chloride ion at the first binding position (Cl−(I) hereafter) is ∼21 Å away from the zinc ion. It is stabilized by a hydrogen-bond network formed by two positively charged residues of R186 and R489 and by a hydrophobic shell formed by W485 and W486 sidechain groups and the D507 backbone. The chloride ion at the second binding position (Cl−(II) hereafter) is ∼10 Å away from the zinc ion and is hydrogen-bonded with Y224, R522, and water molecule. Similarly, a hydrophobic shell formed by P407, P519, and I521 can be found surrounding Cl−(II). A generally accepted opinion is that Cl−(II) stays in the center position in the activation of the reaction based on some kinetic experiments.26,29 It was also suggested that the absence of Cl−(II) in the ACE homologue from Drosophila melanogaster, AnCE, does not seem to affect the binding of small molecule inhibitors, such as lisinopril and enalaprilat.32,33 This article will then be divided into two parts. We will first discuss the detailed hydrolytic processes of Hip-His-Leu

Figure 1. Illustration of the surrounding environment of two chloride ions and the active site of tACE using the crystal structure of the ACE in complex with lisinopril molecule (PDB code 1O86).9

catalyzed by ACE. Then, we will address the function of the chloride ion at the second binding position in the hydrolysis reaction because it has been suggested that the Cl−(II) plays major role in the acceleration of the catalytic rate.

2. COMPUTATIONAL DETAILS 2.1. QM/MM Method. It has been widely accepted that the hybrid quantum mechanical and molecular mechanical approach (QM/MM) is one of the most useful and efficient tools to simulate large systems like enzyme.34−36 The basic idea of QM/MM is based on the divide-and-conquer strategy. In general, the entire system is divided into two parts. The smaller region (QM region), which contains the residues involved in the bond forming and bond-breaking processes, will be treated using a high-level quantum mechanical method, while the surrounding MM region, which contains most of environment atoms, will be described classically using a force field. The covalent interface between the QM and MM regions in this work is approximated using the link atom algorithm.35 The accuracy of the simulation results often depends on the selection of the QM method in the calculation of electronic structure of the QM region. Although first-principle methods are ideal for the QM/MM treatment, their applications are often restricted by the low computational efficiency. In this work, to efficiently calculate the electronic structure of the QM region, a semiempirical approach, namely, the selfconsistent charge-density tight binding (SCC-DFTB) method,37−39 was applied throughout the simulations. This method is derived by a second-order expansion of the DFT total energy functional with respect to the charge density fluctuation around a given reference density. The parameters for the biological zinc B

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

ion were specially developed in the SCC-DFTB framework.40 The combined SCC-DFTB/MM method has successfully applied to investigate several enzyme systems, especially some zinc hydrolases, including carbonic anhydrase,41−43 carboxypeptidases,24,44 and metallo-β-lactamases.45−47 The atoms in MM region were expressed using CHARMM22 all atom force field.48 2.2. QM/MM Model. The model construction has been discussed in detail in our recent work.22 Here only a short description is given. Starting from the tACE-lisinopril complex structure (PDB code 1O869), we first removed the lisinopril molecule from the active site. The substrate of Hip-His-Leu was manually docked into the active site based on a similar backbone molecular structure between substrate and the inhibitor of lisinopril. The resulting system was solvated by a pre-equilibrated TIP3P49 water sphere with a 25 Å radius centered at the zinc ion, followed by 30 ps solvent MD with all protein and substrate molecules fixed. This process was repeated several times to ensure that the system was sufficiently solvated. Stochastic boundary condition50 was applied to reduce computational cost. A group-based switching scheme was applied for nonbonded interactions.51 All calculations reported here were carried out using CHARMM. As depicted in Figure 2, the QM region contains the entire substrate, the zinc ion, side chain groups of H383, E384, H387,

the peak region. The harmonic force constants were used in the simulation windows ranging from 100 to 150 kcal/mol·Å2. A total of 100 ps constrained MD simulation was performed in each window, with first 60 ps equilibration, and following 40 ps simulation to evaluate the distribution function along the reaction coordinate. The final potential of mean force (PMF) was obtained using the weighted histogram analysis method (WHAM).54 For the NA and E steps, a total of 18 and 21 windows were used, respectively. The integration step for MD simulations was 1.0 fs. In the PMF calculation, the SHAKE algorithm55 was applied to keep all covalent bonds involving hydrogen atoms except the hydrogen atoms of the water nucleophile. 2.3. Proton Affinities Calculations. Because the hydrolysis of the Hip-His-Leu catalyzed by ACE has been suggested to be dominated by several proton-transfer processes,22 it would be highly necessary to further carry out proton affinity (PA) calculations for some important molecule species involved in the reaction. In this work, we constructed five models as shown in Figure 3, for which the PAs were calculated

Figure 3. Computational models for proton affinities used in this work: (a) Ac-Glu-OCH3, (b) Ac-His-OCH3, (c) CH3NH−, (d) CH3NH2, and (e) Zn−OH−.

Figure 2. Atom definition and corresponding interactions between substrate and active site residues of tACE.

at both SCC-DFTB and B3LYP/6-311G++(3df, 2p) level of theory. The gas-phase PA of a compound B can be calculated using the following scheme,56 that is, the negative standard reaction enthalpy of protonation at 298.15 K:

and E411, and the zinc-bound water molecule. The total number of atoms in the QM system is 121. To independently examine the accuracy of the SCC-DFTB method, we carried out single-point calculations at the B3LYP/MM//SCC-DFTB/ MM level along the reaction coordinates. The standard 631G(d) basis set was applied for the C, H, O, and N atoms, and the LANL2DZ ECP basis set was used for zinc ion. These calculations were performed using a GAMESS-UK/CHARMM interface.52 As we have discussed in our previous work, the resulting system was finally subjected to 1 ns MD simulation. One of the snapshots was randomly selected from the MD trajectory in the following simulation of enzymatic reaction. Similar to the reaction coordinates definition employed in the hydrolysis of HPA catalyzed by CPA,24 the first step of the reaction is the (NA) by the zinc-bound water nucleophile at the scissile carbon center, and the second step is the elimination step (E) to cleave the amide C−N bond to fulfill the hydrolysis. The corresponding reaction coordinate for the NA step was defined as r1 = dC1···Ow, while for the E step, the reaction coordinate was set as r2 = dH1···Oε2(E384) − dH1···N2. To better depict the reaction, we need the free-energy profiles for each step. The umbrella sampling method53 was applied in both two models to enhance the sampling around

B + H+ → BH+

(1)

PA = −[ESCF(BH+) − ESCF(B) + (Evib(BH+) 5 − Evib(B))] + RT (2) 2 in which ESCF represents the electronic energy obtained from SCF calculations, Evib includes the zero-point energy and temperature corrections to the vibrational enthalpy, and (5/2) RT denotes the contributions from translational energy of proton and the ΔPV term. For SCC-DFTB, the PA can be calculated as:

PA = −[ESCC − DFTB(BH+) − ESCC − DFTB(B) + Erep[ρ0H ] 5 RT (3) 2 In this work, we calculated the proton affinities using the value of −141.9 kcal/mol for Erep[ρH0 ].56 Similar to Zhou et al.,56 the vibration term Evib will not be considered in this work due to the comparative feature of the calculations. + (Evib(BH+) − Evib(B))] +

C

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

3. RESULTS AND DISCUSSION 3.1. Proton Affinities. Models used in this work are based on the possible reactive mechanism (shown in Scheme 1) proposed by us,22 for which the active site E384 serves as the proton acceptor and zinc-bound water serves as the proton donor in the first step of NA. In the second step, this proton will transfer to the new formed amine (R1-NH−) due to cleavage of C−N bond. There are several histidine residues in the active site, for example, H353 and H513 were shown to contribute the binding of the substrate of Hip-His-Leu molecule via hydrogen bonds. It would be also interesting to investigate the PA for the side chain of the histidine residue. Therefore, we constructed five models, as shown in Figure 3. For glutamic acid and histidine residue, we capped the N terminal using an acetyl group and the C terminal using a CH3O− group. For the zinc bound water, methyl imidazole is used to mimic the zinc ligand of histidine, and the propanoic acid is used to replace the glutamic acid. Calculated PAs are listed in Table 1 obtained by both SCCDFTB and density functional theory (DFT). Calculated PAs

neutral water molecule, while the glutamic acid (E384) should be deprotonated. Nevertheless, to fully understand the whole catalytic processes, more calculations are required to get the energy curves along the putative reaction path. 3.2. Model for the Native Enzyme. The structure of ACE in complex with the Hip-His-Leu has been investigated by us using the same combined QM/MM MD simulation. The detailed binding information can be found in our previous published work.22 Starting from a randomly selected structure from the MD trajectory, we performed reaction path calculations along the putative reaction coordinates using the adiabatic mapping approach (also sometimes referred to as reaction coordinate driving method; details can be seen in ref 59). The obtained minimum energy paths are further refined using conjugate peak refinement (CPR) approach.60 It has been suggested that CPR-corrected energy activation barrier height can be sufficiently close to experimental value, especially the proton transfer reactions.61 The reaction coordinate in CPR is then defined as N

λ=

Table 1. PA Values (kcal/mol) Calculated Using Both SCCDFTB and DFT (B3LYP/6-311++G(3df, 2p)) Approachesa,b

i=2

SCC-DFTB

without Evib

with Evib

Zn-OH− CH3-NH2 CH3-NH− Ac-Glu-OCH3 Ac-His-OCH3

283.4 219.6 431.3 359.4 256.6

275.5 219.5 412.4 346.8 252.7

267.4 214.5 401.3 337.9 244.4

3n (4)

where N is the order of the points along the CPR trajectory, x(i) is the coordinate of point i, and n is the number of atoms. Similar to Friedman et al.,61 λ was normalized throughout this work if not otherwise stated, so that 0 ≤ λ ≤ 1. The so-called promoted-water pathway has been established for both CPA24 and TLN25 for the hydrolysis of peptide substrates, in which the first step is the nucleophilic attack by a zinc-bound water molecule at the scissile carbonyl carbon assisted by a nearby glutamic acid residue. Once the tetrahedral intermediate is formed, this protonated Glu residue will serve as the general acid catalyst by transferring the proton to the amide nitrogen atom to facilitate the final cleavage of C−N bond. Such a mechanism was postulated according to the substrate binding mode. The corresponding mechanism is summarized in Scheme 1. Selected key geometric parameters for all five stationary structures are listed in Table 2, in which geometries for two transition states are obtained by CPR calculations. The energetic profiles are displayed in Figure 4A,B for the NA and E steps. Corresponding snapshots along the reaction coordinates are displayed in Figure 5. As the approaching of the nucleophile (Ow) to the scissile carbonyl carbon (C1) atom, the proton of water molecule was fully transferred to E384 Oε2 atom evidenced by the distance of dH1···Oε2(E384) shortening from 1.25 to 1.01 Å at the transition state. This indicates that E384 serves as the general base in the NA step. Only one transition state can be located for the first step of NA. The distance between the scissile carbonyl oxygen and the zinc cofactor is significantly shortened (4.60 Å at ES vs 2.33 Å at TS1), which could further suggest that zinc cofactor will serve as the oxyanion hole to stabilize the polarized carbonyl oxygen atom. Such a phenomenon can be found in many other enzymes belonging to zinc hydrolase subfamily.24,25,45 At the same time, the strong hydrogen bond provided by Y523 to the carbonyl oxygen atom also indicates that Y523 has the function of stabilizing the carbonyl oxygen atom to facilitate the hydrolysis. The NA step leads to the formation of a tetrahedral enzyme−intermediate complex (EI) characterized by the conversion of sp2 to sp3 hybridization of the scissile carbon center (C1), evidenced by the distance of C1−O2 elongated from 1.23 to 1.42 Å. In addition, to maintain the stability of the

DFT models

∑ {|x (⃗ i) − x (⃗ i − 1)|}/

exp.67 214.9 403.2

Zn−OH− represents a model as shown in Figure 3e, in which methyl imidazole is used to mimic the zinc ligand of histidine and the propanoic acid is used to replace the glutamic acid. bAc = acetylprotecting group. a

using DFT method can compare very well with experimental values, as shown in Table 1 and previous calculations.57 The SCC-DFTB seems to be an insufficient level of theory to get accurate PA values for model systems. However, we should emphasize that the overall trend for the calculated models are generally consistent between SCC-DFTB and DFT. More importantly, the computational results could indicate that the glutamic acid has stronger PA than the zinc bound water, which might further confirm that the Glu residue could serve as the general base to activate the zinc-bound water. Moreover, the PAs without contributions of vibrational terms for the model of CH3NH− are calculated to be 431.3 and 412.4 kcal/mol at SCC-DFTB and DFT level of theory, respectively, which is much larger than glutamic acid. This is also an additional evidence that supports our reactive mechanism proposal for the proton migration in the second step. However, it should be noted that for the initial protonation status of those residues in the active site, the values of pKa for those titratable groups might be more direct. Because of nearly the same active site environment among CPA, TLN, and ACE, we could use the pKa values in CPA or TLN for some discussions. As pointed by Lipscomb and Strater,58 for CPA, the pKa for the general base of E270 is assigned to be 4.3, while pKa for the zinc bound water is ∼6.3. For TLN, the pKa of the general base of E143 is assigned to be 5.3. This could tell us that when the substrate is binding in the active site, the water ligand to zinc ion is a D

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

hydrolysis of the hippuryl-L-Phe (HPA) molecule catalyzed by CPA.24 The overall reaction can be grouped into the general acid/general base (GAGB) or the so-called promoted-water mechanism. To ascertain the reliability of the selected reaction coordinates, we further performed B3LYP/MM//SCCDFTB/MM single-point calculations along the reaction path. The corresponding reaction curves for two steps are shown in Figure 4A,B. Clearly, both reaction pathways obtained at two different levels of theory agree with each other qualitatively. Moreover, both methods confirm that the first NA step is more likely the rate-determined step for the whole reaction. It has been shown that the CPR could largely reproduce the experimental barrier height. We have also included CPR refined minimum energy pathway in Figure 4C. It can be that the first step clearly dominates the whole reaction with the barrier height of 24.3 kcal/mol, and the second barrier height is ∼5.2 kcal/mol. The reactive energetic profiles are largely consistent between CPR and minimum-energy pathways calculated at the SCC-DFTB/MM level of theory. Although our reaction path calculations presented above support the promoted water mechanism in the hydrolysis of Hip-His-Leu, it would be better to obtain the free-energy profiles to account for the entropic contributions. The calculated potentials of mean force (PMFs) for both NA and E steps are displayed in Figure 5. As we can see, the overall freeenergy profiles are consistent with our minimum energy path calculations. The first step of NA has the activation energy barrier height of 17.4 ± 0.4 kcal/mol with respect to the sampling time, in which the data analysis in PMF calculations is carried out from the trajectories between 60 to 100 ps. However, only the PMF curve obtained using 100 ps dynamics simulation with a barrier height of 17.8 kcal/mol for the NA step is plotted in Figure 6. It can be comparable to the experimental value63 of 14.9 kcal/mol estimated by observation rate constant of kcat = 99 s−1 using the transition-state theory with the assumption of a unit transmission coefficient. In addition, a relatively shallow free-energy well can be found related to the tetrahedral intermediate; ∼4.9 ± 0.8 kcal/mol needs to be overcome to finally cleave the C−N bond. Overall, the rate-limiting step is still the first step of nucleophilic attack. If we go back to check the reaction mechanism of CPA and TLN, then it can be easily found that all three types of zincdependent hydrolases share an identical reactive mechanism, with an active site glutamic acid residue to serve as both the

Table 2. Selected Geometric Parameters from Stationary Points along the Reaction Coordinatesa QM/MM reaction path distance (Å) angle (deg.)

ES

TS1

EI

TS2

EP

Zn···Ow C1···Ow H1···Ow C1···N2 H1···Oε2 (E384) C1···O2 H1···N2 H2···Oε1 (E384) H2···Ow Zn···Nε2 (H387) Zn···Nε2 (H383) Zn···Oε1 (E411) Zn···O2 O2···Hε2 (H513) O2···HH(Y523) O3···Hε2 (H353) H4···O(A354) O4···HH(Y520) O5···H21(Q281) O4···Hζ1(K511)

2.01 2.82 1.17 1.37 1.25 1.23

2.21 1.85 1.72 1.39 1.01 1.33

2.88 1.42

2.90 1.26

1.99 1.98 2.02 4.60 1.86 2.34 1.72 1.92 1.64 2.00 1.78

2.02 2.02 2.13 2.33 3.89 1.66 1.79 2.44 1.67 1.89 1.76

2.89 1.48 1.80 1.45 1.00 1.42 2.22 2.91 1.01 2.01 1.99 2.06 2.07 4.10 1.69 1.83 2.48 1.61 1.76 1.80

1.79 1.22 1.32 1.30 1.77 1.03 1.99 1.99 2.03 2.07

3.15 2.58 1.31 1.02 1.00 1.77 1.98 1.97 2.04 2.06

1.68 1.86 2.11 1.62 1.98 1.83

1.81 1.89 1.62 2.01 1.79

a

Geometries for the transition states are obtained using conjugate peak refinement calculations along the minimum energy paths.

EI complex, the ligation state of the carbonyl oxygen atom is very important. After the formation of TI, the carbonyl oxygen atom changes its conformation and replaces the water molecule to become the fourth ligand to zinc ion (dO2−Zn = 2.07 Å). Although the hydrogen bond is initially formed between O2 and H513, it is lost due to the conformal change with the formation of EI. In addition, the C-terminal carboxylate group of the peptide molecule is stabilized by several hydrogen bonds formed with Q281, K511, and Y520, respectively. The second step involves the cleavage of the C1−N2 bond to complete the hydrolysis process. Just like E270 in CPA23,24 and E143 in TLN,25,62 E384 in ACE also plays the role of the general acid at the second step to transfer the proton to the amide nitrogen atom to cleave the C−N bond. Another hydrogen atom of the water nucleophile was also finally transferred to E384, which is the same as we observed for the

Figure 4. Minimum energy pathways (black triangle line) for the nucleophilic addition step (panel A) and elimination step (panel B) of the hydrolysis of Hip-His-Leu molecule calculated with the SCC-DFTB/MM model. Also plotted in the Figure are the B3LYP/MM single-point energies (blue dot line) along the SCC-DFTB/MM reaction path. CPR-refined energetic profile for the hydrolysis of Hip-His-Leu is given in panel C with the reaction coordinate (λ) based on the structural difference between points along the CPR trajectory. E

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 6. Potential of mean force for the hydrolysis of Hip-His-Leu catalyzed by ACE calculated using the combined SCC-DFTB/MM method.

activities to hydrolyze both the ester and peptide molecules. Such functions have not been observed for ACE. 3.2. Effect of Cl−(II). A. MD Simulation of Hip-His-LeutACE Complex without Cl−(II). It has been discussed extensively in the literature that the chloride ions in ACE might play a key role in maintaining the overall structural stability or in the enzymatic activity. In our previous MD simulations of the ACE inhibitor and substrate binding, due to the relatively short simulation time, it was not possible to investigate the influence of the chloride ion at the first binding position (∼20 Å away from the zinc ion) to the structure. In this work, we will focus on the effects of Cl−(II), which is ∼10 Å away from the zinc ion according to the X-ray structure9 and our recent MD simulation.22 The QM/MM MD simulation was performed to investigate its influence on the substrate binding and the enzymatic action. To quantitatively understand the contributions of Cl−(II) to the enzymatic reaction, we designed two model systems. The first model is simply to remove the chloride ion (Cl−(II)) from its binding site (Model I), while the second model is to keep the chloride ions in their binding positions but manually scaling the charge of Cl−(II) to zero to examine the effects of the longrange electrostatic interaction (Model II). Here the QM/MM MD simulation was applied to the first model to check the dynamics of the substrate in the active site. A recent classical MD simulation has been reported by Papakyriakou et al. based on the binding pattern of the ACE C-domain complexed with GRnH molecule.21 The effect of Cl−(II) to the reaction might be significant due to direct interaction between enzyme and the substrate molecule. However, for the current case, as our previous MD simulations indicated, there is no direct interaction between the substrate and the Cl−(II). However, about 13-fold activity difference between absence and presence of Cl−(II) has been observed for the hydrolysis of Hip-HisLeu.30 Indeed, it has long been recognized that the fact of the Cl− dependence of hydrolysis is substrate-specific.27,29,30 Therefore, it would be desirable to carry out MD simulations to investigate the substrate binding pattern and subsequent reactive mechanism without the chloride ion at the binding position. A total of 1.2 ns MD simulation was carried out to Model I to investigate the dynamics of the substrate in the active site without Cl−(II). The first 200 ps MD was used for heating and equilibration, and the rest of 1 ns trajectory was saved for final data analysis. The time-averaged geometric parameters for the substrate binding are listed in Table 3, and the corresponding

Figure 5. Snapshots of the stationary points along the reaction path for the hydrolysis of Hippuryl-His-Leu calculated with the SCCDFTB/MM method. Carbon atoms are colored in green, blue for nitrogen atoms, red for oxygen atoms, and purple for the zinc ion. Only the QM region is highlighted to allow a better view.

general base and general acid during the reaction. More importantly, the zinc bound water molecule is the key factor in the reaction, which acts as the nucleophile to attack the scissile carbonyl carbon atom assisted by the glutamic acid residue. The tetrahedral intermediate flanked by two transition states is proposed to connect the whole reaction path. The NA step is more likely the rate-limiting step. Interestingly, the CPA has F

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Y520, respectively. Different from the wild type, the oxygen atoms of carboxylate group switch their position when the hydrogen bonds formed between Q281 and K511. In addition, like the wild-type ACE, the scissile carbonyl oxygen atom (O2) is also hydrogen-bonded with H513 and Y523, but a weaker hydrogen bond between O2 and H513 can be found for the current chloride ion model, as shown in Table 2. Overall, little change was observed for the substrate binding pattern in the presence and absence of Cl−(II), which clearly indicates that Cl−(II) has minor contributions to the substrate binding, at least for the Hip-His-Leu molecule. Current MD simulation without Cl-(II) agrees well with our previous MD simulations of two small inhibitors binding, lisinopril and enalaprilat.22 Such a phenomenon is also consistent with the fact that the absence of Cl−(II) in the ACE homologue from D. melanogaster does not affect the binding of short peptide inhibitors.32,33 B. PMFs for the Hydrolysis Reaction. To further understand the influences on the reactive mechanism, we carried out freeenergy simulations for two designed chloride models, as we described above. Corresponding free-energy profiles are displayed in Figure 8. The removal of Cl− and scaling chloride

Table 3. Selected Geometric Parameters of the Michaelis Complex Obtained from the QM/MM MD Simulations Using the SCCDFTB/MM Method distance (Å) angle (deg.) Zn···Ow Zn···O2 Zn···Nε2 (H387) Zn···Nε2 (H383) Zn···Oε1 (E411) H1···Oε2 (E384) O2···Hε2 (H513) O2···Hη(Y523) C1···Ow O3···Hε2 (H353) O4···Hη (Y520) O4···Hζ1 (K511) O5···H21 (Q281) O4···H21 (Q281) H4···O (A354) Cl···Zn Cl−···Hε (R522) Cl−···Hη (Y224) Cl−···H22 (R522) Nε2(H383)···Zn···Ow Nε2(H387)···Zn···Ow Oε1(E411)···Zn···Ow

ES (without Cl−(II)) 2.02 4.31 2.02 2.00 2.05 1.25 2.41 1.91 2.99 1.90 1.65 2.08

± ± ± ± ± ± ± ± ± ± ± ±

0.06 0.31 0.06 0.06 0.06 0.04 0.35 0.20 0.25 0.18 0.09 0.31

2.17 ± 0.54 1.98 ± 0.16

108.3 ± 5.7 105.7 ± 5.2 128.5 ± 8.7

ref 22 2.04 4.51 2.01 2.00 2.05 1.34 2.01 1.98 3.19 1.89 1.69 1.71 2.02

± ± ± ± ± ± ± ± ± ± ± ± ±

0.06 0.36 0.06 0.05 0.07 0.14 0.36 0.25 0.26 0.17 0.11 0.15 0.47

2.02 ± 0.19 10.33 ± 0.26 1.94 ± 0.14 1.83 ± 0.11 2.00 ± 0.21 100.9 ± 5.4 105.1 ± 6.2 110.1 ± 7.6

geometric parameters for the native enzyme−substrate complex are also included for comparison. We also plotted the overlap between one of snapshot extracted from the MD trajectory of Model I and the wild type in Figure 7. Without the Cl−(II) at

Figure 7. Overlap representation of the snapshots from the MD simulations of the wild type (carbon in green and zinc ion in purple) and the model without the chloride ion at its second binding position (carbon in cyan and zinc ion in yellow).

Figure 8. Potentials of mean force for the model I (panel A) and model II (panel B) calculated using combined SCC-DFTB/ CHARMM method.

the binding position, the tetra-coordination of the zinc ion was still well-maintained. In particular, along with three ligands provided from protein, H383, H387, and E411, the fourth ligand to zinc ion is the water molecule with the dZn···Ow = 2.02 ± 0.06 Å, which can be in good overlap with native enzyme. This water molecule is also hydrogen-bonded with E384 with a strong hydrogen bond of 1.25 ± 0.04 Å for the distance of H1···Oε2, which agrees with 1.34 ± 0.14 Å for this distance in the wild type very well. The substrate molecule is anchored in the active site by a strong hydrogen-bond network; for example, the terminal carboxylate group is stabilized by Q281, K511, and

ion charge to zero do not change the overall shape of the freeenergy profiles of the hydrolysis. However, the energy barrier heights for the first step of NA are calculated to be 19.4 kcal/ mol for Model I and 20.4 kcal/mol for Model II, respectively, while the second step of the cleavage of C−N bond is relatively shallow for both models. Compared with the energy barrier we obtained for the native enzyme, the relative changes of the free energy for the NA step are calculated to be 1.6 (Model I) and 2.6 kcal/mol (Model II), respectively. It should be noted here G

dx.doi.org/10.1021/jp400974n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

that the enzyme used in this article is the testis ACE rather the human somatic ACE. Therefore, the quantitative agreement of the change of the energy barrier height might not be expected compared with the native enzyme. However, we have to emphasize that the removal of the chloride ion from the second binding position does significantly weaken the enzymatic activity according to our simulation, which is consistent with experimental observations. In addition, because of the close activation barrier heights for the first step of two models, we might conclude that the major contribution of the Cl−(II) ion to the reaction might come from the long-range electrostatic interaction with the substrate. Our conclusions are not inconsistent with previous experimental work on the Cl− dependence, in which the authors concluded that the activation of C-domain via the Cl− should not be attributed to an allosteric mechanism but to the anion−substrate interactions.29 It should also be noted that the intrinsic character in the treatment of long-range electrostatic interaction in the stochastic boundary condition might cause some unexpected results. To more accurately predict the extent of how the chloride ion affects the whole enzymatic reaction, it would be better to perform some more calculations using either periodic boundary condition (PBC) or generalized solvent boundary potential (GSBP)64 in the investigation of the contributions of both chloride ions because they have better description of long-range electrostatic interactions. Because the stochastic boundary condition is used in this work, it is thus straightforward to do some more examinations for the contribution of chloride ion using the GSBP method, which has been implemented to CHARMM within SCC-DFTB framework in 2005 by Schaefer and Cui.65 In this work, using the same GSBP setup protocol in ref 65, we then investigated two models, wild-type ACE (wtACE) and the system without Cl−(II) (mACE), to see how extent the chloride ion affects the overall reaction. Only minimum energy path calculation combined with CPR optimization was carried out. In the GSBP simulation, the system is set up with an inner region of an 18 Å radius, which is sufficient to include the electrostatic interaction contribution to the active site from the chloride ion at the second binding position because it is ∼10 Å away from the zinc ion. To obtain the electrostatic potential and the reaction field matrix, the dielectric constant for solvent region is assigned to be ε = 80, while ε = 4 for protein environment. The boundary between protein and solvent was set using the atomic Born radii by Nina et al.66 The Poisson−Boltzmann (PB) equation is solved using a 1113 cubic grid centered on the zinc ion, with a coarse grid spacing of 1.2 Å. The 20th order spherical harmonics were applied to express the inner-region charge density with a total of 400 basis functions. We first obtained the minimum energy paths for both models using the adiabatic mapping approach, as we described above. Then, we carefully optimized the minima with a stringent convergence criterion of the gradient