Long-Range Electrostatics in Serine Proteases: Machine Learning-Driven Reaction Sampling Yields Insights for Enzyme Design

Abstract

Computational enzyme design is a promising technique for producing novel enzymes for industrial and clinical needs. A key challenge that this technique faces is to consistently achieve the desired activity. Fundamental studies of natural enzymes revealed critical contributions from second-shell - and even more distant - residues to their remarkable efficiency. In particular, such residues organize the internal electrostatic field to promote the reaction. Engineering such fields computationally proved to be a promising strategy, which, however, has some limitations. Charged residues necessarily form specific patterns of local interactions that may be exploited for structural integrity. As a result, it is impossible to probe the electrostatic field alone by substituting amino acids. We hypothesize that an approach that isolates the influences of residues' charges from other influences could yield deeper insights. We use molecular modeling with AI-enhanced QM/MM reaction sampling to implement such an approach and apply it to a model serine protease subtilisin. We find that the negative charge 8 Å away from the catalytic site is crucial to achieving the enzyme's catalytic efficiency, contributing more than 2 kcal/mol to lowering the barrier. In contrast, a positive charge from the second-closest charged residue opposes the efficiency of the reaction by raising the barrier by 0.8 kcal/mol. This result invites discussion into the role of this residue and trade-offs that might have taken place in the evolution of such enzymes. Our approach is transferable and can help investigate the evolution of electrostatic preorganization in other enzymes. We believe that the study and engineering of electrostatic fields in enzymes is a promising direction to advance both fundamental and applied enzymology and lead to the design of new powerful biocatalysts.

Figure 1

Conserved charged residues in subtilisin-like enzymes and their identity and interactions in subtilisin Carlsberg. (A) Bioinformatic analysis yields conserved locations of charged groups in structurally similar enzymes. (B) Conserved charge centers in subtilisin Carlsberg closest to the catalytic triad correspond to residues Asp60 and Lys94. (C) An overview of the enzyme’s active site. Asp60 and Lys94 (yellow) participate in extensive h-bonding networks. The catalytic triad is colored orange, and a substrate peptide is shown in blue. Nonpolar hydrogen atoms are omitted for clarity. (D) Conservation of positions 60 and 94 in the S8A family. Position numbering corresponding to PDB entry 1R0R is employed herein and throughout the paper.

Figure 2

DeepTDA-powered reaction sampling quantifies the influence of conserved remote charged residues closest to the active site of subtilisin Carlsberg. (A) Reaction mechanism of the first (acylation) stage of the reaction catalyzed by subtilisin. (B) Reaction mechanism of the second (deacylation) stage. (C) Free energy profiles of the acylation stage in the intact model and with partial charges for Asp60 and Lys94 side chains shifted to yield net zero charge (Table S1). (D) Free energy profiles of the deacylation stage for the same models as in (C). Standard errors of means are shown but are often smaller than line widths. For their values refer to Figure S1.

Figure 3

Physical meaning and TS nondegeneracy of the learned DeepTDA variables. Shown are reweighted profiles for the intact enzyme. (A–C) First step of the acylation reaction stage. (D,E) Second step of the acylation reaction stage. (G–I) First step of the deacylation reaction stage. (J–L) Second step of the deacylation reaction stage. Ser = catalytic Ser221, His = catalytic His64, wat = attacking water molecule, and pept = substrate peptide. C (pept) = carbonyl oxygen of the substrate peptide’s Leu (P1′) residue, N(pept) = peptide group nitrogen of the substrate peptide’s Ala (P1) residue. State labels correspond to those from Figure 2A,B.

Figure 4

Disruption in the h-bonding network centered at the residue 60 affects subtilisin function through changes in the geometry of the catalytic Asp32. (A) Catalytic Asp32 geometry and h-bonding changes when the net charge on Asp60, but not on Lys94, is zero. (Upper graph) Distributions of the side chain rotation angles chi1 (N-CA-CB-CG) and chi2 (CA-CB-CG-OD1). (Lower graph) Distributions of the distance from Thr33 polar side chain hydrogen (HG1) to two potential h-bond acceptors (Asp32 OD1 and Asp60 OD2). Shown are distributions from the last 100 ns of each replica of enzyme–substrate complex simulations. (B) Original interaction geometry observed in the intact enzyme–substrate complex. Thr33 h-bonds Asp60 (yellow). (C) Altered interaction geometry was observed with the modified net charge for Asp60. Thr33 h-bonds catalytic Asp32. As in (B), nonpolar hydrogen atoms are omitted for clarity, the catalytic triad is shown in orange, the substrate peptide in blue. (D) Influence of the Thr33-Asp32 h-bond on the reaction efficiency. Profiles for systems with the presence of this h-bond are shown as dashed. (Upper graph) Free energy profiles for the acylation stage. (Lower graph) Free energy profiles for the deacylation stage. Standard errors of means are shown, but are often smaller than line widths.

Figure 5

Disruptions in h-bonding networks around residues 60 and 94 affect the substrate-binding loop. (A) Deviations of backbone positions of loop regions that harbor Asp60 and Lys94 and interact with them. Substrate-binding region 100–103 (lower graph) is of particular interest since it forms a beta-sheet with the substrate (C) Aggregate RMSD values from the last 100 ns of each replica are shown. (B) Distribution of the width of the substrate-binding groove in free enzyme models. Shown are measures from the last 100 ns. The K94M system was modeled for 1.5 μs, others for 0.5 μs per replica. (C) Example of changes in the substrate-binding gorge in the K94M variant (shown in yellow) as compared to the WT (gray). Substrate peptide is shown for illustrative purposes and colored blue. Dashed lines correspond to the distance plotted in panel B.

Figure 6

Changes in the outer-shell electrostatics might influence the catalytic efficiency through (de)stabilization of the catalytic His64. Shown is the difference between free energy profiles computed for the system with the “neutral” Asp60 (red) or Lys94 (blue) and the profile for the unperturbed system (Figure 2C,D). Shaded gray areas represent regions of the CV that correspond to the protonated cationic form of catalytic His64.