AsiteDesign: a Semirational Algorithm for an Automated Enzyme Design

AsiteDesign: a Semirational Algorithm for an Automated Enzyme Design Sergi Roda, ⊥ Henrik Terholsen, ⊥ Jule Ruth Heike Meyer, Albert Canellas-Solé, Victor Guallar,* Uwe Bornscheuer,*and Masoud Kazemi* Cite This: J. Phys. Chem. B 2023, 127, 2661−2670 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: With advances in protein structure predictions, the number of available high-quality structures has increased dramatically. In light of these advances, structure-based enzyme engineering is expected to become increasingly important for optimizing biocatalysts for industrial processes. Here, we present AsiteDesign, a Monte Carlo-based protocol for structure-based engineering of active sites. AsiteDesign provides a framework for introducing new catalytic residues in a given binding pocket to either create a new catalytic activity or alter the existing one. AsiteDesign is implemented using pyRosetta and incorporates enhanced sampling techniques to efficiently explore the search space. The protocol was tested by designing an alternative catalytic triad in the active site of Pseudomonas fluorescens esterase (PFE). The designed variant was experimentally verified to be active, demonstrating that AsiteDesign can find alternative catalytic triads. Additionally, the AsiteDesign protocol was employed to enhance the hydrolysis of a bulky chiral substrate (1-phenyl-2-pentyl acetate) by PFE. The experimental verification of the designed variants demonstrated that F158L/F198A and F125A/F158L mutations increased the hydrolysis of 1-phenyl-2-pentyl acetate from 8.9 to 66.7 and 23.4%, respectively, and reversed the enantioselectivity of the enzyme from (R) to (S)-enantiopreference, with 32 and 55% enantiomeric excess (ee), respectively. ■INTRODUCTION Structure-based enzyme engineering is widely used in the development of biocatalysts for industrial purposes. 1−5 These approaches have been employed to engineer protein thermostability, enzyme activity, or substrate selectivity. FRESCO, 6 FoldX, 6,7 and Rosetta 8 are among the noteworthy methods for enhancing enzyme thermostability, but many more also have been developed. 9−15 Similarly, a wide range of computational methods are available for engineering enzyme activity or for changing the selectivity, for instance, Rosetta, 8 nAPOLI, 16 EnzymeMiner, 17 HotSpot Wizard, 18 Caver, 19 FireProt-ASR, 20 LoopGrafter, 21 and DaReUS-Loop. 22 However, Rosetta and its derived methods 23,24 have become one of the most widely used tools in this area. Current advances in deep learning structural prediction methods increased the number of available protein structures dramatically. 25,26 Additionally, de novo protein design techniques are becoming increasingly accurate, 24,27,28 which makes it possible to create tailor-made protein scaffolds. Combining statistical energy derived from homologous sequences with physics-based simulations is also shown to be a powerful approach for uncovering the origin of enzyme stability and catalysis. 29 Considering these developments, structure-based enzyme engineering is expected to play an even more critical role in optimizing enzymes for industrial applications. Structure-based enzyme engineering could be used to either optimize enzyme activity and selectivity for a given substrate or to introduce new functionality into a protein cavity. 30−36 The latter approach has the potential to create enzymes that are capable of catalyzing new chemical reactions. 37 For example, by mutating the catalytic glutamate, a glycosidase enzyme was converted to a glycosynthase. 34,35 Alternatively, new catalytic residues can be introduced in a protein cavity. This approach, for instance, was used to create a second active site with hydrolysis activity in a transaminase. The engineered multifunctional enzyme combines transaminase and hydrolase activities in a single protein scaffold, allowing the conversion of β-keto esters into β-amino acids, which can be used for the synthesis of a key precursor of a family of antidiabetic drugs. 33 These applications require (re)designing catalytic residues in a given binding pocket. One way to approach such tasks is by Received: October 8, 2022 Revised: February 15, 2023 Published: March 21, 2023 Articlepubs.acs.org/JPCB © 2023 The Authors. Published by American Chemical Society 2661 https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 Downloaded via 5.67.84.229 on November 23, 2025 at 13:09:41 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. grafting an active site into a protein scaffold, which has been performed, e.g., by Rosetta Match. 38,39 This method is based on identifying a suitable cavity by searching many protein scaffolds. Here we present an alternative method, AsiteDesign, for this task. AsiteDesign is capable of identifying the best positions for a set of predefined catalytic residues in a given active site without the need for searching different protein scaffolds. The method is based on a Monte Carlo (MC) simulation and has been implemented using pyRosetta. Additionally, the protocol employs enhanced sampling techniques to improve the simulation convergence. It also includes the sampling of rotatable substrate bonds, which can potentially improve the identification of design solutions that may not be found otherwise. Furthermore, the substrate sampling can easily be restricted to predefined intervals, allowing to perform partial optimizations for large molecules. Here, we demonstrate the application of the method by designing an alternative catalytic triad in the active site of the Pseudomonas fluorescens esterase (PFE). 40 Since enantioselective hydrolases enable the synthesis of chiral building blocks for drugs by the (dynamic) kinetic resolution of racemic esters, 5 the protocol was further employed to alter the enzyme selectivity for a bulky chiral substrate. ■COMPUTATIONAL METHODS AsiteDesign employs MC sampling to explore both the active site and the ligand degrees of freedom. As mentioned above, the protocol was employed to engineer a new catalytic triad in PFE and also to redesign the enzyme binding pocket. The details of these design approaches are presented in the following. Design of Catalytic Residues. For a given enzymatic chemical reaction, the amino acid identity of the catalytic residues is known a priori. Considering PFE as an example, the catalytic triad consists of a nucleophile residue (Ser) that attacks the substrate, a general base (His) that accepts the proton, and an acidic residue (Asp/Glu) that activates the base. The residues of the catalytic triad are common to all serine hydrolases. As such, the main task here is to identify a set of positions in the binding pocket that can accommodate the catalytic residues at the correct distances relative to each other. The catalytic design protocol is composed of three stages (Figure 1A): (1) initial assignment, (2) Monte Carlo sampling, and (3) re-spawning of the next epoch by adaptive reinforcement learning protocol. In the initiation stage, the catalytic residues are assigned to the user-defined positions randomly. The initiation stage is performed for each explorer independently. The Monte Carlo sampling stage involves the sampling of both the positions of catalytic residues and the ligand conformation. The active site and ligand sampling are performed in series. During the active site sampling, in each iteration, one catalytic residue is assigned to a new random position by mutating it, and the previously assigned position is mutated back to the wild-type (WT) amino acids (Figure 1B). The mutations are performed by the fast-relax algorithm of the rosetta library. 8 It should be highlighted that the previously assigned position can be mutated to any set of user-defined amino acids. In this work, the wild-type amino acids are used for recovery to minimize disturbing the active site. During this stage, the correct distances of the catalytic residues are enforced by imposing distance restraints. The ligand sampling is performed using a similar approach as PELE 41,42 by introducing random perturbations to the ligand rotation and translation degrees of freedom, as well as dihedral angles. The main difference in the AsiteDesign protocol is that the sampling of the ligand-rotatable bonds is performed by a partial minimization in a random angle interval, rather than by assigning a random value to the dihedral angle. That is, the ligand-rotatable bonds are chosen randomly and an interval in either a positive or negative direction is assigned to each selected dihedral angle at random (Figure 1C). The number of rotatable bonds that are selected and the magnitude Figure 1. Schematic representation of the main steps of the MC algorithm in the AsiteDesign (A). Sampling of the catalytic residues (B). Sampling of ligand-rotatable bonds using an iterative grid search (C). The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2662 of the angle intervals are defined by the user. The rotatable bonds that are not selected retain their current interval. Additionally, for a given dihedral angle, the selection of the new intervals is restricted to the immediate neighborhood of its current interval to avoid introducing large conformational changes to the ligand. For example, a user-defined new interval in the positive direction is assigned to X1 in Figure 1C, whereas X2 maintains its current interval. The ligand dihedral angles are then minimized in the assigned angle intervals (Figure 1C). The minimization is performed using a two-step process. In the first step, the ligand’s lowest-energy conformation is identified by sampling its rotatable bonds systematically in a grid constructed using the assigned intervals and evaluating the ligand energy. This step is repeated three times where, in each iteration, the resolution of the grid search is increased and the search interval is reduced. The final ligand conformation is then minimized with a gradient descent optimization in cartesian coordinates. This approach thus corresponds to a partial minimization in a randomly chosen sampling interval. Compared to the conventional method of choosing a random angle, this approach can identify the most suitable conformation for a given interval. The algorithm is fully customizable and individual rotatable bonds of the ligand can be excluded by the user to speed up the AsiteDesign simulation. Adequate sampling could be challenging as the number of design elements increases. The adaptive reinforcement learning protocol 43 was incorporated into the simulation to overcome this issue. In this scheme, the simulation is performed in epochs in which the MC sampling is performed by multiple explorers in parallel using a distributed memory parallelization scheme. At the end of each epoch, the results from explorers are collected and ranked based on a given objective function (total energy, ligand energy, restraint energy, etc.), and the next epoch is spawned by the top-ranking results (Figure 1A). To further improve the sampling, a simulated annealing scheme was also used during the simulation. While simulated annealing helps the trajectories to escape local minima, the explorers often converge to different solutions (mutations) depending on the conformation of the ligand. The adaptive reinforcement learning algorithm here is used to identify the solutions that are more promising and thereby to allocate computational resources to the part of the search space that is of interest. Binding Pocket Redesign. The redesign of the binding pocket follows a similar approach to the design of catalytic residues. The main difference here is that the mutation of the noncatalytic residues is selected by the fast-relax algorithm of the rosetta library using a set of user-defined amino acids. This corresponds to a minimization step, in which the amino acid with the most favorable potential energy is selected for a given position. The number of positions included in the fast-relax minimization is defined by the user, and during this step, the positions with the most unfavorable rosetta energies are subjected to minimization. Therefore, in this approach, the noncatalytic residues are minimized to improve the interaction between the ligand and active site residues. The ligand sampling is performed as described above. As such, the noncatalytic residues of the active site are evolved in response to the ligand conformation. To ensure that the software can run consistently across different environments and architectures, we have developed a container of AsiteDesign, allowing an easy installation and distribution of the software reliably and securely. The container installation can be done via a script based on an already existing container with all the packages and dependencies and installing PyRosetta with the user credentials. A license for PyRosetta is needed. Molecular Dynamics (MD). Four replicas of 100 ns of molecular dynamics (MD) simulations with OPENMM 44 were performed to analyze the stability of the newly designed catalytic triads. A water cubic box (distance of 8 Å between the closest protein atom and the edge of the box) was created around the system using the TIP3P water model, and the charge of the system was stabilized using monovalent ions (Na+and Cl−). The protein system was parameterized with the AMBER99SB force field. Andersen thermostat and MC barostat were applied for the NPT ensemble (constant pressure and temperature, being 1 bar and 300 K, respectively). An NVT equilibration phase lasted 400 ps using a constraint of 10 kcal/(mol·Å2) to the whole solute system, followed by a 1 ns NPT equilibration with a milder constraint of 5 kcal/(mol·Å2); the production run only included constraints between H and heavy atoms. The Verlet integrator with a 2 fs time step was used with an 8 Å nonbonded long-range interactions cutoff. Protein Energy Landscape Exploration Simulations. Protein Energy Landscape Exploration (PELE) was used to analyze the substrate binding of the evolved variants using AsiteDesign. PELE is a heuristic MC-based algorithm coupled with protein structure prediction methods. 41,42 The software begins by sampling the different microstates of the ligand through small rotations and translations. Applying normal modes through the anisotropic network model (ANM) approach, 45 the protein’s flexibility is also considered. Once the whole system has been perturbed, side chains of the residues close to the ligand are sampled to avoid steric clashes. Last, a truncated Newton minimization with the OPLS2005 force field is performed, 46 and the new microstate is accepted or rejected based on the Metropolis criterion. The Variable Dielectric Generalized Born Non-Polar (VDGBNP) implicit solvent model 47 was used to mimic the effect of water molecules around the protein. The exploration of the substrate was enhanced with Adaptive-PELE 43 to improve the exploration of the search space. 48 ■EXPERIMENTAL DESCRIPTION Material. The chemicals rac-1-phenyl-2-pentanol (≥99%) and rac-1-phenylethyl acetate (≥98%) were ordered from Sigma-Aldrich. All other chemicals and solvents were purchased from Sigma-Aldrich, VWR, or Carl Roth and were used without further treatment. Synthesis of 1-Phenyl-2-pentyl Acetate. Five hundred microliters of acetic anhydride and 100 μLrac-1-phenyl-2pentanol were added in a 1.5 mL tube. The reaction was started by adding 10 μL pyridine to the mixture. The reaction was shaken at 25 °C and 500 rpm until complete conversion was achieved. Samples of 2 μL were withdrawn and diluted in 198 μL ethyl acetate for gas chromatographic (GC) analysis. The reaction was quenched by adding the mixture to a 15 mL tube containing 2 mL ddH2O. The product rac-1-phenyl-2pentyl acetate (1) formed a second phase and was separated and dried over anhydrous sodium sulfate. The oily rac-1phenyl-2-pentyl acetate was obtained in 51% yield. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2663 Plasmid Construction and Site-Directed Mutagenesis. Synthetic genes of the PFE variants 1, 4, 8, 11, 12 in pET28a were ordered from BioCat (Heidelberg, Germany) using seamless cloning with the flanking regions 5′ aaggagatatacc 3′(5′flanking) and 5′CACCACCACCACCACCACTGAGATCCGG 3′(3′flanking). The mutants are based on the sequence of the PFE wild-type (GeneBank: WP_120448209.1). The sequences were extended by a Cterminal linker (GS) and a His6-tag, as is the case of the sequence used for the 1VA4 crystal structure. 49 PFE variants 2 and 3 were constructed based on PFE_1, PFE_5-7 were based on PFE_4, and PFE_9-10 were based on PFE_8 using the Q5 Site-Directed Mutagenesis Kit (New England Biolabs GmbH, Ipswich, U.K.). Nonoverlapping DNA-oligonucleotides were designed using the online NEBaseChanger tool for the mutations. The list of primers used for mutagenesis is given in Table S5. The annealing temperatures suggested by NEBaseChanger online tool (https://nebasechanger.neb.com/) were used for the polymerase chain reaction (PCR), which was performed according to the manufacturer’s protocol. The obtained constructs were amplified in Escherichia coli Top 10 and used for heat-shock transformation of E. coli BL21 (DE3). Protein Preparation. Precultures (4 mL Luria−Bertani (LB) containing kanamycin) of E. coli BL21 (DE3) colonies harboring the constructs for the expression of the PFE variants were incubated overnight (37 °C, 180 rpm). LB medium (containing 50 mg/L kanamycin) was inoculated with 1% (v/ v) of the preculture and incubated (37 °C, 180 rpm) until it reached an OD600 of 0.6. Protein expression was induced by the addition of isopropyl-β-D-thiogalactopyranoside (IPTG) to a final concentration of 0.5 mM followed by incubation for ∼20 h at 20 °C at 180 rpm. Cells were harvested by centrifugation at 10,000gand 4 °C for 3 min, and the cell pellets were resuspended with 4 mL equilibration buffer (50 mM potassium phosphate, 300 mM sodium chloride, 10 mM imidazole, pH 8.0). Cells were disrupted by sonication on ice (five cycles of 1 min sonication at 30% intensity, and 50% pulsed cycle) using a SONOPULS HD 2070 (BANDELIN Electronic GmbH & Co. KG, Berlin, Germany), and the lysates were clarified by centrifugation at 10,000gand 4 °C for 30 min. For purification, the crude lysates were applied to 1.5 mL Roti Garose-His/Ni Beads (Carl Roth, Karlsruhe, Germany). The resins were washed with 15 mL washing buffer (50 mM sodium phosphate, 300 mM sodium chloride, 20 mM imidazole, pH 8.0) before target proteins were eluted with elution buffer (50 mM sodium phosphate, 300 mM sodium chloride, 250 mM imidazole, pH 8.0). Protein-containing fractions were pooled and re-buffered in 50 mM KPi pH 7.5 using PD10 columns (GE Healthcare, Buckinghamshire, U.K.). PFE wild-type was expressed as previously reported. 50 Activity Assays. The activity of the PFE variants toward the hydrolysis of para-nitrophenyl acetate (pNPA) was analyzed. For this purpose, 20 μL of a 100 mM pNPA solution in dimethyl sulfoxide (DMSO) was added to a 96-well plate and 180 μL of a PFE solution of known concentration in 50 mM KPi pH 7.5 was added to start the reaction. The absorbance was followed at 405 nm, and the initial slope was calculated. The reaction was carried out at 25 °C. Autohydrolysis was determined by adding 180 μL of the 50 mM KPi pH 7.5 buffer and subtracting the value from the hydrolysis rate of the PFE mutants. Specific activity was calculated using a standard curve of para-nitrophenol (0− 1 mM). The hydrolysis of rac-1-phenyl-2-pentyl acetate (1) and rac-1-phenylethyl acetate (2) was analyzed for all PFE mutants. For this purpose, 975 μL of a 90 μg/mL PFE solution was added to a 1.5 mL tube. The reaction was started by adding 25 μL of a 200 mM substrate solution in acetonitrile (final concentration 5 mM) and was run for 24 h at 37 °C and 1000 rpm. Time samples of 200 μL were taken after 1, 2, 4, and 24 h and extracted with 200 μL of ethyl acetate twice. The organic phases were pooled, dried over anhydrous sodium sulfate, and analyzed by GC. Gas Chromatography (GC) Analysis. Analysis was performed by gas chromatography with a flame ionization detector (GC-2010, Shimadzu, Kyoto, Japan) equipped with a Hydrodex β3P column (25.0 m ×0.25 mm, 0.25 μm film thickness, Macherey−Nagel, Duren, Germany). For the detection of the synthesis or hydrolysis of 1column temperature was held at 95 °C for 30 min, increased to 110 °C with 5 °C/min, and held for 45 min. Retention times: (S)-1-phenyl-2-pentyl acetate 50.5 min, (R)-1-phenyl-2-pentyl acetate 51.5 min, (S)-1-phenyl-2-pentanol 66.7 min, (R)-1phenyl-2-pentanol 69.2 min. For the detection of the hydrolysis of 2column temperature was held at 110 °C for 30 min. Retention times: (S)-1-phenyl-ethyl acetate: 4.3 min, (R)-1-phenyl-ethyl acetate: 5.9 min, (S)-1-phenylethanol: 8.0 min, (R)-1-phenylethanol: 7.1 min. ■RESULTS AND DISCUSSION Catalytic Residue Redesign. The esterase I from P. fluorescens (PFE) was chosen as the test scaffold protein. PFE has been extensively studied by mutagenesis, and its activity has been characterized for multiple esters. 50−52 It has been shown previously that the wild-type (WT) enzyme and its studied variants exhibit low activity for substrate 1(Scheme 1), which makes it a good candidate for testing AsiteDesign. This Scheme 1. Kinetic Resolution of Substrates 1 and 2 Studied Using PFE and Its Variants The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2664 enzyme hydrolyzes small aliphatic esters, and its active site contains the typical Ser-His-Asp catalytic triad (Figure 2). To test the performance of AsiteDesign in identifying optimum positions for placing the catalytic residues, the amino acids forming the esterase catalytic triad in the WT enzyme (Ser94, His251, and Asp222) were mutated to Ala. Using this mutated structure, an MC search was performed, employing ethyl acetate as the probing substrate. In the simulation, all residues forming the first shell of the active site (Table S1) were allowed to be mutated to one of the residues of the catalytic triad. As mentioned in the Methods Section, during the simulation, once a new position is accepted for a given catalytic residue, the previous position of the catalytic residue is mutated back to the WT amino acid (Figure 1B). Encouragingly, the catalytic residues of the WT enzyme were recovered as the best solution (Table 1). This result demonstrates that the protocol can indeed identify the optimal positions for the catalytic residues. In addition to the native catalytic triad, the second-best variant contains a catalytic triad at positions S28, H29 and D191. It is interesting to highlight that, in this variant, the catalytic residues are the mirror image of the WT enzyme (Figure 3). Therefore, this variant is expected to exhibit opposite enantioselectivity relative to the WT enzyme. Interestingly, the simulation resulted in multiple variants (PFE_1, PFE_1*+ I155D, PFE_1*+ A183D) in which both positions 28 and 29 are assigned to the Ser and His residues, respectively (Table 1), indicating that the probability of sampling these mutations is high. The only difference between these variants is the location of the acid residue (Table 1 and Figure S1). To test the stability of the designed variants, 100 ns MD simulations were performed for the WT enzyme and PFE_1. These simulations indicated that the designed variant is less stable compared to the WT enzyme. The average distances of the catalytic triad of PFE_1 (Ser28-His29: 5.30 ±0.38 Å, His29-Asp191: 2.74 ±0.76 Å) were found to be higher than that of the WT enzyme (Ser94-His251: 2.86 ±0.27 Å, His251Asp222: 1.85 ±0.06 Å). Based on the visual inspection of the MD trajectories, variants PFE_1 + C194T (PFE_2) and PFE_1 + V195M (PFE_3) were created to improve the enzyme stability (MD simulation results in the Supporting Information (SI);Figures S2−S4). These variants, however, exhibited similar MD metrics as PFE_1. Although PFE_1 and its derivatives appear to be less stable than the WT enzyme, the catalytic distances of these variants were in acceptable ranges, and thus, they were chosen for the experimental characterization. The computationally designed variants, recombinantly expressed in E. coli and purified, were then verified experimentally to characterize the enzymes’ activities (Table 2). In these variants, the native catalytic machinery was disabled by mutating the nucleophilic Ser94 to Ala. Because the main objective of the experimental characterization was to test the activity of the identified alternative catalytic triad without the interference of extra mutations, His251 and Asp222 were not mutated. The experimental results showed that the designed variants exhibit hydrolysis activity and PFE_1 is indeed active in the hydrolysis of pNPA and the racemic compounds 1and 2(Table 2). Its activity, however, was lower than the WT enzyme. This could be due to the destabilization effect of mutations, as it can be seen from the decreased melting temperature of the designed variants (Table 2). Alternatively, the lower observed activity could be because of less optimum catalytic distances and less stability of the active site as indicated by the MD simulations. These observations imply that this variant has a less organized catalytic geometry. This phenomenon has been observed in other designed (or natural) hydrolase active sites, where improving these distances gave better overall activities. 31,53 Figure 2. PFE and its catalytic residues. The catalytic triad residues are colored in red and labeled (PDB code: 1VA4). 49 Table 1. Top 10 Catalytic Designs Given at the End of the AsiteDesign Simulation with PFE’s structure a total energy mutations −3170.1 A94S/A251H/A222D (WT) −3160.1 W28S/L29H/T191D (PFE_1) −3159.6 A94S/V225H/A222D −3157.8 A94S/A251H/F162D −3156.8 A94S/A251H/I224D −3156.1 W28S/L29H/I155D (PFE_1*+ I155D) −3150.8 A94S/V225H/F125D −3150.1 W28S/L29H/A183D (PFE_1*+ A183D) −3148.5 W28S/V195H/T191D −3147.7 W28S/M95H/V121D a The reported energies are rosetta potential energies. PFE_1*stands for W28S/L29H/S94A. Figure 3. PFE WT and the newly designed active site. The catalytic triad residues of the WT are colored in red, while the ones from the PFE_1 design are shown in yellow. The labels are based on 1VA4 structure. The figure displays that the designed variant is the mirror image of the WT active site in the same protein cavity. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2665 This was predicted by the MD simulations of the mutants. Interestingly, the PFE_1 variant exhibited an inverse selectivity for the bulky compound 1. Variants PFE_2 and PFE_3 were both inactive toward 1but showed activity toward pNPA. It is known that imidazole used for protein purification can also hydrolyze the reactive substrate pNPA. 54 The hydrolysis of substrate 2was analyzed to rule out the possibility that the detected activity in the pNPA assay was caused by imidazole impurities, which should no longer be present after protein purification since buffer exchange was performed. Substrate 2 is not imidazole hydrolyzable (data not shown) and less challenging PFE substrate than bulky substrate 1and therefore ideal to measure even low enzymatic hydrolysis activities. Since all purified PFE variants showed activity and selectivity in the hydrolysis of 2(Table 2), and autohydrolysis was not observed, the experimental data clearly confirm enzymatic hydrolysis. These results suggest that, for a given binding pocket, the protocol is able to identify multiple viable solutions for designing catalytic residues, which can be used as a starting point for further optimization. Binding Pocket Redesign. To test the performance of AsiteDesign for noncatalytic residues, 1was chosen as the substrate (Scheme 1). The WT enzyme exhibits low activity and enantioselectivity for 1, which makes it a good candidate for improvement. Additionally, the previous site-directed mutagenesis of this enzyme did not yield any variants with high activity for 1. 50 The binding pocket design simulations were performed by including 31 residues of the active site (Table S2; design domain, notice that catalytic residues were excluded). In these simulations, no assumptions were made for the positions of mutations and all residues that are present in the first shell of the active site, 11 residues (Figure 4; highlighted in yellow), were allowed to mutate while the rest were only repacked. Since the enzyme was expected to hydrolyze a hydrophobic substrate, the allowed mutations were limited to hydrophobic amino acids (A, I, L, F, P, W, V, Y). In addition, sequence restraints were imposed on all mutable residues (i.e., introduction of any non-WT amino acid is penalized) to avoid large divergence from the WT enzyme, thereby favoring sequence conservation with an energy penalty. The substrate was placed in the active site manually, and MC simulations were performed while imposing distance restraint between the carbonyl carbon of the substrate and Ser94. Two separate MC simulations were performed for the (R)- and (S)-enantiomers of 1, hence evolving the active site for each enantiomer independently. For each enantiomer, the 50 variants with the overall best energies (protein and substrate binding) were selected. These structures were then clustered based on the binding mode of the substrate, and from each cluster, one variant was chosen (Table S3). Encouragingly, the simulations targeted many of the active site positions that were previously suggested to be important for enantioselectivity (F125, F158, and I224) 50 in addition to some new positions (W28, V121, and F198). However, the predicted mutations for these positions may differ from the previous study. To test the predicted variants, substrate binding was simulated by PELE software. 42 In these simulations, the Table 2. Experimentally Measured Activities for the Catalytic Designs in the H ydrolysis of Substrates; pNPA, 1 and 2 c , d PFE variants mutations pNPA activity (U/mg) substrate 1substrate 2 Tm (°C) predicted selectivity WT 162.2 8.9% (8% ee (R), E 1) a 48.7% (91% ee (R), E 59) b 71.8 PFE_1 W28S/L29H/T191D/S94A 0.9 2.3% (13% ee (S), E 1) a 17.7% (80% ee (R), E 11) a 44.9 (S) PFE_2 W28S/L29H/T191D/S94A/C194T 0.2 not detectable 2.1% (43% ee (R), E 3) a 50.9 (S) PFE_3 W28S/L29H/T191D/S94A/V195M 1.3 not detectable 0.9% (11% ee (R), E 1) a 44.8 (S) a After 24 h. b After 1 h. c pNPA activity is reported according to specific activity, while for substrates 1and 2the conversion is reported. The residue numbering corresponds to the 1VA4 structure. The melting points (Tm) of the PFE variants were determined by nanodifferential scanning fluorimetry (NanoDSF). 55 d Evalues were calculated according to Chen et al. 56 Figure 4. Active site of PFE and the used design domain. The catalytic triad residues are colored in red, the mutable residues in yellow, and the only repackable residues in violet. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2666 substrate was placed outside the active site, and the binding was monitored by counting near-attack conformations (NAC) 57 and computing the average energy of the ligand in the active site (Figure 5). These simulations offer a qualitative measure to identify variants that yield a productive binding mode. The primary goal of this filtering step was to limit the number of variants to be tested experimentally. A NAC is defined as a conformation where the distance between the carbonyl C of the substrate and the alcoholic O of the catalytic Ser residue is within 4 Å, while the H-bonds of the catalytic triad are within reasonable distances (≤3.5 Å). 32 Also, the distance between the catalytic His residue and the ether O of the substrate is less than 6.5 Å (as the protonated His residue will give a proton to this atom later on in the reaction to release the alcohol product 58 ). All the thresholds used for identifying NACs are based on previous studies. 32,33,59 Overall, the predicted variants exhibited a higher number of NACs compared to the WT (Table S3). Moreover, the distribution of the average interaction energy (Figure S5) is better in many variants compared to the WT, and key catalytic distances have good values as well (Figures S6−S9). Based on the in silico analysis, eight variants with the highest number of NACs relative to the WT enzyme were selected for experimental verification (Tables 3 and S4). The experimental results show that three predicted variants (PFE_5, PFE_8, and PFE_10) exhibited significant improvement over the WT enzyme in the hydrolysis of 1and, in contrast to the WT enzyme, they are selective for the (S) substrate (Table 3). These variants contain F158L, F125A, and F198A substitutions. The active site cavity of PFE is surrounded by bulky residues that hinder the access of large substrates. Thus, the mutation of these voluminous residues to smaller hydrophobic ones opens up the active site to accommodate the bulky substrates (Figure 6). The results demonstrate that AsiteDesign not only identified the positions of residues that needed to be changed but also predicted beneficial mutations as well. The developed method resulted in double and triple mutants with enhanced activity toward substrate 1. Hence, AsiteDesign can both guide the engineering of active sites of enzymes by identifying the hot spots and suggest variants with enhanced activity toward the hydrolysis of a specific substrate. Additionally, the suggested variants can serve as a starting point for further optimization, e.g., by directed evolution. Directed evolution has proven to be a powerful tool for optimizing newly introduced activities in protein scaffolds. 60 However, the simulations were not able to predict the variants’ enantioselectivity accurately. The main reason for this could be that the design/predictions were performed based on Figure 5. Initial setup for PELE simulations and the representation of a NAC. The substrate is placed outside the active site and allowed to explore around the drawn box (top). The NAC is represented with every key distance highlighted in a different color (blue for serine− histidine, beige for acid−histidine, violet for serine−substrate, and green for histidine−substrate) (bottom). The catalytic triad residues are colored in red, and the substrate in cyan. Table 3. Experimental Measured Activities for the Binding Pocket Redesigns in the Hydrolysis of Substrate 1 a PFE variants mutations Tm(°C) substrate 1predicted selectivity WT 71.8 8.9% (8% ee (R), E 1) PFE_4 W28A/F158L/F198A 57.0 6.0% (3% ee (S), E 1) PFE_5 F158L/F198A 69.5 66.7% (32% ee (S), E 4) (S) PFE_6 W28A/F125A/F158L/F198A 58.2 3.3% (9% ee (S)) (S) PFE_7 W28A/F158L/F198A/I224L 57.3 5.2% (8% ee (R), E 1) (S) PFE_8 F125A/F158L 63.6 23.4% (55% ee (S), E 4) (R) PFE_9 F125A/F158L/I224L 62.0 1.4% (29% ee (R)) (R) PFE_10 F125A/F158L/F198A 63.6 16.3% (60% ee (S), E 4) (R) PFE_11 V121A/F125A/I224L 59.7 2.9% (38% ee (S)) (R) PFE_12 V121A/F158A/F198V 60.4 1.6% (100% ee (S)) (R) a The activity is reported as conversion after 24 h. The residue numbering corresponds to the 1VA4 structure. The “predicted selectivity” column is based on whether the mutant was obtained from the AsiteDesign simulation with (R)- or (S)-enantiomer. The melting points (Tm) of the PFE variants were determined by nanodifferential scanning fluorimetry (NanoDSF). 55 The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2667 ligand binding energies, which do not necessarily correlate with enantioselectivity. This can potentially be improved by incorporating a transition state analog as the probing substrate, which is the main deciding factor for enantioselectivity. Another option could be to measure the energy barrier of each enantiomer using relatively inexpensive quantum mechanics/molecular mechanics (QM/MM) methods such as the empirical valence bond method. 61 Nevertheless, an accurate prediction of the enantioselectivity is notoriously challenging as the energy difference between the activation energies of enantiomers is often very small. The binding pocket-designed variants also exhibited lower melting temperatures (Table 3). The main reason for this is that the simulations are driven by improving the ligand binding energies, at the expense of protein stability. This issue can be alleviated by downstream enzyme stability optimization of the designed variants either computationally or experimentally. 6−10 ■CONCLUSIONS This work presents the AsiteDesign protocol, which aims at engineering active sites of enzymes to either introduce new catalytic residues or modify an existing active site in silico. The protocol is implemented using the pyRosetta library and combines MC sampling of the active site residues with enhanced sampling techniques to identify the most suitable positions of catalytic residues for a given active site. The ligand sampling is also included in the simulation, which is necessary to determine the optimal solutions. To demonstrate the performance of the protocol, a new catalytic triad was designed in the active site of the esterase I from P. fluorescens (PFE). The experimental characterization demonstrated that the designed variants are not only active biocatalysts, but they also exhibited inverse enantioselectivity for the bulky chiral substrate 1. Thus, the binding pocket of the enzyme was also successfully engineered to improve the activity for 1. Overall, these examples demonstrate that the AsiteDesign protocol can identify multiple viable solutions for designing active site residues for a given active site. This approach, thus, can be used in the engineering of multifunctional catalysts or in designing new catalytic residues in a given putative binding pocket. ■ASSOCIATED CONTENT Data Availability Statement The code used in the study is available at https://github.com/ masoudk/AsiteDesign. The code is also available as a container at https://github.com/BSC-CNS-EAPM/AsiteDesigncontainer. * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcb.2c07091. Recompilation of all the residues in the design domain of the catalytic residues redesign’s experiment; recompilation of all the residues in the design domain of the binding redesign pocket’s experiment; PELE simulation results of selected mutants with 1; experimental measured activities for the binding pocket redesigns with the other tested substrates; primers used for the construction of the designed variants; the potential positions for the design of an alternative catalytic triad in the PFE esterase; short description of the results from the MD simulations followed by the violin plots of the distribution of different metrics of interest; violin plots of the distribution of different metrics of interest obtained in the PELE simulations (PDF) ■AUTHOR INFORMATION Corresponding Authors Victor Guallar −Barcelona Supercomputing Center (BSC), Barcelona 08034, Spain; InstitucióCatalana de Recerca i Estudis Avancats (ICREA), Barcelona 08010, Spain; orcid.org/0000-0002-4580-1114; Email: victor.guallar@ bsc.es Uwe Bornscheuer −Department of Biotechnology &Enzyme Catalysis, Institute of Biochemistry, University of Greifswald, D-17487 Greifswald, Germany; orcid.org/0000-00030685-2696; Email: [email protected] Masoud Kazemi −Barcelona Supercomputing Center (BSC), Barcelona 08034, Spain; Biomatter Designs, Vilnius 09120, Lithuania; orcid.org/0000-0002-0750-8865; Email: [email protected] Authors Sergi Roda −Barcelona Supercomputing Center (BSC), Barcelona 08034, Spain; orcid.org/0000-0002-01747435 Henrik Terholsen −Department of Biotechnology &Enzyme Catalysis, Institute of Biochemistry, University of Greifswald, D-17487 Greifswald, Germany Jule Ruth Heike Meyer −Department of Biotechnology & Enzyme Catalysis, Institute of Biochemistry, University of Greifswald, D-17487 Greifswald, Germany; orcid.org/ 0000-0003-0293-4686 Albert Canellas-Solé −Barcelona Supercomputing Center (BSC), Barcelona 08034, Spain Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jpcb.2c07091 Author Contributions ⊥ S.R. and H.T. contributed equally to this work. Figure 6. Representative catalytic pose of the WT enzyme and the successful in silico evolved variants. The catalytic triad residues are colored in red, the substrate in cyan, and the mutated residues in yellow. Possible π−πinteractions between the substrate and Phe residues are shown with a dashed green line. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2668 Notes The authors declare no competing financial interest. ■ACKNOWLEDGMENTS The authors thank Barcelona Supercomputing Center (BSC) for providing the computational resources. S.R. thanks the Spanish Ministry of Science and Innovation for Ph.D. fellowship (FPU19/00608). S.R. and V.G. were funded by the FuturEnzyme Project of the European Union’s Horizon 2020 Research and Innovation Program (Grant Agreement No. 101000327). M.K. thanks Juan De La Cierva-Formación for their support (FJC1018-038089). M.K., V.G., and A.C.-S. received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement 101000607 (Project OXIPRO). H.T. was funded by the Leibniz Association’s strategic networking funding program Leibniz ScienceCampus ComBioCat. ■REFERENCES (1) Bornscheuer, U. T.; Huisman, G. W.; Kazlauskas, R. J.; Lutz, S.; Moore, J. C.; Robins, K. Engineering the Third Wave of Biocatalysis. Nature 2012,485, 185−194. (2) Yi, D.; Bayer, T.; Badenhorst, C. P. S.; Wu, S.; Doerr, M.; Höhne, M.; Bornscheuer, U. T. Recent Trends in Biocatalysis. Chem. Soc. Rev. 2021,50, 8003−8049. (3) Lovelock, S. L.; Crawshaw, R.; Basler, S.; Levy, C.; Baker, D.; Hilvert, D.; Green, A. P. The Road to Fully Programmable Protein Catalysis. Nature 2022,606, 49−58. (4) Bell, E. L.; Finnigan, W.; France, S. P.; Green, A. P.; Hayes, M. A.; Hepworth, L. J.; Lovelock, S. L.; Niikura, H.; Osuna, S.; Romero, E.; et al. Biocatalysis. Nat. Rev. Methods Primers 2021,1, No. 46. (5) Wu, S.; Snajdrova, R.; Moore, J. C.; Baldenius, K.; Bornscheuer, U. T. Biocatalysis: Enzymatic Synthesis for Industrial Applications. Angew. Chem., Int. Ed. 2021,60, 88−119. (6) Wijma, H. J.; Floor, R. J.; Jekel, P. A.; Baker, D.; Marrink, S. J.; Janssen, D. B. Computationally Designed Libraries for Rapid Enzyme Stabilization. Protein Eng. Des. Sel. 2014,27, 49−58. (7) Schymkowitz, J.; Borg, J.; Stricher, F.; Nys, R.; Rousseau, F.; Serrano, L. The FoldX Web Server: An Online Force Field. Nucleic Acids Res. 2005,33, W382−W388. (8) Rohl, C. A.; Strauss, C. E. M.; Misura, K. M. S.; Baker, D.Protein Structure Prediction Using Rosetta. Methods in Enzymology; Elsevier B.V., 2004; Vol. 383, pp 66−93. (9) Marques, S. M.; Planas-Iglesias, J.; Damborsky, J. Web-Based Tools for Computational Enzyme Design. Curr. Opin. Struct. Biol. 2021,69, 19−34. (10) Kulshreshtha, S.; Chaudhary, V.; Goswami, G. K.; Mathur, N. Computational Approaches for Predicting Mutant Protein Stability. J. Comput.-Aided Mol. Des. 2016,30, 401−412. (11) Musil, M.; Konegger, H.; Hon, J.; Bednar, D.; Damborsky, J. Computational Design of Stable and Soluble Biocatalysts. ACS Catal. 2019,9, 1033−1054. (12) Goldenzweig, A.; Fleishman, S. J. Principles of Protein Stability and Their Application in Computational Design. Annu. Rev. Biochem. 2018,87, 105−129. (13) Shirke, A. N.; Basore, D.; Butterfoss, G. L.; Bonneau, R.; Bystroff, C.; Gross, R. A. Toward Rational Thermostabilization of Aspergillus oyzae Cutinase: Insights into Catalytic and Structural Stability. Proteins 2016,84, 60−72. (14) Nguyen, V.; Wilson, C.; Hoemberger, M.; Stiller, J. B.; Agafonov, R. V.; Kutter, S.; English, J.; Theobald, D. L.; Kern, D. Evolutionary Drivers of Thermoadaptation in Enzyme Catalysis. Science 2017,355, 289−294. (15) Bednar, D.; Beerens, K.; Sebestova, E.; Bendl, J.; Khare, S.; Chaloupkova, R.; Prokop, Z.; Brezovsky, J.; Baker, D.; Damborsky, J. FireProt: Energyand Evolution-Based Computational Design of Thermostable Multiple-Point Mutants. PLoS Comput. Biol. 2015,11, No. e1004556. (16) Fassio, A. V.; Santos, L. H.; Silveira, S. A.; Ferreira, R. S.; de Melo-Minardi, R. C. nAPOLI: A Graph-Based Strategy to Detect and Visualize Conserved Protein-Ligand Interactions in Large-Scale. IEEE/ACM Trans. Comput. Biol. Bioinf. 2020,17, 1317−1328. (17) Hon, J.; Borko, S.; Stourac, J.; Prokop, Z.; Zendulka, J.; Bednar, D.; Martinek, T.; Damborsky, J. EnzymeMiner: Automated Mining of Soluble Enzymes with Diverse Structures, Catalytic Properties and Stabilities. Nucleic Acids Res. 2020,48, W104−W109. (18) Sumbalova, L.; Stourac, J.; Martinek, T.; Bednar, D.; Damborsky, J. HotSpot Wizard 3.0: Web Server for Automated Design of Mutations and Smart Libraries Based on Sequence Input Information. Nucleic Acids Res. 2018,46, W356−W362. (19) Pavelka, A.; Sebestova, E.; Kozlikova, B.; Brezovsky, J.; Sochor, J.; Damborsky, J. CAVER: Algorithms for Analyzing Dynamics of Tunnels in Macromolecules. IEEE/ACM Trans. Comput. Biol. Bioinf. 2016,13, 505−517. (20) Musil, M.; Khan, R. T.; Beier, A.; Stourac, J.; Konegger, H.; Damborsky, J.; Bednar, D. FireProtASR: A Web Server for Fully Automated Ancestral Sequence Reconstruction. Briefings Bioinf. 2021, 22, No. bbaa337. (21) Planas-Iglesias, J.; Opaleny, F.; Ulbrich, P.; Stourac, J.; Sanusi, Z.; Pinto, G. P.; Schenkmayerova, A.; Byska, J.; Damborsky, J.; Kozlikova, B.; Bednar, D. LoopGrafter: A Web Tool for Transplanting Dynamical Loops for Protein Engineering. Nucleic Acids Res. 2022,50, W465−W473. (22) Karami, Y.; Rey, J.; Postic, G.; Murail, S.; Tufféry, P.; de Vries, S. J. DaReUS-Loop: A Web Server to Model Multiple Loops in Homology Models. Nucleic Acids Res. 2019,47, W423−W428. (23) Khersonsky, O.; Lipsh, R.; Avizemer, Z.; Ashani, Y.; Goldsmith, M.; Leader, H.; Dym, O.; Rogotner, S.; Trudeau, D. L.; Prilusky, J.; et al.et al. Automated Design of Efficient and Functionally Diverse Enzyme Repertoires. Mol. Cell 2018,72, 178e5−186e5. (24) Anishchenko, I.; Pellock, S. J.; Chidyausiku, T. M.; Ramelot, T. A.; Ovchinnikov, S.; Hao, J.; Bafna, K.; Norn, C.; Kang, A.; Bera, A. K.; et al.et al. De Novo Protein Design by Deep Network Hallucination. Nature 2021,600, 547−552. (25) Jumper, J.; Evans, R.; Pritzel, A.; Green, T.; Figurnov, M.; Ronneberger, O.; Tunyasuvunakool, K.; Bates, R.; Zídek, A.; Potapenko, A.; et al.et al. Highly Accurate Protein Structure Prediction with AlphaFold. Nature 2021,596, 583−589. (26) Baek, M.; DiMaio, F.; Anishchenko, I.; Dauparas, J.; Ovchinnikov, S.; Lee, G. R.; Wang, J.; Cong, Q.; Kinch, L. N.; Schaeffer, R. D.; et al.et al. Accurate Prediction of Protein Structures and Interactions Using a Three-Track Neural Network. Science 2021, 373, 871−876. (27) Wang, J.; Lisanza, S.; Juergens, D.; Tischer, D.; Watson, J. L.; Castro, K. M.; Ragotte, R.; Saragovi, A.; Milles, L. F.; Baek, M.; et al.et al. Scaffolding Protein Functional Sites Using Deep Learning. Science 2022,377, 387−394. (28) Dauparas, J.; Anishchenko, I.; Bennett, N.; Bai, H.; Ragotte, R. J.; Milles, L. F.; Wicky, B. I. M.; Courbet, A.; de Haas, R. J.; Bethel, N.; et al. Robust Deep Learning-Based Protein Sequence Design Using ProteinMPNN. Science 2022,378, 49−56. (29) Xie, W. J.; Asadi, M.; Warshel, A. Enhancing Computational Enzyme Design by a Maximum Entropy Strategy. Proc. Natl. Acad. Sci. U.S.A. 2022,119, No. e2122355119. (30) Roda, S.; Santiago, G.; Guallar, V.Mapping Enzyme-Substrate Interactions: Its Potential to Study the Mechanism of Enzymes. Advances in Protein Chemistry and Structural Biology; Elsevier B.V., 2020; Vol. 122, pp 1−31. (31) Alonso, S.; Santiago, G.; Cea-Rama, I.; Fernandez-Lopez, L.; Coscolín, C.; Modregger, J.; Ressmann, A. K.; Martínez-Martínez, M.; Marrero, H.; Bargiela, R.; et al.et al. Genetically Engineered Proteins with Two Active Sites for Enhanced Biocatalysis and Synergistic Chemoand Biocatalysis. Nat. Catal. 2020,3, 319−328. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article https://doi.org/10.1021/acs.jpcb.2c07091 J. Phys. Chem. B 2023, 127, 2661−2670 2669