Theoretical Study on the Mechanism of the Acylate Reaction of β-Lactamase.

Using density functional theory and a cluster approach, we study the reaction potential surface and compute Gibbs free energies for the acylate reaction of β-lactamase with penicillin G, where the solvent effect is important and taken into consideration. Two reaction paths are investigated: one is a multi-step process with a rate-limit energy barrier of 19.1 kcal/mol, which is relatively small, and the reaction can easily occur; the other is a one-step process with a barrier of 45.0 kcal/mol, which is large and thus makes the reaction hard to occur. The reason why the two paths have different barriers is explained.

Introduction

β-Lactam antibiotics (penicillin, cephalosporin, natural and synthetic monocyclic β-lactam, carbapenems, oxapenams, carbacephems, and oxacephems)[1] are widely used antimicrobials because of their pharmacological advantages such as potent activity, low toxicity, ease of delivery, and low production costs.[2] The bacterial production of β-lactamases is the most common mechanism of resistance to β-lactam antibiotics.[1] As a group of proteins that hydrolyze antibiotics, β-lactamases play a key role in antibiotic resistance.[3] On the basis of primary amino sequence homology, β-lactamases can be divided into four classes (A, B, C, and D).[4,5] Class A β-lactamases are the major source of bacterial resistance to the β-lactam antibiotics, and they have been extensively studied mechanistically.[6] The most commonly encountered and best studied class A β-lactamase is TEM-1,[7] whose structure and potential catalytic mechanisms have been investigated extensively.[3,6,8−20] Using penicillin as an example, the overall hydrolysis process of TEM-1 can be divided into acylation of a serine at the active site of the enzyme by the β-lactam carbonyl group, generation of an acyl-enzyme intermediate, and deacylation of the enzyme by a solvent.[3,17]

There is a consensus on the deacylation step of the mechanism. However, some ongoing controversies on the details of the acylation mechanism have remained.[3,17] Thus, more research studies are focused on the understanding of the acylation mechanism. Early in 1987, Herzberg and Moult[21] determined the crystal structure of β-lactamase from the Gram-positive bacterium Staphylococcus aureus PC1 at a 2.5 Å resolution and revealed that the active site is located at the interface between two closely associated domains with the key catalytic residue Ser70 at the amino terminus of a buried helix. Three years later, Gibson et al.[22] constructed two single mutants and the corresponding double mutant of β-lactamase I from Bacillus cereus 569/H and investigated their kinetics. The result showed that the rate constants of both acylation and deacylation for the hydrolysis of benzylpenicillin are decreased about 2000-fold in this mutant. In 1992, Strynadka et al.[23] studied the X-ray crystal structure of the molecular complex of penicillin G with a diacylation-defective mutant of the RTEM-1 β-lactamase from Escherichia coli. The crystallographic analysis showed that Penicillin G is covalently bound to Ser70 Oγ as an acyl-enzyme intermediate. Strynadka et al. deduced that the catalytic mechanism uses Ser70 Oγ as the attacking nucleophile during acylation. A year later, the result of a kinetic study carried out by Knox et al.[24] suggested that Glu166 plays an important role in diacylation but is minimal or not involved in acylation.

In 1997, the molecular modeling (AMBER) studies done by Guillaume et al.[25] confirmed the crucial roles of Glu166 and of the “catalytic” water molecule in both the acylation and deacylation processes. Using a mean-field approach (AMBER) in 2000, Atanasov et al.[26] calculated the pH dependence of the pKa values of all ionizable groups and of the electrostatic potential at grid points corresponding to catalytically important atoms in the active site of TEM-1 β-lactamase, and they concluded that the catalytic action of TEM-1 β-lactamase results from a trigger effect of N-protonation through electrophilic general-acid catalysis to facilitate general-base-catalyzed nucleophilic attack by Ser70 to form the acylenzyme. Wang et al.[27] studied the clinically isolated mutants of the β-lactamase TEM-1 in 2002 and determined the X-ray crystallographic structures of three mutant enzymes. The results indicated that activity gain and stability loss are related to an enlarged active site cavity in the mutant enzymes. By molecular dynamics simulations (AMBER) of the TEM-1 β-lactamase in aqueous solution, Díaz et al.[14] proposed that either the Glu166 carboxylate-Water1 or substrate carboxylate-Ser130 moieties could abstract a proton from the nucleophilic Ser70. In 2005, the ab initio quantum mechanics/molecular mechanical (QM/MM) study of class A β-lactamase acylation performed by Meroueh et al.,[28] where MP2/6-31+G* was applied to calculate single-point energy based on HF/3-21G geometry optimization for the QM part and the AMBER force field was used for the MM region, showed that the pathway involving Glu166 as the general base promoting Ser70 through a conserved water molecule exists in competition with the Lys73 process. A recent study of the reaction between the antibiotic cefotaxime and the CTX-14 class A serine hydrolase by Lizana and Delgado with QM/MM (B3LYP-D3/6-31+G(d,p)/CHARMM36) calculations predicted that the reaction should occur via a concerted mechanism where the acylation of carbonyl carbon of the lactam, protonation of the N lactam atom, and the opening of the β-lactam ring take place simultaneously.[29]

Hence, different proposals have been postulated on the mechanism of the acylation of benzylpenicillanic acid by TEM-1 β-lactamase. Just as mentioned by Lizana and Delgado, the mechanism of acylation is still uncertain and controversial, and further information is of interest for the design of novel inhibitors and antibiotics of β-lactamase enzymes.[29] In this work, we study the potential energy surface (PES) of the acylate reaction of β-lactamase using density functional theory (DFT) calculations and seek to elucidate its mechanism.

Computational Details

Generally, when scientists theoretically study the enzyme reaction with tens of thousands of atoms, they apply the QM/MM method, where the active part is chosen to be calculated using quantum chemistry computation and the other part is treated with molecular mechanics. With the development of computers, the active part becomes larger, and the model with 250–300 atoms can be computed. The cluster approach is to choose the active part with many atoms and neglect the other part, which is utilized to investigate many enzyme-catalyzed reactions successfully.[30−43] Here, we apply this method and DFT computations to discuss the reaction mechanism for the acylate reaction of β-lactamase with the substrate penicillin G.

The starting cluster structure (R) is chosen based on a crystal structure (the RTEM-1 ß-lactamase with PDB ID code 1FQG(23)) corresponding to the acylate reaction of β-lactamase. The cluster includes the substrate penicillin G and the residues Ser70, Lys73, Lys234, Ser130, and Asn166. Some atoms are fixed at their crystallographic position to avoid the running away of the substrate and residues without a chemical bond when we optimize the cluster structure (see Figure ).

Figure 1

Computational model of the reaction center of the acylate reaction of β-lactamase and the optimized reactant (R) structure. The atoms fixed are marked with asterisks.

DFT calculations are performed with Gaussian 16 software.[44] Stationary points on the PESs are optimized with the M08HX functional and the 6-31G(d,p) basis set. The M08HX functional can investigate the intermolecular and intramolecular nonbonding interactions well.[45] In our previous investigations, the M08HX functional is utilized to study polypeptides and organometallic compounds, and good results are obtained.[46−48]

The corresponding zero-point vibrational energies (ZPVEs, in kcal/mol) and entropy (S, in cal mol–1 K–1) are obtained at the same level of theory. In the following discussion, the transition state, which has only one imaginary frequency, and the reaction intermediate are denoted as TS and IM, respectively. Relative Gibbs free energies (G, in kcal/mol) obtained at 298.15 K and 1 atm[6] are applied because the entropy effect may significantly change the PES. Generally, there are two methods of calculating the free energy in solution. The first method is to optimize the molecular structure in vacuum and then compute the free energy under the influence of a solvent based on the optimization, where ΔG(solvent) = ΔG(gas) + ΔG(solvation). This method is fast and can save computing resources. Its disadvantage is that the influence of solvation on the molecular structure is not considered. The second method is to optimize the molecular structure in the solvent using the polarizable continuum model (PCM).[49] After the optimization, the free energy is also computed with PCM. Here, we use the second method.

We compute the cluster in vacuum. Furthermore, to account for the dielectric effect of the environment, the solvation effect is calculated by PCM using the integral equation formalism variant (IEFPCM)[49] with three different dielectric constants ε (4, 24, and 78). Generally, the dielectric constant of protein is about 2–4.[50] We use a value of 4 for the protein environment. The dielectric constant of ethanol is about 24, and that of water is about 78.

Results and Discussion

Path A

The proposed mechanism of this reaction is presented in Figure . The potential energy diagram is depicted in Figure . The relative computational electronic energies and Gibbs free energies are listed in Table .

Figure 2

(a) Schematic representation of the mechanism of the acylate reaction of β-lactamase and (b) geometric parameters of TS3 using PCM with the dielectric constant of 4. Bond lengths are in angstroms.

Figure 3

PESs for the acylate reaction of β-lactamase in vacuum and using PCM with different dielectric constants of 4 (brackets), 24 (square brackets), and 78 (curly brackets). The relative Gibbs free energies are in kcal/mol.

Table 1

Entropy S (in cal mol–1 K–1), Relative Electronic Energies (RE, in kcal/mol) and Gibbs Free Energies (ΔG, in kcal/mol) at 298.15 K and 1 atm of the Reactant, Intermediates, Transition States, and the Product Calculated by M08HX/6-31G(d,p) in Vacuum and Using PCM with Three Different Dielectric Constants ε (4, 24, and 78)

ε value	species	S	RE	ΔG
vacuum	R	325.0	0.0	0.0
	IM1	322.8	12.2	11.6
	IM2	327.8	–17.4	–17.8
	TS1	322.8	25.5	23.5
	TS2	320.6	11.7	8.8
	TS3	319.1	53.5	51.6
4	R	326.4	0.0	0.0
	IM1	324.7	6.4	6.5
	IM2	329.1	–18.7	–18.8
	TS1	325.1	21.3	19.1
	TS2	320.1	7.8	5.7
	TS3	324.9	46.7	45.0
24	R	330.1	0.0	0.0
	IM1	328.2	3.8	4.2
	IM2	326.8	–19.1	–17.8
	TS1	321.2	19.7	19.2
	TS2	322.7	6.2	4.7
	TS3	326.2	45.9	45.3
water (78)	R	330.2	0.0	0.0
	IM1	324.2	3.3	4.5
	IM2	327.7	–19.1	–18.1
	TS1	329.7	19.4	17.1
	TS2	323.0	5.9	4.4
	TS3	326.3	45.7	45.3

Our results show that the title reaction is divided into two steps in this path. First, via a H-transfer process, R rearranges to IM1 over a transition state TS1 with an energy barrier of 19.1 kcal/mol when the dielectric constant is taken to be 4 (if not specified, the following dielectric constant is 4 because of its similarity to the protein environment). When we compare the free energies in the solvents with those in vacuum, we find that the solvent effect influences the free energies and the energy barriers. However, the differences of the free energies caused by different dielectric constants are small (see Figure ). For example, the reaction barriers are 23.5, 19.2, and 17.1 kcal/mol when the dielectric constant is taken as 1, 24, and 78, respectively. Comparing the calculated barriers in protein and water with that in vacuum, the reaction barriers decrease much, which makes the corresponding reaction to easily occur. From Figure , we know that this barrier is the largest in the title reaction. Thus, this transition state determines the reaction rate. In ref (51), Pitarch et al. applied the QM/MM (AM1/CHARMM) method to study the acylation reaction of penicillanate with TEM1 ß-lactamase. They found that this step forming an initial tetrahedral adduct is a rate-limiting step with a barrier energy of 18.3 kcal/mol. Our results are in good agreement with theirs. Although we only treat the enzyme reaction with a cluster, the reasonable results can be achieved. Furthermore, in ref (17), Xie et al. used the nudged elastic band method (NEB) and the QM/MM (AM1/AMBER) method to investigate this reaction, and the corresponding reaction barrier for this step is 26.0 kcal/mol, which is larger than the results calculated by us and Pitarch et al.[51] Just as pointed out by Xie et al., the crystallographic waters were taken into consideration in ref (51); in their study, only waters to solvate the active site were included and outer waters were considered with the charge scaling technique, which makes the electrostatic interaction between the QM part and the MM part weaker and the energy barrier larger.

In this process, the O2 and H3 atoms of Ser70 in the structure of R are added to the C1 atom of penicillin G and the N4 atom of Lys73 in the structure of IM1 to form the C1–O2 and H3–N4 bonds, respectively, along with the breaking of the O2–H3 and C1–N5 bonds in R. The Cartesian coordinates of the species on the PES are listed in the Supporting Information. In the geometry of TS1, the bond-forming C1···O2 and H3···N4 distances are shortened to 2.525 and 1.123 Å, respectively, which are 2.387 and 1.497 Å in ref (17), while the breaking O2–H3 and C1–N5 bonds are lengthened to 1.659 and 1.383 Å, respectively. The length of the breaking O2–H3 bond is 1.915 Å in ref (17). We find that the structure of the transition state TS1 obtained by our optimization is somewhat different from that in ref (17).

Secondly, IM1 isomerizes to IM2 through concerted H-transfer via the transition state TS2 with a barrier of −0.8 kcal/mol. The Gibbs free energy of the barrier of the transition state is abnormal, and its value is smaller than that of the reactant IM1. When the molecular structure is optimized, Gaussian software applies the forces on atoms, which are the derivatives of the electronic energy with respect to molecular coordinates. When the dielectric constant is 4, the relative electronic energies of IM1 and TS2 are 6.4 and 7.8 kcal/mol (see Table ), respectively. Thus, there is a transition state TS2 connecting the intermediate IM1 and the product IM2. However, what determines the reaction is not the electronic energy but the Gibbs free energy. The free energy of IM1 is 6.5 kcal/mol, which is larger than that of TS2 with a value of 5.7 kcal/mol (see Table ). Hence, IM1 is unstable, and the reaction will go directly to the product without a barrier. This indicates that it is very easy for IM1 to generate the product.

This phenomenon has appeared in many literature studies when scientists investigated molecular transition states.[47,52−55] The reason may be that the reaction energy barrier of the transition state is very small.

During this process, the H6 atom in IM1 transfers from N4 to O7, while the H8 atom shifts from O7 to N5. In the structure of TS2, the bond-forming H6···O7 and N5···H8 distances are reduced to 1.611 and 1.286 Å, respectively, while the breaking N4–H6 and O7–H8 bonds are elongated to 1.090 and 1.208 Å, respectively. The whole process is also exothermic by 18.8 kcal/mol. Our mechanism of this process is somewhat different from those in refs (17) and (51), where there is the two-step process: in the first step, the H8 atom of Ser130 transfers from O7 to N5 of penicillin G with a relatively low energy barrier of 5.9[17] and 2.8 kcal/mol,[51] respectively; in the second step, the H6 atom of Lys 73 in IM1 shifts from N4 to O7 of Ser 130 with a relatively low energy barrier of 1.7[17] and 3.6 kcal/mol,[51] respectively. All the computational results indicate that the energy barriers of this process are very small, which makes the reaction easy to occur.

Path B

In ref (29), Lizana and Delgado performed the computational study on the reaction between cefotaxime and the CTX-14 class A serine hydrolase. B3LYP-D3/6-31+G(d,p) was used to investigate the QM subsystem, where Grimme corrections were included to consider the dispersive (vdW) interactions, and the CHARMM36 force field was applied to study the MM region. They suggested a concerted mechanism where only one step is involved. In the reaction, the proton of the hydroxyl group of Ser70 is transferred to the N atom of cefotaxime, and the oxygen atom of the hydroxyl group attacks the C atom of the carbonyl group of the β-lactam ring. When the above two distances are about 2.0 and 1.2 Å, respectively, the reaction reaches a saddle point with an activation barrier of 20.0 kcal/mol.[29] They concluded that the reaction path is feasible and may compete with other mechanisms.[29]

In order to know whether our title reaction can occur through the concerted mechanism, we also perform the computational study. In this path, R transfers to IM2 directly through concerted forming and breaking of several bonds via the transition state TS3. The geometry of TS3 is presented in Figure b. In the structure of TS3, the bond-forming C1···O2 and H3···N5 distances are reduced to 1.687 and 1.361 Å, respectively, while the breaking C1–N5 and O2–H3 bonds are elongated to 1.558 and 1.229 Å, respectively. The structure of TS3 is quite different from that of the transition state in ref (29). The energy barrier of this process is 45.0 kcal/mol, which is 25.0 kcal/mol larger than the barrier in ref (29) and 25.9 kcal/mol higher than that of Path A (imposed by TS1, see Figure ). The substrate and enzyme in ref (29) are different from ours, and different structures of the transition states result in the different barriers. Thus, this path is difficult to occur.

We also perform computations on the electronic energies with M08HX/6-311+G(d,p), M08HX/6-311++G(d,p), and B3LYP-D3/6-311+G(d,p)[56] on the structures optimized by M08HX/6-31G(d,p) using PCM with the dielectric constant of 4 (see Table S2). When we compare the energies with those in Table , we find that the results calculated by M08HX with large basis sets agree with those by M08HX/6-31G(d,p). The influence of the basis set on the electronic energies is small. When the dispersion correction is included, the energies calculated by the B3LYP-D3 functional are in good agreement with those by the M08HX functional. Furthermore, we do the energies with AM1 (see Table S3). The semi-empirical molecular orbital theory approach is not as good as DFT. The electronic energies calculated by AM1 are quite different from those by M08HX. For example, the electronic energy of TS1 is 21.3 kcal/mol calculated by M08HX/6-31G(d,p) and 38.1 kcal/mol by AM1.

In order to know the forces on the frozen atoms, we list all of them and their absolute values when the M08HX/6-31G(d,p) functional is applied using PCM with the dielectric constant of 4 (see Table S4). The convergence criteria of optimizing the molecular structure are that the maximum force is 0.002500 hartrees/bohr and the root-mean-square (rms) force is 0.001667 hartrees/bohr. When the forces on the frozen atoms are smaller than the maximum force, the influences of these forces can be neglected. However, the forces on some frozen atoms are larger than the maximum force. For the species having different orientations, their X, Y, and Z components of the forces are somewhat different (see Table S4). When we compare the absolute values of the forces, we find that not all the forces on the same atoms of the reactant, intermediates, transition states, and the product are close. Thus, the frozen atoms may have a certain influence on the calculational results. However, the computational results show that the cluster approach may not result in large errors.

Why does Path A easily occur and Path B react with difficulty? The reason may be that the distances between the different heavy atoms, which connect the hydrogen atoms, are different. In Path A, the distances between the O2 and N4 atoms, the N4 and O7 atoms, and the O7 and N5 atoms are 2.893, 2.805, and 3.047 Å, respectively. In Path B, the distance between the O2 and N5 atoms is 3.093 Å, which is larger than the distances in Path A. The longer the distance between heavy atoms is, the larger the energy barrier may be required for H-atom transfer.

In our computations, there are 134 atoms when we use a cluster approach. The number is larger than that of the QM part but much smaller than that of the total atoms including the QM and MM parts in ref (17). Although the approach is simple, the reasonable results can be obtained when the solvent effect is taken into consideration. Our theoretical investigations and the previous computational studies demonstrate that the cluster approach can be applied to treat the complex enzyme reactions well.

Conclusions

Using quantum chemistry computations, M08HX/6-31G(d,p), and a cluster model, we present the reaction potential surface for the acylate reaction of β-lactamase with penicillin G; calculate the corresponding Gibbs free energies of the reactant, transition states, intermediates, and the product including the solvent effect that has an important influence on the structures and energies; and discuss the reaction mechanisms. There are two reaction paths that we study: one is a multi-step process, which can easily occur, and the energy barrier of the rate-limit step is 19.1 kcal/mol; the other is a one-step process, which is difficult to react with a reaction barrier of 45.0 kcal/mol. The reason why the two paths have different barriers is explained.