The global health crisis posed by bacterial resistance to β-lactam antibiotics is largely due to β-lactamases, enzymes that hydrolyze these crucial drugs. β-lactamases are categorized into four classes (A, B, C, and D), with Class A serine-based β-lactamases (SβLs) being the most prevalent and concerning. TEM-1, a representative Class A SβL, is commonly found in Gram-negative bacteria. Extensive research has illuminated the hydrolysis mechanism, with Glu166 proposed as a general base in benzylpenicillin acylation. The widely accepted mechanism involves Ser70 attacking the β-lactam carbonyl, forming a tetrahedral intermediate. Proton transfer occurs to Glu166 via bridging water, and Lys73 activates Ser130 to protonate the β-lactam nitrogen, cleaving the β-lactam bond and completing acylation. Other residues like Asn170 and Ser235 are also crucial for Michaelis complex formation. Computational methods, particularly QM/MM and molecular dynamics (MD), have deepened our understanding of the TEM-1 catalytic mechanism. However, challenges remain in addressing reactions with high degrees of freedom, as a single pathway might not fully represent the reaction. Free energy simulations have been successfully applied to similar systems, and chain-of-states (CoS) pathway optimization offers a cost-effective approach for finding minimum energy pathways (MEPs). This involves re-optimizing pathways under various MM potential fields sampled from MD simulations. Machine learning methods provide a solution to analyze the resulting massive data and complex correlations between geometrical degrees of freedom. These techniques have been used successfully in protein analysis, drug discovery, and accelerating QM/MM calculations. This study utilizes machine-learning regression to predict reaction pathway energetic profiles, quantifying the importance of chemical properties associated with structural descriptors.

Prior computational studies using QM/MM methods have investigated the TEM-1 catalytic mechanism, focusing on the acylation reaction with benzylpenicillin. These studies provided reaction pathways and potential energy surfaces (PESs), but limitations persisted in adequately representing the system's high degrees of freedom. Some studies concluded that the formation of the tetrahedral intermediate is the rate-limiting step, while others reported different findings. For example, Hermann et al. (2005) reported the intermediate as lower in energy, while Meroueh et al. (2005) indicated it as higher in energy compared to the reactant. These discrepancies highlight the need for a more comprehensive approach that considers multiple pathways and accounts for the complexities of the enzymatic environment. The role of key residues, especially Glu166, Lys73, and Ser130, has been extensively studied experimentally and computationally. However, there is still debate regarding the concertedness of proton transfer steps and the relative energetic importance of individual steps.

This study employs a combined QM/MM and machine-learning approach to analyze the TEM-1/benzylpenicillin acylation reaction. First, multiple minimum energy pathways (MEPs) were generated using the chain-of-states (CoS) method. This involved sampling various protein environments through constrained molecular dynamics simulations. Three configurations from an initial pathway were used as starting points for independent 200 ns MD simulations, during which the QM atoms were fixed, and the MM atoms were allowed to move freely. Eighteen representative conformations were selected from the trajectories for subsequent QM/MM geometry optimizations. Two QM levels of theory (DFTB3/mio and B3LYP/6-31G) were used for geometry and pathway optimizations. Single-point energies were further refined with larger basis sets (6-31+G*, 6-31++G**, 6-311++G**) and a dispersion-corrected B3LYP functional (B3LYP-D3), resulting in seven QM levels in total. Machine-learning regression models were then trained on these optimized pathways to predict the potential energy surface (PES). The input features were 105 pairwise distances between bonded atoms within the QM region. Three non-linear regression models—support vector regression (SVR), Gaussian process regression (GPR), and kernel ridge regression (KRR)—were employed. Feature selection was performed using recursive feature elimination (RFE) to identify the most relevant features for prediction. The prediction quality was assessed using root mean squared error (RMSE). To understand the energetic contributions of individual chemical events, two model-independent criteria were developed: intrinsic energy contribution and dynamic energy contribution. The intrinsic energy contribution quantifies the overall energetic contribution of a feature subset (representing a chemical event) to the reaction profile. This is measured as the RMSE difference between a model trained with the full feature set and a model with the target feature subset set to zero. The dynamic energy contribution examines the contribution of each chemical event along the reaction progress, using partial derivatives weighted by a factor that accounts for the variable space's domain size and feature correlation. This involved perturbing the input feature vectors and assessing the change in predicted energies.

The study yielded several key findings: 1. **Pathway Optimization and Prediction:** The CoS method generated 18 reaction pathways under various MM environments. Machine learning models accurately predicted these pathways, with RMSE values below 2.0 kcal/mol for B3LYP pathways. The DFTB3 pathways showed lower prediction accuracy due to more flexible and diverse conformational changes along the reaction coordinate. 2. **Tetrahedral Intermediate:** Unlike previous studies, the B3LYP optimized pathways showed that the tetrahedral intermediate is not always lower in energy than the reactant. This highlights the impact of the chosen QM methodology and conformational sampling. DFTB3 showed lower intermediates. 3. **Rate-limiting Step:** Contrary to previous conclusions, the analysis revealed that the rate-limiting step is not the formation of the tetrahedral intermediate. Both intrinsic and dynamic energy contributions showed that the concerted dual proton transfer from Lys73 to the β-lactam nitrogen via Ser130, along with β-lactam bond cleavage, are the most crucial events during acylation. 4. **Intrinsic Energy Contribution:** Analysis showed proton transfer between Lys73 and Ser130, and protonation of the thiazolidine nitrogen as the most energetically significant. Glu166 proton acceptance was less crucial. This aligns with experimental mutagenesis studies showing that Lys73 mutations deactivate the enzyme, while Ser130 mutations decrease but don't eliminate activity. The importance of Glu166 is supported by the fact that its mutations can alter the enzyme's function. 5. **Dynamic Energy Contribution:** This analysis confirms the concerted nature of the proton transfer processes and identifies the key events during both tetrahedral intermediate formation and acyl-enzyme product formation. The rate-limiting events were the proton transfer from Lys73 to the β-lactam nitrogen and the β-lactam scissile bond cleavage. 6. **Acylation Mechanism:** The analysis supports a concerted four-proton transfer mechanism for the acylation reaction, challenging prior interpretations that highlighted separate steps.

The findings address the research question by providing a quantitative assessment of the energetic contributions of individual chemical events in the TEM-1/benzylpenicillin acylation reaction. The machine-learning approach provides a novel way to analyze high-dimensional reaction pathways, overcoming limitations of previous methods. The identification of the rate-limiting steps as concerted dual-proton transfer events refines our understanding of the mechanism, challenging previous interpretations that emphasized the tetrahedral intermediate formation. This study demonstrates the power of combining QM/MM simulations with machine learning techniques for understanding complex enzymatic reactions. The agreement between intrinsic and dynamic energy contributions adds robustness to the conclusions, supporting the identification of the key chemical events. The findings are relevant to the field of drug design, offering insights that can inform the development of new β-lactamase inhibitors.

This study introduces novel machine-learning-based models for quantifying the energetic contributions and correlations of chemical events during enzyme catalysis. The model-independent criteria, intrinsic and dynamic energy contributions, successfully identified the rate-limiting steps in the TEM-1/benzylpenicillin acylation reaction. Future research could extend this methodology to other enzymes, explore different machine-learning models and feature selection techniques, and incorporate entropy effects for a more complete understanding of enzymatic catalysis.

While the study provides a comprehensive analysis of the TEM-1/benzylpenicillin acylation reaction, several limitations exist. The accuracy of the machine-learning predictions depends on the quality of the QM/MM data. The choice of QM level of theory affects the results, and a higher level might yield more precise data but increase computational cost. The current model doesn't explicitly include entropy effects, which can play a significant role in enzymatic reactions. The study focused on a specific enzyme-substrate pair; the findings might not be directly transferable to all β-lactamases or other enzyme systems.