Free Energy Contribution Analysis Using Response Kernel

Article pubs.acs.org/JPCB

Free Energy Contribution Analysis Using Response Kernel Approximation: Insights into the Acylation Reaction of a Beta-Lactamase Toshio Asada,*,†,‡ Kanta Ando,† Pradipta Bandyopadhyay,§ and Shiro Koseki†,‡ †

Department of Chemistry, Graduate School of Science, Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai 599-8531, Japan ‡ The Research Institute for Molecular Electronic Devices (RIMED), Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai 599-8531, Japan § School of Computational and Integrative Sciences, Jawaharlal Nehru University, New Delhi, India 110067 S Supporting Information *

ABSTRACT: A widely applicable free energy contribution analysis (FECA) method based on the quantum mechanical/molecular mechanical (QM/MM) approximation using response kernel approaches has been proposed to investigate the influences of environmental residues and/or atoms in the QM region on the free energy profile. This method can evaluate atomic contributions to the free energy along the reaction path including polarization effects on the QM region within a dramatically reduced computational time. The rate-limiting step in the deactivation of the β-lactam antibiotic cefalotin (CLS) by β-lactamase was studied using this method. The experimentally observed activation barrier was successfully reproduced by free energy perturbation calculations along the optimized reaction path that involved activation by the carboxylate moiety in CLS. It was found that the free energy profile in the QM region was slightly higher than the isolated energy and that two residues, Lys67 and Lys315, as well as water molecules deeply influenced the QM atoms associated with the bond alternation reaction in the acylenzyme intermediate. These facts suggested that the surrounding residues are favorable for the reactant complex and prevent the intermediate from being too stabilized to proceed to the following deacylation reaction. We have demonstrated that the free energy contribution analysis should be a useful method to investigate enzyme catalysis and to facilitate intelligent molecular design.

1. INTRODUCTION

approaches, referred to as the FEG-NEB method, was successfully applied to the chemical reaction in triosephosphate9,10 or the isomerization reaction of glycine in water.6 It is widely known that high-quality potential energy functions are essential for the reliable calculation of a given property of a system using quantum mechanical calculations. Although the QM/MM method6,11−13 dramatically reduces the computational time required for large systems, it is still extremely time-consuming to estimate FEGs in the QM region if ab initio molecular orbital (MO) calculations are employed for QM calculations. Consequently, cheaper

Enzymes are widely distributed in the living body to accelerate bond-breaking or bond-forming reactions1−3 under mild conditions. To investigate enzymatic reaction mechanisms, the reaction pathway on the free energy surface (FES) should be obtained taking into consideration both solvation and entropy effects at a given temperature. Computer simulation techniques have recently been expanded to treat explicit solvent molecules with their thermal fluctuations to understand the fundamental nature of complex systems. The free energy gradient (FEG) technique proposed by Nagaoka et al.4 and the nudged elastic band (NEB) method,5−8 which is usually used to optimize the reaction path on the energy surface, enable the optimization of the reaction pathway on a multidimensional FES. A combination of these © XXXX American Chemical Society

Received: June 16, 2016 Revised: August 7, 2016

A

DOI: 10.1021/acs.jpcb.6b06104 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Scheme 1. Schematic Representation of the Proposed Mechanism of Enzymatic Deactivation of a β-Lactam Antibiotic by a β-Lactamase

on the hydrolysis mechanism is necessary to design better antibiotics that are not susceptible to β-lactamases. Which amino acid(s) actually activate(s) Ser64 is a matter of debate. Recently, one group24 investigated the role of Lys67, Tyr150, and the β-lactam substrate by performing pKa calculations and molecular dynamics (MD) simulations, which suggested that hydrogen in Ser64 of the β-lactamase is activated by the nearby carboxylate group in CLS rather than the amino acid residues. In general, the acylation reaction has been considered to be the rate-limiting step in the case of the CLS substrate.26 We have therefore applied the FECA method to the carboxylate-activated acylation reaction of CLS by β-lactamase.

approximations such as semiempirical or tight-binding density functional theory calculations14 have been used for the QM regions, even in QM/MM calculations. Morita et al.15−18 proposed mathematical expressions to represent the polarized charges of molecules based on charge response kernel approximations19−21 instead of conventional ab initio MO calculations. In these expressions, polarized charges are described by the change in the electrostatic potentials on atomic sites affected by environmental charge distributions. Unfortunately, there exist potential disadvantages for some molecules, especially linear or planar ones.7,8 To overcome these difficulties, we have recently proposed an efficient computational approach to obtain reliable FEGs using charge and atom dipole response kernel (CDRK) approximations for the polarizable QM region in the context of a QM/MM hybrid computational scheme,7,8 which is denoted as the QM(CDRK)/MM method in this manuscript. In the QM(CDRK)/MM method, the whole system is divided into QM and MM regions similar to the conventional QM/MM calculations; however, energies and forces in the QM region are evaluated using prepared response kernels that are constructed using ab initio MO calculations prior to the simulation instead of expensive MO calculations. Consequently, the computational time required for reaction path optimizations can be reduced, and the FEG-NEB method, which leads to detailed reaction processes in molecular assemblies from a microscopic point of view, can be rendered feasible for larger molecular systems. In this manuscript, the free energy contribution analysis (FECA) method is formulated on the basis of the CDRK approximation so as to address important surrounding molecules or QM atoms affecting the free energy profile. To demonstrate the efficiency of the FECA method, it has been applied to the enzymatic reaction of a class C β-lactamase.22−27 β-Lactamases are produced by bacteria and protect bacterial cell wall synthesis from antibiotics such as penicillin or CLS by hydrolyzing the β-lactam ring. Drug resistance resulting from class C β-lactamase inhibitors has become a major concern. A general mechanism for the hydrolysis of the β-lactam ring has been proposed, as illustrated in Scheme 1;28,29 however, a detailed understanding of the effect of the surrounding residues

2. THEORY AND COMPUTATIONAL METHODS 2.1. The Reaction Coordinates on the Free Energy Surface. The reaction coordinates in the QM region can be obtained by optimizing structures using the FEG-NEB method. Since the detailed techniques used to optimize reaction paths on the FES with the QM(CDRK)/MM method were described in our previous paper,7,8 only a brief outline is described. The FEG4 on QM atoms can be expressed as ∂A = ∂rQM

∂E ∂rQM

rMM

(1)

where A is the free energy of the system, E is the total energy of the system, and ⟨X⟩rMM means the ensemble average of X over MM configurations at a given temperature. In order to evaluate QM polarizations induced by MM charge distributions, MO calculations should usually be carried out to achieve selfconsistent polarized MOs. Unfortunately, this is not practicable for most FEG calculations unless the QM regions are quite small because reliable FEGs require a large number of computations of QM forces at different MM configurations. In contrast, by adapting CDRK approximations to the QM region, reliable electrostatic polarizations can be obtained in a drastically reduced computational time by evaluating only the external electrostatic potentials and fields on QM atoms.7 In the NEB approach, the reaction path is represented by a chain of conformations connecting the reactant to the product. B


Article

The Journal of Physical Chemistry B Each intermediate conformation is optimized using force components perpendicular to the vector that defines the reaction path. The reaction coordinates,6 ds, may be written as follows:

∑

ds =

This type of evaluation must be used because polarization interactions are not simple pairwise additive properties but rather many-body interactions.29−31 In eq 7, Q 0a and μ0a mean the charge and the atom dipole of the ath QM atom in a given reference state, respectively, while v(ra) and E⃗ (ra) are the external electrostatic potentials and fields, respectively, on the ath QM atom induced by all the MM residues. Here, response kernels in eq 7 can be written as the following expressions:

|d R j|2 (2)

j ∈ QM

d R j = ( mj dxj ,

mj dyj ,

mj dzj)

(3)

where mj and dRj denote the mass and mass-weighted coordinates of the jth QM atom, respectively. The highest-free energy transition state structure was chosen as the origin of the reaction coordinates. 2.2. Free Energy Contribution Analysis along the Reaction Coordinate. Once the reaction coordinate on the FES is optimized, the free energy profile can be obtained by integrating the FEGs along the reaction coordinate. The free energy difference between two adjacent structures in the NEB method can be approximated by the following equation if the QM structural displacements are very small ΔA =

∂A ·ΔrQM ∂rQM

χa , b = ζa , b = κa , b = γa , b =

(4)

∂EQM ∂rQM

= ΔEQM +

n

+

∑ i ∈ MM QM ∂Epol

∂rQM

∂Ei ∂rQM

·ΔrQM rMM n

+

∑ i ∈ MM

∂Ei ∂rQM

·ΔrQM (5)

rMM

∂v(rb)

∂μa ∂ 2⟨ψ |Heff |ψ ⟩ = ∂E ⃗(ra)∂E ⃗(rb) ∂E ⃗(rb) ∂Q a ∂rb ∂μa ∂rb

The second term in the right-hand side, is the polarization energy in the QM region caused by the MM molecules. The third term is the sum of the free energy contributions originating from the ith MM molecule. Here, Ei means the interaction energies, except for the polarization interaction between QM and the ith MM molecule, which is represented by

∂ 2⟨ψ |Heff |ψ ⟩ ∂E ⃗(ra)∂rb

∑

=

∂rQM − +

∑ i ∈ MM

{μa0

− MM ∂EiQM ,est

∂rQM

= − {Q a0 +

∑

∑ −

∑

κb , av(rb) −

∑ b ∈ QM

∂rQM

∑

QM − MM ∂(EjQM ) ,pol + E j ,est

∂rQM [(Q a0 +

a ∈ QM

−

∑

χab Δv(rb))Ei⃗ (ra)

b ∈ QM

∑

χab Δvi(rb){E ⃗(ra) − Ei⃗ (ra)}

b ∈ QM

− (μa0 +

∑

ζa , bΔE ⃗(rb))

b ∈ QM

−

χab Δv(rb)}E ⃗(ra)

+

∂E ⃗(ra) + ∑ ζa , bΔE ⃗(rb)} ∂ra b ∈ QM

b ∈ QM

(j ≠ i) ∈ MM

∑

=

ΔAjest + pol

∑

QM − MM ∂(EjQM ) ,pol + E j ,est

j≠i

(6)

b ∈ QM

γb , aE ⃗(rb)

ΔAjest + pol −

j

The interaction energy, Ei, is the sum of the bonding energy, the van der Waals interaction, EQM−MM , and the i,vdW electrostatic interactions, EQM−MM , between the QM and the ith i,est MM molecule. While EQM−MM /∂rQM and EQM−MM /∂rQM are i,bond i,vdW easily evaluated from MM potentials because EQM−MM and i,bond EQM−MM are represented by conventional MM interactions i,vdW similar to the conventional QM/MM scheme, ∂EQM pol /∂rQM and EQM−MM /∂rQM should be evaluated using CDRK approximations i,est as follows: +

(8)

j ∈ MM

EQM−MM , i,bond

QM ∂Epol

∂ 2⟨ψ |Heff |ψ ⟩ ∂v(ra)∂rb

=

=

ΔAiest + pol =

EQM pol ,

− MM − MM − MM Ei = EiQM + EiQM + EiQM ,bond ,vdW ,est

∂ 2⟨ψ |Heff |ψ ⟩ ∂v(ra)∂v(rb)

=

The right-hand side includes the polarized wave function ψ and the effective Hamiltonian Heff in the QM region with MM charges; thus, the free energy component between the QM and MM atoms can be obtained from eqs 5, 6, and 7. It is noteworthy that these terms can be evaluated within a short duration, as shown in our previous paper,7,8 using site−site CDRK approximations and MM charge distributions. Carrying out these free energy decompositions using conventional QM/MM calculations is difficult, since atomic contributions to the polarization energies in the QM region are difficult to separate from delocalized MOs. In this paper, the sum of the free energy contributions of both electrostatic and polarization interactions between the QM and ith MM molecules is denoted as ΔAest+pol , which can be derived i from eq 7 as

where A and rQM represent the free energy and QM coordinates, respectively. Since the evaluated free energy profiles will be somewhat varied between forward integrations or backward integrations, mean free energy profiles will be estimated with error bars indicating deviations from mean values for forward and backward integrations. If there are n MM molecules surrounding the QM region, according to eq 1, the free energy difference can be written as follows: ΔA =

∂Q a

∂Ei⃗ (ra) ∂ra

⎛ ∂E ⃗(r ) ∂E ⃗ (r ) ⎞ a − i a ⎟ ζa , bΔEi⃗ (rb)⎜ ∂ra ⎠ ⎝ ∂ra b ∈ QM

∑

∑ b ∈ QM

κb , avi(rb) −

∑ b ∈ QM

γb , aEi⃗ (rb)]

·Δra rMM

(9)

where Δvi(rb) and ΔE⃗ i(rb) represent the change in v(rb) and E⃗ (rb), respectively, on the bth QM atom originating from charges

(a ∈ QM) (7) C


Article


Figure 1. Structure of the complex formed between β-lactamase and cefalotin. (a) β-Lactamase and cefalotin are drawn as ribbon- and ball-and-stickstyle representations, respectively. (b) Atoms shown in ball-style were assigned to the QM region, while remaining atoms were assigned to the MM region.

Figure 2. Comparison of the fluctuation in the energy gradient on QM atoms between QM(CDRK)/MM (blue) and QM/MM (red) methods for different MM configurations, which were obtained using a MD simulation of the MM region using the CDRK model. Only the three largest fluctuations are presented: (a) Fx on C7, (b) Fy on C10, and (c) Fz on O4 in the reactant QM structure.

in the ith MM molecule. This contribution analysis is described as the FECA method, which can be used to find important molecules in the MM region and/or important atoms in the QM region for a given QM structure on the free energy profile. As polarization interactions are not pairwise additive properties, over the ith MM molecule is not exactly the sum of ΔAest+pol i the same as the sum of the whole electrostatic and polarization contributions to the free energy. For example, the polarized charge on the ath QM atom induced by the ith MM charge influences the electrostatic interactions between the ath QM and the jth (j ≠ i) MM atoms, in addition to those between the ath QM and the ith MM atoms. In the actual calculation, free energy components between the ath QM and the ith MM atoms were obtained by zeroing out the ith MM charges to satisfy eq 9. If a pairwise additive decomposition analysis is preferred, QM charges should not be affected by the modification of MM charges. In this manuscript, the pairwise additive component analysis is also performed with fixed QM charges, which are fully polarized charges under electrostatic fields for all MM atoms. 2.3. Model Structure Preparation. The X-ray crystal structure of AmpC β-lactamase complexed with the CLS substrate at a resolution of 1.53 Å, whose PDB ID code is 1kvl,32 was used as the initial structure. Since 1kvl is a S64G mutant, Gly64 was replaced with Ser64 for modeling the wild-type structure; hydrogen atoms were added to the heavy atoms in the crystal structure using the Leap module of the Amber9 package.33 The prepared structure was surrounded by TIP3P34 water molecules (11,484 molecules) within 8.0 Å of a bounding box around the

complex, and the net charge on the system was neutralized by adding chloride anions. The prepared structure was equilibrated by MD simulation at room temperature, in which AMBER99 force fields35 and GAFF parameters36 were used for amino acids and CLS, respectively. In the prepared structure, the active sites, which included Ser64 in the β-lactamase, the β-lactam ring moiety, and the carboxylate group in CLS as indicated by the ball-and-stick model in Figure 1b, were assigned to the QM region. The remainder of the system was treated as the classical MM region. 2.4. Molecular Dynamics Simulation. To obtain FEGs in the QM region, QM(CDRK)/MM MD simulations were performed at a temperature of 300 K and a pressure of 1 atm using the Amber-Gaussian interface program37 based on the Amber 933 and Gaussian 0938 packages. The response kernels of the QM region were prepared at the M0639/6-311++G(d) level of theory. For each optimization step in the QM region, 10 ps QM(CDRK)/MM MD simulations with a frozen QM structure were carried out for the thermal equilibration using a time step of 1.0 fs, and then 50 ps MM MD simulations followed to obtain FEGs on QM atoms using eq 1. Long-range electrostatic interactions were evaluated by means of the particle mesh Ewald method, where Lennard-Jones forces and real-space electrostatics were cut off at 8.0 Å.

3. RESULTS AND DISCUSSION Figure 2 shows the time development of the calculated energy gradients on the QM atoms obtained using ab initio MO calculations D


Article

The Journal of Physical Chemistry B for 1 ns. The values obtained using the QM(CDRK)/MM approach accurately reproduced those obtained using the more expensive QM/MM calculations. Since FEG calculations used the frozen QM structure, the fluctuations in QM forces come from the differences between the MM configurations. These facts suggest that the evaluated FEGs should be sufficiently reliable for use in long time frame simulations despite including a large number of environmental molecules. The computational times required to evaluate gradients for 100 000 configurations on 8-core CPUs using an Intel Xeon processor (X5570 2.93 GHz) were 2434 h and 2.6 h using the QM/MM and QM(CDRK)/MM methods, respectively. By applying both the QM(CDRK)/MM and FEG-NEB methods, chains of conformations representing the acylation reaction path have been successfully optimized. All optimized structures were provided in the Supporting Information. The root-mean-square deviations of the FEGs for the optimized structures were less than 0.2 kcal/(mol Å), which indicates that the QM(CDRK)/MM method is a powerful technique to optimize reaction paths within a reduced computational time. Since the reaction path on the FES could be optimized, FEGs on the QM atoms were integrated along the reaction path according to eq 4 to obtain the free energy profile of the reaction (Figure 3).

Figure 4. Interatomic distances in the QM region along the optimized reaction path. Atomic labels are given in Figure 1b.

Figure 5. Natural population of QM atoms along the reaction path.

In the subsequent step, a small energy barrier of only 0.2 kcal/mol could be found for the ring opening reaction of the β-lactam ring, where R(N1−C10) increased from 1.65 Å (3.4 amu1/2 bohr) to 2.43 Å (17.0 amu1/2 bohr). The structure indicated as TS2′ in Figure 3 was not characterized as the real TS, since there is only one structure between TS1 and TS2′ in our NEB calculations with a small free energy barrier. At this step, the charges on N1, C2, and S1 decreased, and the electron moved from the N1−C10 bond to the N1 atom, corresponding to the bond-breaking reaction. The electron localization on N1 caused a resonance in the charge density in the thiazine ring that decreased the charges on C2 and S1, in accordance with the ortho−para orientation rule. Following this step, HG moved from O3 in the carboxylate moiety in CLS to N1 through the transition state TS3 to complete the acylation reaction, leading to the stabilized product with a relative free energy of −14.6 ± 2.4 kcal/mol based on the reactant energy. As the proton transfer reaction proceeded, the natural charge at N1, as well as C2 and S1, increased as expected from the ortho−para orientation rule. Although TS3 had a relatively high activation barrier, the reverse reaction from the intermediate structure at 17.0 amu1/2 bohr to the reactant is apparently more unfavorable than the forward reaction. Thus, the acylation reaction was confirmed to proceed once Ser64 was activated by the nearby carboxylate moiety in CLS. To determine the important residues in the β-lactamase, FECAs were performed, as shown in Figure 6. These include free energy components for nine important residues and water, in addition to the isolated energies of the QM region. The contribution of water indicates the solvation effects coming from all surrounding water molecules. Although the shape of the free energy profile was mainly explained by the isolated QM energy,

Figure 3. Relative free energy profile of the acylation reaction in water along the optimized reaction path on the free energy surface. The optimized reaction path was represented by 16 structures in the NEB method. Error bars were evaluated from the difference between the forward and reverse thermal integration results.

The free energy barrier was predicted to be 26.0 ± 0.5 kcal/mol, which was slightly higher than the experimental barrier of 24.7 kcal/mol.40 Detailed structural changes along the reaction process were derived from the interatomic distances given in Figure 4. The initial process involving the highest transition state, TS1, was proton transfer from OG(Ser64) to O3 of the CLS carboxylate moiety (see Figure 1b) corresponding to the reaction coordinates from −18.0 to 3.4 amu1/2 bohr. Simultaneously, R(OG−C10) shortened from 2.98 to 1.46 Å to form the acyl bond between CLS and the β-lactamase. Since R(HG−O3) at 0.0 amu1/2 bohr was 1.05 Å, which was equivalent to a standard OH bond length in a strong hydrogen bond,41,42 the proton HG was confirmed to be tightly bound to the O3 atom. As shown in Figure 5, charges at both O3 and O4 were increased at TS1, since the positively charged HG moved from OG to O3. On the contrary, the negative charge increased on O5 in the β-lactam ring, as was suggested for intermediate 1 in the previously proposed mechanism illustrated in Scheme 1. E


Article


atoms that are strongly affected by the surrounded residues. Figure 8 shows QM atomic contributions affected by water and the three most important residues illustrated in Figure 6. In order to ensure that the sum of the atomic contributions could reproduce the total free energy profile, pairwise free energy contributions were evaluated. For hydration free energies, Figure 8a indicates that destabilization effects come from the HG atom in the QM region, while the largest stabilization effect was attributed to S1, which was placed at the outside of CLS and was exposed to the solvent. In contrast, HG was located at the interface region of the active site between CLS and β-lactamase. As a result of the large structural displacement of HG from OG to N1 in the QM region, the contribution of HG was influenced not only by water (Figure 8a) but also by the nearest charged residues such as Lys67 and Lys315. As the reaction proceeds, HG moves to escape from the positively charged Lys67 and Lys315 residues, resulting in the large stabilization of the free energy that could be explained by a decrease in the electrostatic repulsive interactions (see Figures 7 and 8b). It is noteworthy that relative free energies were evaluated by multiplying FEGs by atomic displacements, as given in eq 4, which implies that the amount of each QM atomic displacement during the reaction drastically affects the relative free energies as well as the strengths of the FEGs on the QM atoms. The unfavorable contribution from Lys67 was emphasized at C10 after the β-lactam ring was opened. This is because the atomic charge on C10 was large and increased after the ring opening reaction to reach a final value of ca. +0.87e, as seen in Figure 5. The interatomic distance between C10 and Lys67 was only 5 Å, which is sensitive to the interaction energy gradients. Similarly, positive contributions from Lys315 were obtained at the OG and C10 atoms, in which the OG charges, in particular, increased during the HG transfer reaction. Unfavorable contributions at C10 by both Lys67 and Lys315 suggested that, on average, C10 moved toward these residues during the simulation as the reaction proceeded. In our free energy component analysis, the negatively charged carbonyl group was stabilized by −3.3 kcal/mol at 35.2 amu1/2 bohr and −2.7 kcal/mol at 13.3 amu1/2 bohr by formation of a hydrogen bond with the backbone NH group in Ser64 and Ala318, respectively, which indicated that the carbonyl O5 atom was situated in the oxyanion hole, as suggested by Tripathi et al.22,23 This stabilization effect has been commonly observed in enzymatic acylation processes in a wide range of serine proteases.

Figure 6. Free energy contributions of nine important residues in βlactamase and water molecules. The relative energies of the QM region are also indicated.

especially for the free energy barrier height at TS1, this does not mean that the surrounding residues and water do not affect the free energy profile. For example, there were large contributions from water at the active site after the β-lactam ring was opened to form the acyl−enzyme intermediate. The hydration effects lowered the relative free energies by −12.5 kcal/mol for 19.8− 23.7 amu1/2 bohr. While positively charged Lys67 and Lys315 largely destabilized the free energies by ca. 11.6 and 8.5 kcal/mol, respectively, the negatively charged Glu272 stabilized the free energy by ca. −5.6 kcal/mol after the HG(Ser64) atom was transferred from OG(Ser64) to O3 in CLS. A stereo view of the location of important residues within 6.9 Å of the closest QM atoms in CLS is illustrated in Figure 7. Lys67, Lys315, and Glu272 are all located a similar distance from the CLS substrate. On the whole, the free energy barrier at TS1 was slightly higher than the isolated QM energy, and the free energy of the resultant acyl−enzyme intermediate was also higher than the QM energy, which suggested that the surrounding residues prevent the produced intermediate structure from being too stabilized to easily proceed to the following deacylation reactions. Unfortunately, it is difficult to modify residues in the β-lactamase structure, since the sequence alignments are coded in the bacterium genome as vital functions. On the other hand, if detailed contributions of QM atoms to the free energy profile could be obtained, we could design structures in antibiotics such as CLS or penicillin that are not easily hydrolyzed by β-lactamases, thereby producing improved antibiotics. We have therefore decomposed the free energy components into QM atomic contributions using the FECA method and eqs 5, 6, and 9 to find important QM

Figure 7. Stereo view of the relative position of important residues. F


Article


Figure 8. Relative free energy contributions of QM atoms along the reaction path affected by water and important residues: (a) water, (b) Lys67, (c) Lys315, and (d) Glu272.

and S1 atoms affected by Lys67 and water molecules, respectively; in contrast, unfavorable effects come from OG and HG affected by Lys67 and water, respectively. Molecular design focused on these QM atoms will be able to effectively influence the free energy profile while taking into account the QM atomic displacements and electric fields arising from environmental residues. We have successfully demonstrated that the FECA approach, in combination with the QM(CDRK)/MM model, should be a useful and widely applicable method to derive important atoms in the free energy profile and facilitate intelligent molecular design.

In spite of the difficulties encountered in decomposing the free energy profile into pair additive components, important QM atoms could be addressed through the FECA approach proposed in this manuscript. Actually, the C10, S1, OG, and HG atoms were controlled by closely located Lys67, Lys315, and water molecules to determine the relative free energy profile in the acylation reaction of the acyl−enzyme intermediate. Accordingly, molecular design focused on these QM atoms will be able to effectively influence the free energy profile along the reaction based on QM atomic displacements and electric fields from environmental residues. In this manuscript, we have successfully demonstrated that the FECA approach in combination with the QM(CDRK)/MM model should be a useful method to derive important atoms that affect the free energy profile and facilitate intelligent molecular design taking into consideration polarization and atomic displacements. The FECA approach proposed here is quite general and applicable for a wide range of large molecular systems.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b06104. Figure showing optimized structures of the QM region along the reaction coordinate on the free energy surface by the FEG-NEB method (PDF)

4. CONCLUDING REMARKS The FECA method in the context of the QM/MM approach was proposed on the basis of CDRK approximations to derive the free energy components of surrounding residues and QM atoms affecting the free energy profile, while still considering QM polarization effects and atomic displacements. The reaction path on the FES can be optimized by applying the QM(CDRK)/MM approach in conjunction with the FEG and NEB techniques. Since the CDRK approach enables a dramatic reduction in the computational time required to evaluate reliable FEGs on QM atoms, reaction path optimizations of the acylation reaction in a β-lactamase−CLS complex have been performed to obtain information about the drug resistance mechanism of class C β-lactamase inhibitors. The energy barrier of the acylation reaction was calculated to be 26.0 ± 0.5 kcal/mol, which is slightly higher than the experimental value of 24.7 kcal/mol. The FECA results confirmed that the relative free energies along the optimized reaction path were more unfavorable than the isolated QM energies except for the reactant. This fact suggests that the surrounding residues prevent the reaction intermediate from being too stabilized to proceed to the following deacylation reactions. Large stabilization components come from the HG

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by a Grant-in-Aid for Science Research from the Ministry of Education, Culture, Sport, Science and Technology (MEXT) in Japan (No. 16K05664), and also by the FLAGSHIP2020, MEXT within the priority study 5 (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use) by the MEXT program “The Strategic Program for Innovation Research (SPIRE)” using computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp160225). G


Article


■

the Charge Sensitivity of Aromatic and Nonaromatic Species. J. Am. Chem. Soc. 1997, 119, 4021−4032. (16) Ishiyama, T.; Morita, A. Analysis of Anisotropic Local Field in Sum Frequency Generation Spectroscopy with the Charge Response Kernel Water Model. J. Chem. Phys. 2009, 131, 244714. (17) Ishiyama, T.; Morita, A. Molecular Dynamics Study of Gas− Liquid Aqueous Sodium Halide Interfaces. I. Flexible and Polarizable Molecular Modeling and Interfacial Properties. J. Phys. Chem. C 2007, 111, 721−737. (18) Ishiyama, T.; Imamura, T.; Morita, A. Theoretical Studies of Structures and Vibrational Sum Frequency Generation Spectra at Aqueous Interfaces. Chem. Rev. 2014, 114, 8447−8470. (19) Isegawa, M.; Kato, S. Polarizable Force Field for Protein with Charge Response Kernel. J. Chem. Theory Comput. 2009, 5, 2809−2821. (20) Nakamura, H.; Ohto, T.; Nagata, Y. Polarizable Site Charge Model at Liquid/Solid Interfaces for Describing Surface Polarity: Application to Structure and Molecular Dynamics of Water/Rutile TiO2(110) Interface. J. Chem. Theory Comput. 2013, 9, 1193−1201. (21) Nagata, Y. Polarizable Atomistic Calculation of Site Energy Disorder in Amorphous Alq3. ChemPhysChem 2010, 11, 474−479. (22) Tripathi, R.; Nair, N. N. Mechanism of Acyl−Enzyme Complex Formation from the Henry−Michaelis Complex of Class C βLactamases with β-Lactam Antibiotics. J. Am. Chem. Soc. 2013, 135, 14679−14690. (23) Tripathi, R.; Nair, N. N. Deacylation Mechanism and Kinetics of Acyl−Enzyme Complex of Class C β-Lactamase and Cephalothin. J. Phys. Chem. B 2016, 120, 2681−2690. (24) Sharma, S.; Bandyopadhyay, P. Investigation of the Acylation Mechanism of Class C beta-lactamase: pKa Calculation, Molecular Dynamics Simulation and Quantum Mechanical Calculation. J. Mol. Model. 2012, 18, 481−492. (25) Díaz, N.; Suárez, D.; Sordo, T. L. Molecular Dynamics Simulations of Class C beta-Lactamase from Citrobacter freundii: Insights into the Base Catalyst for Acylation. Biochemistry 2006, 45, 439−451. (26) Lietz, E. J.; Truher, H.; Kahn, D.; Hokenson, M. J.; Fink, A. L. Lysine-73 Is Involved in the Acylation and Deacylation of β-Lactamase. Biochemistry 2000, 39, 4971−4981. (27) Chen, Y.; McReynolds, A.; Shoichet, B. K. Re-examining the role of Lys67 in Class C β-Lactamase Catalysis. Protein Sci. 2009, 18, 662− 669. (28) Goldberg, S. D.; Iannuccilli, W.; Nguyen, T.; Ju, J.; Cornish, V. W. Identification of Residues Critical for Catalysis in a Class C β-lactamase by Combinatorial Scanning Mutagenesis. Protein Sci. 2003, 12, 1633− 1645. (29) Asada, T.; Nishimoto, K.; Kitaura, K. Theoretical Study on Binding Enthalpies and Populations of Isomers of Cl−(H2O)n Clusters at Room Temperature. J. Phys. Chem. 1993, 97, 7724−7729. (30) Kitaura, K.; Sawai, T.; Asada, T.; Nakano, T.; Uebayashi, M. Pair Interaction Molecular Orbital Method: An Approximate Computational Method for Molecular Interactions. Chem. Phys. Lett. 1999, 312, 319− 324. (31) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayashi, M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecule. Chem. Phys. Lett. 1999, 313, 701−706. (32) Beadle, B. M.; Trehan, I.; Focia, P. J.; Shoichet, B. K. Structural Milestones in the Reaction Pathway of an amide hydrolase: Substrate, acyl, and Product Complexes of Cephalothin with AmpC betalactamase. Structure 2002, 10, 413−424. (33) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R. K.; Merz, M.; Pearlman, D. A.; Crowley, M.; et al. AMBER 9; University of California: San Francisco, CA, 2006. (34) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (35) Wang, J.; Cieplak, P.; Kollman, P. A. How Well does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating

ABBREVIATIONS FES, free energy surface; FEG, free energy gradient; NEB, nudged-elastic band; QM/MM, quantum mechanical/molecular mechanical; CDRK, charge and atom dipole response kernels; MO, molecular orbital; FECA, free energy contribution analysis; CLS, cefalotin; MD, molecular dynamics

■

REFERENCES

(1) Abe, T.; Hashimoto, Y.; Zhuang, Y.; Ge, Y.; Kumano, T.; Kobayashi, M. Peptide Bond Synthesis by a Mechanism Involving an Enzymatic Reaction and a Subsequent Chemical Reaction. J. Biol. Chem. 2016, 291, 1735−1750. (2) Isaksen, G. V.; Hopmann, K. H.; Åqvist, J.; Brandsdal, B. O. Computer Simulations Reveal Substrate Specificity of Glycosidic Bond Cleavage in Native and Mutant Human Purine Nucleoside Phosphorylase. Biochemistry 2016, 55, 2153−2162. (3) Raich, L.; Borodkin, V.; Fang, W.; Castro-López, J.; van Aalten, D. M. F.; Hurtado-Guerrero, R.; Rovira, C. A Trapped Covalent Intermediate of a Glycoside Hydrolase on the Pathway to Transglycosylation. Insights from Experiments and Quantum Mechanics/ Molecular Mechanics Simulations. J. Am. Chem. Soc. 2016, 138, 3325− 3332. (4) Nagaoka, M.; Okuyama-Yoshida, M.; Yamabe, T. Origin of the Transition State on the Free Energy Surface: Intramolecular Proton Transfer Reaction of Glycine in Aqueous Solution. J. Phys. Chem. A 1998, 102, 8202−8208. (5) Mills, G.; Jónsson, H. Quantum and Thermal Eff’ects in H2 Dissociative Adsorption: Evaluation of Free Energy Barriers in Multidimensional Quantum Systems. Phys. Rev. Lett. 1994, 72, 1124− 1127. (6) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Asada, T.; Nagaoka, M. Reaction Path Optimization and Vibrational Frequency Analysis via ab initio QM/MM Free Energy Gradient (FEG) method: Application to Isomerization Process of Glycine in Aqueous Solution. Theor. Chem. Acc. 2011, 130, 215−226. (7) Asada, T.; Ando, K.; Sakurai, K.; Koseki, S.; Nagaoka, M. Efficient Approach to include Molecular Polarizations using Charge and Atom Dipole Response Kernels to Calculate Free Energy Gradients in QM/ MM Scheme. Phys. Chem. Chem. Phys. 2015, 17, 26955−26968. (8) Asada, T.; Ando, K.; Koseki, S. Efficient Approach to obtain Free Energy Gradient using QM/MM MD Simulation. AIP Conf. Proc. 2015, 1702, 090005. (9) Lu, Z.; Yang, W. Reaction Path Potential for Complex Systems derived from Combined ab initio Quantum Mechanical and Molecular Mechanical Calculations. J. Chem. Phys. 2004, 121, 89−100. (10) Hu, H.; Lu, Z.; Parks, J. M.; Burger, S. K.; Yang, W. Quantum Mechanics/Molecular Mechanics Minimum Free-energy Path for Accurate Reaction Energetics in Solution and Enzymes: Sequential Sampling and Optimization on the Potential of Mean Force Surface. J. Chem. Phys. 2008, 128, 034105. (11) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO + MM Method for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels−Alder Reactions and Pt(P(t-Bu)3)2 + H2 Oxidative Addition. J. Phys. Chem. 1996, 100, 19357−19363. (12) Asada, T.; Ohta, K.; Matsushita, T.; Koseki, S. QM/MM Investigation on the Degradation Mechanism of the Electron Transporting Layer. Theor. Chem. Acc. 2011, 130, 439−448. (13) Xu, J.; Zhang, J. Z. H.; Xiang, Y. Ab Initio QM/MM Free Energy Simulations of Peptide Bond Formation in the Ribosome Support an Eight-Membered Ring Reaction Mechanism. J. Am. Chem. Soc. 2012, 134, 16424−16429. (14) Seifert, G. Tight-Binding Density Functional Theory: An Approximate Kohn−Sham DFT Scheme. J. Phys. Chem. A 2007, 111, 5609−5613. (15) Morita, A.; Kato, S. Ab Initio Molecular Orbital Theory on Intramolecular Charge Polarization: Effect of Hydrogen Abstraction on H


Article

The Journal of Physical Chemistry B Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (36) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General AMBER Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (37) Okamoto, T.; Yamada, K.; Koyano, Y.; Asada, T.; Koga, N.; Nagaoka, M. A Minimal Implementation of the AMBER−GAUSSIAN Interface for Ab Initio QM/MM−MD Simulation. J. Comput. Chem. 2011, 32, 932−942. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision E.01; Gaussian, Inc.: Wallingford, CT, 2009. (39) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−41. (40) Galleni, M.; Amicosante, G.; Frère, J. M. A Survey of the Kinetic Parameters of Class C Beta-lactamases. Cephalosporins and Other Betalactam Compounds. Biochem. J. 1988, 255, 123−129. (41) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond. Structural Chemistry and Biology; Oxford University Press: New York, 1999. (42) Asada, T.; Koseki, S.; Nagaoka, M. Ab initio Electron Correlation Studies of the Intracluster Reaction NO+(H2O)n → H3O+(H2O)n−2(HONO) at n = 4 and 5. Phys. Chem. Chem. Phys. 2011, 13, 1590−1596.

I


Free Energy Contribution Analysis Using Response Kernel

Recommend Documents