Computational study on the polymerization reaction of D-aminopeptidase for the synthesis of D-peptides

Abstract

Almost all natural proteins are composed exclusively of L-amino acids, and this chirality influences their properties, functions, and selectivity. Proteases can recognize proteins composed of L-amino acids but display lower selectivity for their stereoisomers, D-amino acids. Taking this as an advantage, D-amino acids can be used to develop polypeptides or biobased materials with higher biostability. Chemoenzymatic peptide synthesis is a technique that uses proteases as biocatalysts to synthesize polypeptides, and D-stereospecific proteases can be used to synthesize polypeptides incorporating D-amino acids. However, engineered proteases with modified catalytic activities are required to allow the incorporation of D-amino acids with increased efficiency. To understand the stereospecificity presented by proteases and their involvement in polymerization reactions, we studied D-aminopeptidase. This enzyme displays the ability to efficiently synthesize poly D-alanine-based peptides under mild conditions. To elucidate the mechanisms involved in the unique specificity of D-aminopeptidase, we performed quantum mechanics/molecular mechanics simulations of its polymerization reaction and determined the energy barriers presented by the chiral substrates. The enzyme faces higher activation barriers for the acylation and aminolysis reactions with the L-stereoisomer than with the D-substrate (10.7 and 17.7 kcal mol⁻¹ higher, respectively). The simulation results suggest that changes in the interaction of the substrate with Asn155 influence the stereospecificity of the polymerization reaction.

---

Introduction

Of the 20 canonical proteinogenic amino acids, all except Gly have a chiral centre at the Cα backbone atom, resulting in the presence of L- and D-stereoisomers. Every living organism exclusively uses L-amino acids for constructing proteins. Although the biological occurrence of D-amino acids is rare, their presence has been reported in bacterial cell walls as components of peptidoglycan and other periplasmic extracellular polymers,¹ in antibacterial and antifungal peptides,^2,3 in the cellular fluids of some invertebrate marine worms and shellfish,⁴ in the venom of some spiders⁵ and platypuses^6,7 and in the skin secretions of some amphibians acting as homologs of mammal neurotransmitters and hormones.^4,8 The discovery of these peptides and their compelling biological functions has drawn increasing interest to the incorporation of D-amino acids into polypeptides.

The interest in polypeptides and peptides in various fields is increasing due to their diverse physical properties, biological functions and promising applications, such as in bio-based materials and as green alternatives to petroleum-derived plastics, pharmaceutical agents, therapeutic drugs, drug delivery agents, and antimicrobial peptides.^9–11 In particular, interest in peptides as delivery agents and inhibitors has arisen because they are highly specific for their targets and have high biocompatibility and low toxicity.^10,12 However, peptides have low biostability because they are prone to hydrolysis by proteases, reducing their half-life in circulating plasma, and they can be immunogenic, which severely limits their application as therapeutic agents in vivo.^12,13

Using peptides with D-amino acids in their structure is a highly promising strategy for overcoming these limitations. D-Amino acids can alter the secondary structure relative to that generate with all L-amino acids, which is the structure widely recognized by proteases. They can adopt different backbone φ and ψ angles that are disfavoured by L-amino acids, stabilizing specific structural conformations that cannot be assumed by L-amino acids. The different configurations available with D-amino acids can confer favourable biological and chemical proprieties to peptides and proteins. D-Amino acids can stabilize α-helices,¹⁴ β-turns and β-hairpins^15,16 and promote novel peptide topologies. The α-helices in proteins made entirely of L-amino acids are nearly always right-handed. On the other hand, it has been observed that proteins formed exclusively from D-amino acids are the mirror image of their all L-amino acid counterpart.^17,18 Thus, artificial proteins made of D-amino acids form left-handed α-helices, changing the conformation of the secondary structure.¹⁹ Therefore, D-peptides are highly resistant to proteolysis, increasing their biological half-life and biostability^20–23 and resulting in low immunogenicity, as they are not recognized by other proteins.²⁴ Furthermore, short D-peptides can be administered orally, and they are systemically absorbed; in contrast, L-peptides have to be injected to avoid digestion.²⁵ For these reasons, the incorporation of D-amino acids into peptides and proteins has many potentially benefits.

Although the incorporation of D-amino acids into polypeptides presents several advantages, neither chemical synthesis nor chemoenzymatic peptide synthesis have been efficient enough for their synthesis or incorporation into peptides. Chemical synthesis strategies face problems regarding the size of the peptides they can produce, as these strategies require multiple difficult and complex protection and deprotection steps²⁶ as well as toxic chemicals and solvents.²⁷ On the other hand, chemoenzymatic peptide synthesis is based on the use of protease enzymes as biocatalysts to synthesize polypeptides via the aminolysis of amino acids (Scheme 1). Under physiological conditions, proteases do not catalyze polymerizations since the thermodynamic equilibrium of the reactions strongly favours hydrolysis. Thus, in chemoenzymatic synthesis, polymerization is strongly dependent on fundamental reaction parameters, especially substrate concentration.²⁸ This technique affords high yields, is atom economical, can be performed under mild conditions, and does not require organic solvents.²⁹ However, one major drawback is that the synthesis is restricted by the natural substrate specificity of proteases and they generally only recognize L-amino acids. Thus, the use of proteases as biocatalysts is limited by their overall lack of D-stereospecificity and inability to recognize D-amino acids. As recently explored in our previous work with quantum mechanics/molecular mechanics (QM/MM) simulations,³⁰ the stereospecificity showed by papain, an enzyme that only polymerizes L-amino acids, was mainly influenced by the chiral conformation adopted by the acyl intermediate with the substrate (L- or D-Ala-OEt). The simulations revealed that the aminolysis of L-Ala-OEt by papain had lower energetic barriers than that of the reaction of D-Ala-OEt,³⁰ its energetic barriers for aminolysis and hydrolysis are comparable, and therefore, the reactions are competitive. However, the aminolysis reaction with D-acyl intermediates has higher energetic barriers, and because hydrolysis is much more energetically favourable, the polymerization of the D-amino acids into D-peptide does not occur.

Scheme 1 DAP-mediated polymerization showing the acylation, aminolysis and hydrolysis reactions.

Here, we investigate one of the few known D-stereospecific proteases, D-aminopeptidase (DAP) (EC 3.4.11.19, MEROPS Family S12), to develop chemoenzymatic polymerization reactions of D-amino acids. DAP is a serine protease isolated and purified from Achromobacter anthropic by Asano et al., and it has a molecular weight of ≈110 kDa and consists of 2 identical subunits.³¹ Based on the conserved elements in its amino acid sequence (the Ser–X–X–Lys and Tyr–X–Asn conserved motifs, where X can be any amino acid), its primary structure is similar to β-lactamases, R61 D-carboxypeptidase and penicillin-binding proteins, and its activity is inhibited by β-lactam antibiotics, and hence DAP is considered a penicillin-recognizing protein.^32–34 Indeed, its active site presents high structural homology with the active sites of class A and C β-lactamases and penicillin-binding proteins, and the residues Ser62, Lys65, Tyr153, Asn155, His287 and Gly289 are conserved. Among these residues, Ser62 and Lys65 are essential for catalytic activity as confirmed by site-directed mutagenesis,³¹ whereas Tyr153, located near these residues, stabilizes a hydrogen-bond network.^33,35 Regarding the peptidase activity, DAP presents remarkable stereospecificity towards amide bonds of D-amino acids, especially low-molecular-weight D-peptides, D-amino acid esters, and peptides with D-alanine at the N terminus.³¹ Due to this specificity for D-amino acids and because it is inhibited by β-lactam antibiotics, this enzyme may have a physiological role in the biosynthesis or degradation of peptidoglycans.^31,36,37 Additionally, DAP was reported to be able to synthesize D-alanine oligopeptides from D-alanine methyl esters in organic solvents.³⁸

Using DAP as a model enzyme of D-stereospecific proteases, we tried to identify the elements important to its stereospecificity. We hope that these results can be used to rationally modify the active sites of enzymes to increase their specificity and usefulness, which will facilitate the synthesis of more biostable peptides through the incorporation of D-amino acids, the creation of new peptide-derived biomaterials,^20,21 and the synthesis of fully D-peptides for testing as new pharmaceuticals.^13,39 Additionally, having access to D-stereospecific proteases could increase the number of accessible peptide structures for researchers, reduce their synthesis costs, and make peptides more easily available for exploration and use in new biotechnologies, such as in gene or protein delivery systems targeting cells or organelles,^40,41 in cell-penetrating peptides or in peptide-based antibiotics.^42–44

Experimental

Materials

L-Alanine ethyl ester hydrochloride (L-Ala-OEt) and D-alanine ethyl ester hydrochloride (D-Ala-OEt) were purchased from Sigma-Aldrich (St. Louis, MO, U.S.A.). L-Proline benzyl ester (Pro-OBzl) was purchased from Tokyo Chemical Industry (Japan). Deuterium oxide was purchased from Cambridge Isotope Laboratories, Inc. Standard chemicals were purchased from Wako Chemical Co. (Kanagawa, Japan).

Synthesis of D-alanine dipeptide ethyl ester as standard

To a dispersion of D-Ala-OEt (0.768 g, 5 mmol), N-tert-butoxycarbonyl-D-alanine (Boc-D-alanine) (0.946 g, 5 mmol), and 1-hydroxybenzotriazole monohydrate (HOBt, 0.766 g, 5 mmol) in chloroform (10 mL) in a flask equipped with an addition funnel and a stopcock was added triethylamine (0.70 mL, 5 mmol) at 0 °C under nitrogen (Fig. S1†). A solution of 1-ethyl-3-(3-dimethylaminopropyl)carbodiimide hydrochloride (0.959 g, 5 mmol) in chloroform (10 mL) was added dropwise to this mixture, and the solution was stirred at −10 °C for 30 minutes and then at 25 °C for 18 h. After the reaction, the solution was washed with 5% sodium hydrogen carbonate aq. solution and brine. The organic layer was dried with magnesium sulfate and concentrated with a rotary evaporator. After drying under vacuum, the Boc-protected dipeptide ester Boc-D-Ala-D-Ala-OEt was obtained as a white solid. Then, the solid was dissolved in dichloromethane (5 mL), and trifluoroacetic acid (1.9 mL, 25 mmol) was added at 0 °C under nitrogen, and the solution was stirred at 25 °C for 6 h. After the reaction, the solvent was removed by vacuum distillation. The crude viscous liquid was mixed with dioxane/HCl (4 M, 2.5 mL) and poured in excess diethyl ether. The white precipitate was separated by filtration, washed with diethyl ether, and dried under vacuum to give D-Ala-D-Ala-OEt as the hydrochloride salt. The yield was 0.67 g (60%).

D-Aminopeptidase expression and purification

D-Aminopeptidase from Ochrobactrum anthropi SCRC C1-38 cDNA was kindly provided by Prof. Asano. DAP was synthesized and purified using a reported previously cell-free protein expression system.⁴⁵ A cell-free dialysis was performed on a large scale (36 mL) using a dialysis membrane with a molecular weight cutoff (MWCO) of 15 kDa (Pierce, Rockford, IL) using the reaction conditions reported previously.⁴⁶ The internal solution was dialyzed in a dialysis tube (Spectra/Por 7, MWCO of 15 kDa, Spectrum) against the external solution at 30 °C for 16 h with shaking.⁴⁷ The internal solution with the tagged protein was purified by AKTA Xpress (GE Healthcare, Little Chalfont, U.K.) using a previously reported procedure.⁴⁶ Briefly, the protein solution was purified on a HisTrap column (5 mL, nickel–nitrilotriacetic acid (Ni–NTA) column, GE Healthcare). Sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) was performed using 15 to 20% precast Tris–HCl gels (DRC Co. Ltd, Kyoto, Japan). The gel was stained with Coomassie brilliant blue (Fig. S1c†).

Chemoenzymatic polymerization assays

The polymerization reactions with L- or D-Ala-OEt were based on previously reported reaction conditions.⁴⁸ The reactions were conducted at 25 °C in a final volume of 1 mL of 1 M PBS buffer (pH 8.0), containing 0.01 mg mL⁻¹ DAP and 0.5 M of L- or D-Ala-OEt monomer as the substrate. The reactions were allowed to proceed for 1 minute under shaking. In the control reactions, no enzyme was added to the assay. The reactions were immediately transferred to an ice-cold bath and centrifuged at 10 000g for 10 minutes in an Amicon 10K tube at 4 °C to remove the enzyme. The samples were analysed by RP-HPLC.

RP-HPLC analysis

RP-HPLC system consisted of an auto sampler (AS-2055, JASCO, Tokyo, Japan), a gradient pump (PU2089, JASCO) a column oven (CO-4060, JASCO), and a C18 column (YMC-Triart C18, particle size 5 μm, 150 × 4.6 mm i. d., YMC, Kyoto, Japan). The samples were injected into the HPLC system in the mobile phase, which was composed of 0.1 M PBS buffer (pH 8) as eluent A and acetonitrile as eluent B. The composition of the mobile phase was linearly changed from 98% A and 2% B to 90% A to 10% B over 30 minutes and was used at a flow rate of 1 mL min⁻¹. The elution of the various components was monitored by UV absorbance at 220 nm. The peak areas and retention times were quantified with the chromatography software (ChromNAV, JASCO, Tokyo, Japan) and corrected relative to an internal standard (L-proline benzyl ester (H-Pro-Bzl)). The concentration of the product was calculated by comparison to a calibration curve that was prepared by plotting the peak areas of the standard D-alanine dipeptide ethyl ester against the concentration.

Nuclear magnetic resonance

¹H NMR spectra were recorded on a Varian NMR 500 (500 MHz) spectrometer (Varian Medical Systems) at 25 °C and acquisition was controlled with the VnmrJ software. Sixty-four scans were taken during each NMR experiment. Data were processed and analyzed with ACD/NMR Processor Academic Edition, version 12.01 (Advanced Chemistry Development, Inc.).

Model preparation and classical molecular dynamics simulations

The crystallographic structure of DAP (PDB ID 1EI5)³³ was used for the simulations. The substrates L-Ala-OEt and D-Ala-OEt were first docked into DAP with the docking program AutoDock Vina version 1.1.2.⁴⁹ Before docking, all ligands and water molecules were removed from the PDB file. Polar hydrogens were added to ligands and receptor by using the Hydrogen module in AutoDock Tools version 1.5.6,⁵⁰ then Gasteiger united atom partial charges and atom types were assigned and PDBQT files generated. The protein molecule was kept rigid, while all the torsional bonds in the substrates were set free to rotate. A 15 Å docking box around the oxygen atom of Ser62 was defined. The best docking conformations of the substrates were used as initial configuration for the MD. The protonation state of the amino acid residues for the MD was set according to the environment and previous studies; for modelling the acylation reaction, Lys65 was set as neutral, and nearby residue His287 was kept neutral and ε-protonated.^35,51 All molecular dynamic simulations were performed using AMBER 16.⁵² The atoms in the protein were described using the AMBER ff14SB force field,⁵³ L- and D-Ala-OEt substrates were described with a GAFF,⁵⁴ and the system was solvated with TIP3P water molecules.⁵⁵ A time step of 2 fs along with the SHAKE algorithm was used,⁵⁶ while the particle-mesh Ewald (PME) method⁵⁷ was used to calculate long-range interactions. The systems were equilibrated by a series of gradual-step simulations in which distance restraints were set between the substrate molecules and the protein to maintain the important interactions. The simulations were extended until the system was considered equilibrated according the root mean square deviation (RMSD) (Fig. S2†). The relaxed structures obtained from the last part of the trajectories were used as the starting structures for hybrid QM/MM calculations. The trajectories were analysed using CPPTRAJ⁵⁸ and VMD software.⁵⁹

QM/MM molecular dynamics simulations

Hybrid QM/MM molecular dynamics simulations were carried out as previously described.³⁰ Briefly, we used density functional theory (DFT) B3LYP/6-31G* basis set⁶⁰ to describe the atoms in the QM system. The QM calculations were carried out using Gaussian09 (ref. 61) together with the AMBER/Gaussian interface.⁶² The integration time step for the QM calculations was 0.2 fs. The explicit link atom approach, implemented in sander, the molecular dynamics program of AMBER, was used to separate the QM and the MM regions when their boundary crossed covalent bonds.^63,64 The electrostatic interactions between the MM and QM regions were truncated at a cutoff of 8 Å.⁶²

QM/MM adaptively biased molecular dynamics

The free-energy landscapes (FELs) of the acylation and aminolysis reactions were determined using adaptively biased molecular dynamics (ABMD).^65,66 This enhanced sampling technique based on metadynamics^67,68 introduces a time-dependent biasing potential in the simulation, enabling the exploration of energy surfaces of selected collective variables (CV). In our study, the flooding timescale of the bias deposition was set at 30 fs (150 MD steps), and the resolution was set at 1 kcal mol⁻¹. The first crossing criterion was used to determine the end of the simulation, as recommended for chemical reactions.⁶⁹ Walls at 3.5 Å between substrates/products and Ser62 were used to reduce the FEL space of the chemical event and reduce computational effort.

QM/MM ABMD simulations of the acylation reaction

A total of 64 atoms (including 3 H link atoms) were treated quantum mechanically in the acylation reaction with the putative systems while 86 300 atoms were treated with MM. Three CVs (CV₁, CV₂ and CV₃) were used to describe the acylation reaction. CV₁ is the number of bonds between the O atom of Ser62 and the H atoms of O_Tyr153 and N_Lys65 (Fig. 1a). CV₂ was defined as the distance between the O1_Ala-OEt–H_Tyr153 atoms. CV₃ is defined as a linear combination of distances (LCOD) in the form of the distance between bond breaking minus bond forming in the reaction: the distance between the C1_Ala-OEt–O1_Ala-OEt atoms (a) minus the distance between the C1_Ala-OEt–O_Ser62 atoms (b) (CV₁ = d(a) − d(b)). CV₁ was selected because it samples the proton transfer between Ser62 to any nearby base. CV₂ was chosen to model the proton transfer from Tyr153 to O1 of Ala-OEt to form ethanol. CV₃ was selected to accelerate the nucleophilic attack of the Ser62 O to Ala-OEt and the breaking of the ester bond between the C1 and O1 atoms of Ala-OEt (Fig. 1a).

Fig. 1 Schematic of the atoms treated at the B3LYP/6-31G* QM level of theory. (a) Acylation reaction; the catalytic residues Ser62, Lys65 and Tyr153 together with the L- or D-Ala-OEt substrates were subjected to QM modelling. Green arrows indicate the atoms used in multiple bonds interaction in CV₁. The orange arrows indicate CV₂ and the blue arrows show the distances a and b in used as LCOD in CV₃. (b) Aminolysis reaction; all the represented atoms were subjected to QM modelling. CV₁ LCOD are indicated by blue arrows, the interactions present in CV₂ are shown with orange arrows and the atoms involved in interactions in CV₃ are indicated with green arrows. Link atoms are marked with curved lines.

QM/MM ABMD simulations of the aminolysis reaction

For the aminolysis reaction, a total of 94 atoms (including 5 H link atoms and a water molecule) were treated at the QM level of theory (B3LYP/6-31G* basis set) while 86 411 atoms with MM. The following CVs were used to model the reaction. CV₁ was the LCOD between the C1_Acyl-interm. and O_Ser62 atoms (a) minus the distance between the C1_Acyl-interm and N_Ala-OEt atoms (b) (CV₁ = d(a) − d(b)) (Fig. 1b). CV₂ was defined as the number of bonds between the N atom of the attacking Ala-OEt and the H atoms of N_Ala-OEt, O_Tyr153 and N_Lys65. CV₃ was defined as the number of bonds from O_Ser62 to the H atoms of N_Ala-OEt, O_Tyr153 and N_Lys65. CV₁ was chosen to sample the nucleophilic attack of the N atom of Ala-OEt to the C1 atom of the acyl intermediate and its cleavage from the O atom of Ser62. CV₂ was selected to model the proton transfer from one of the protons of the attacking nucleophile to any nearby group. CV₃ was selected to model the transfer of a proton to the O atom of Ser62 (Fig. 1b).

Results and discussion

DAP-mediated polymerization activity in aqueous buffer

DAP was synthesized using a 36 mL-scale cell-free protein synthesis method and was purified by affinity chromatography, resulting in a final protein yield of 48.3 mg (Fig. S1c†). We polymerized the peptides with DAP under mild conditions, namely, 1 M PBS buffer at pH 8.0 and 25 °C in the absence of any organic solvent based on our previous procedure for other enzymes.^29,48 Reactions were carried out with L- or D-Ala-OEt, and the soluble products of these reactions were quantified and analyzed by RP-HPLC. To identify the retention time of the product peak and quantify the rate of the reaction, D-alanine dipeptide ethyl ester was chemically synthesized and used as a standard (Fig. S1a†). The reactions were conducted using Boc-D-alanine and D-Ala-OEt in chloroform, and the residue was purified to afford a white powder in 60% yield. The synthesis of D-alanine dipeptide ethyl ester was confirmed by ¹H NMR (Fig. S1b†). D-Alanine dipeptide ethyl ester, together with L-Ala-OEt and D-Ala-OEt, were used to determine the retention time of each substance in the HPLC chromatograms (Fig. 2a). DAP is able to polymerize L- and D-Ala-OEt under aqueous conditions (1 M PBS, pH 8.0) without a tertiary amine or organic solvent (Fig. 2b).³⁸ Similarly, as reported previously for other enzymatic activities,^31,32 DAP displayed a higher polymerization efficiency for D-substrates (D-Ala-OEt ∼4000 U mg⁻¹) than for L-substrates (L-Ala-OEt activity ∼150 U mg⁻¹). In addition, additional peaks appeared in the HPLC chromatograms of the reactions with D-Ala-OEt, conceivably corresponding to the trimer and tetramer (Fig. 2b).

Fig. 2 Representative RP-HPLC chromatograms of the polymerization reaction. (a) RP-HPLC chromatograms of the standards used to characterize the reactions. The chromatograms of the L-Ala-OEt, D-Ala-OEt and D-alanine-OEt dipeptide standard are shown in black, red and blue lines, respectively. (1) The peak for L- or D-Ala-OEt (the substrate) at a retention time of 11–12 minutes. (2) The peak of the D-alanine dipeptide product appeared at a retention time of 16–17 minutes. (b) Representative RP-HPLC chromatograms of the reactions performed with DAP with L- or D-Ala-OEt, shown in black and red lines, respectively. (1) The peak of the L- or D-Ala-Ala-OEt products formed in the reaction. (2 and 3) Possible peaks of D-alanine ethyl ester tri- and tetrapeptides, respectively. (4) Peak of the internal standard used in the reaction.

QM/MM acylation reaction

To understand the unique stereospecificity of DAP and confirm its proposed mechanism,³⁵ we carried out QM/MM ABMD simulations since this simulation method is a powerful tool for elucidating the stereospecificity of protease enzymes in the polymerization of amino acids.³⁰ The equilibrated structures of DAP and the substrates were obtained after performing classical MD simulations. The catalytic triad is formed by the residues Ser62, Lys65 and Tyr153, which are in close proximity. The catalytic Ser62 oxygen forms a hydrogen bond with the amide group of Lys65, which in turn forms a hydrogen bond with Tyr153.^33,35 Moreover, based on QM studies with closely related enzymes of the same family with similar catalytic residues (penicillin-binding proteins and β-lactamases),^51,70,71 Lys65 in the catalytic triad was set as neutral, as it is considered to act as the general base for extracting the proton of Ser62. After equilibration, the amino group of the D-Ala-OEt substrate interacts with Asn155 and the oxygen group in the main chain of Ala480, while in the simulation with L-Ala-OEt, this residue interacts with the oxygen of the main chain of Ala288. Additionally, the methyl group of D-Ala-OEt is located in a hydrophobic pocket close to the residues Ala288, Ala480, and Trp220 (Fig. S3†).

The acylation reaction was modeled using three CVs (Fig. 1a) to drive the reaction from the initial reactant state to the acyl intermediate. CV₁ includes the number of bonds from the Ser62 oxygen to the hydrogens of Ser62, Lys65 and Tyr153, which accelerate proton transfer from Ser62 (Fig. 1a in green). CV₂ takes into account the transfer of a proton from Tyr153 to the ester group of Ala-OEt to form ethanol (Fig. 1a in orange). CV₃ samples the nucleophilic attack of the Ser62 O on C1 of the substrate and the cleavage of the ester bond of the substrate (Fig. 1a in blue). Attempts to model the acylation reaction with similar CVs used previously to model the same reaction³⁰ were not successful.

The 3D FELs obtained from the QM/MM simulations comparing the acylation reactions with L- and D-Ala-OEt are shown in Fig. 3a and b, respectively. The shapes of the FELs of the reaction with both substrates indicate a concerted mechanism with two energy minima. One is situated at the bottom of the FEL corresponding to the reactants (R), and the other is located at the top of the acyl intermediate (AI), and they are separated by a transition state (TS) (Fig. 3). The acylation reaction with D-Ala-OEt presents an energy barrier of 26.5 kcal mol⁻¹ (Table 1) and a R minimum with an energy of −30.2 kcal mol⁻¹. On the other hand, the DAP-mediated acylation with L-Ala-OEt presents a free energy barrier of 37.2 kcal mol⁻¹, which is 10.7 kcal mol⁻¹ higher than that with the D-amino acid (Table 1). The collective variables used to model the acylation reactions L- and D-Ala-OEt changed from CV₁ ∼1 bonds, CV₂ of approximately 3 Å and CV₃ of ∼−1.5 Å in the R, to be placed in the AI at approximately 0.2 bonds for CV₁, ∼1 Å for CV₂ and approximately 1.5 Å for CV₃. These changes take place because in the AI, the Ser62 oxygen has transferred its proton and is bound to C1, the proton of Tyr153 is covalently bound to the O1 atom making ethanol (CV₂), the distance between the oxygen of Ser62 to C1 is shorter to allow formation of the acyl intermediate (CV₃a), and the C1–O1 distance of Ala-OEt was increased due to the cleavage of the ester bond (CV₃b).

Fig. 3 FEL of the DAP-catalyzed acylation reaction visualized as three-dimensional surfaces for (a) L-Ala-OEt and (b) D-Ala-OEt. CV₁ (number of bonds), CV₂ and CV₃ are represented in x, y and z axes, respectively. R, reactants; TS, transition state; AI, acyl intermediate.

Table 1 Calculated free energy barriers (kcal mol⁻¹) for the acylation and aminolysis reactions, and the energy difference between the barriers for the L- and D-Ala-OEt isomers

DAP	L-Ala-OEt	D-Ala-OEt	ΔG L–D
Acylation	37.2	26.5	10.7
Aminolysis	44.5	28.8	17.7

---

The reactions with both substrates follow the same mechanism; they start with the transfer of the proton of Ser62 to neutral Lys65, protonating this residue (Fig. 4a–c, S4 green line and S5a–c†). This proton extraction is followed by the attack of the Ser62 oxygen on C1 atom of Ala-OEt and the shortening of the distance between the Tyr153 proton and the oxygen of the ester of Ala-OEt, bringing them quite close (Fig. 4d, e, S4b and d blue dotted line and S5d and e†). Afterward, the ester bonds of the tetrameric species elongate to the point of cleavage, and the proton of Tyr153 is transferred to the leaving group, releasing the newly formed ethanol molecule (Fig. 4e, f, S5e and f†). Finally, now-protonated Lys65 comes closer to deprotonated Tyr153, and a proton is transferred to negatively charged Tyr153, neutralizing it (Fig. 4g, h and S5†). The acylation mechanism observed in this work is analogous to the mechanism proposed by Khaliullin and collaborators.³⁵ One difference observed in the CVs as a function of time (Fig. S4†) is that, in the reaction with D-Ala-OEt, CV₃b monotonically decreases from 3 Å to ∼1.4 Å. That is, the distance between the C1⋯O1 atoms decreases until a covalent bond is formed (Fig. S4d†), and the distances in CV₂ are short. CV₁ (the number of bonds to the Ser62 oxygen) only changes after CV₃b begins to decrease. With L-Ala-OEt, the CV₃b distance is already approximately 2 Å and only decreases after the sharp drop of CV₁ (Fig. S4b†) from ∼1 to 0 bonds, followed also by a sharp drop in CV₂. This suggests that with D-Ala-OEt, the reaction occurs more cooperatively than with L-Ala-OEt, in which the driving event seems to be the attack of Ser62 to the C1 atom of the substrate. Previous QM studies with class-C β-lactamase enzymes, which have similar active sites, showed that the barriers to the acylation of its inhibitors were 22 kcal mol⁻¹ and 17 kcal mol⁻¹ for aztreonam and cephalothin,⁷² respectively, and 25 kcal mol⁻¹ for avibactam.⁷¹ The free energy barrier of the acylation of DAP with D-Ala-OEt is similar to those reported for inhibitors. However, L-Ala-OEt presents a higher energy barrier for acylation than D-Ala-OEt and the inhibitors modeled in other studies. The difference in the energy barrier between the L- and D-Ala-OEt substrates could be explained by the interaction between the amino group of D-Ala-OEt and Asn155. On the other hand, L-Ala-OEt does not form an interaction with Asn155, and its amino group is located near the main-chain oxygen of Ala288.

Fig. 4 Mechanism of the acylation reaction of DAP and L-Ala-OEt obtained from QM/MM ABMD simulations. Time-course results are shown from (a)–(h). DAP is represented in grey cartoon, atoms treated QM (Ser62, Lys65, Tyr153 and L-Ala-OEt) are shown in ball and stick representation. The carbon atoms of the L-Ala-OEt are shown in purple color. AI, acyl-intermediate; TS, transition state.

QM/MM aminolysis reaction

After modeling the acylation reaction, we simulated the aminolysis reaction starting from the acyl intermediate enzyme. We treated the same atoms used to model the acylation, that is, the catalytic triad, at the QM level of theory, and added the formed acyl intermediate and a new Ala-OEt monomer acting as a nucleophile. However, all attempts to model the reaction were unsuccessful or presented unreasonably high energetic barriers (>60 kcal mol⁻¹). During the equilibration of the system with classical MD, we observed that the loop on which Asp479 is located moved closer to the acyl intermediate and that Asp479 stabilized the acyl intermediate through a H-bond to its amino group. Asn155 has also been reported to be conserved in the active sites of proteins in the same family as DAP.³⁴ Additionally, we observed that a water molecule was located in the catalytic site close to Asp479, Asn155 and the amino groups of the attacking Ala-OEt nucleophiles and the acyl intermediate (Fig. S6†). Thus, we incorporated the Asn155 and Asp479 residues and the water near these residues into the QM system to model the aminolysis reaction.

The only CVs we found that were able to carry out the aminolysis reaction are shown in Fig. 1b. CV₁ increased the sampling of the nucleophilic attack of Ala-OEt to the C1 atom of the acyl intermediate, and the separation of the acyl intermediate from Ser62 is given by an LCOD (Fig. 1b in blue). CV₂, as the number of bonds to the N atom of Ala-OEt, is correlated with the acceleration of the transfer of one of the protons of the amino group of Ala-OEt to a nearby acceptor atom (Fig. 1b in orange). CV₃, as in CV₁ in the acylation reaction, represents the number of bonds between the Ser62 oxygen and the nearby H atoms, and these bonds accelerate the transfer of a proton to Ser62 to regenerate the enzyme (Fig. 1b in green).

The free energy landscapes of the aminolysis reactions with L- and D-Ala-OEt as substrates are presented in Fig. 5. The aminolysis with L-Ala-OEt presents an energy barrier of 44.5 kcal mol⁻¹, which is higher than the energy barrier with D-Ala-OEt (28.8 kcal mol⁻¹, Table 1). Thus, there is a greater difference between the energy barriers of the aminolysis reactions with the different substrates, 17.7 kcal mol⁻¹, than there was between the acylation reactions, for which the difference between the barriers with the Ala-OEt isomers was 10.7 kcal mol⁻¹. Thus, we concluded that aminolysis is the limiting step of the polymerization reaction for both substrates.

Fig. 5 FEL of the DAP-catalyzed aminolysis reaction visualized as three-dimensional surfaces for (a) L-Ala-OEt and (b) D-Ala-OEt. CV₁ is a LCOD and CV₂ and CV₃ are number of bonds, represented on the x, y and z axes, respectively. AI, acyl intermediate; TS, transition state; P, product. For clarity, only points with an energy lower than −10 kcal mol⁻¹ are shown.

The aminolysis reactions for D-Ala-OEt and L-Ala-OEt follow the same process. First, the nucleophilic attack by the N atom of the D- or L-Ala-OEt to the C1 atom of the acyl intermediate occurs, reducing its distance from approximately 2 Å to ∼1.6 Å, resulting in the formation of a covalent bond (Fig. 6a, b, S7, S8a and b†). However, to form the TS and advance the reaction, one of the H atoms of the amino group of the Ala-OEt nucleophile must be near the oxygen of Tyr153. At the same time, the proton of Tyr153 must be close to Lys65 (Fig. 6c and S8b†). Notably, Asn155 interacts with Lys65, helping it accept an additional proton from Tyr153. Then, proton exchange occurs; Lys65 extract a proton from Tyr153, and immediately following this, Tyr153 receives a proton from the amino group of the Ala-OEt nucleophile (Fig. 6c, d, S8c and d†). The deprotonated nucleophile then binds to the acyl intermediate, forming the tetrahedral structure seen in the TS (Fig. 6e, f and S8e†). If the oxygen of Tyr153 is not close enough to the nucleophile to extract a proton, the TS would be unstable, the reaction would not take place and the N of the Ala-OEt will separate from the C1 atom of the acyl intermediate. After the formation of this tetrahedral species, protonated Lys65 moves closer to Ser62, and the covalent bond between the Ser62 oxygen and the C1 atom of the acyl intermediate stretches until it breaks (Fig. 6f, g, S8f and g†). Then, Ser62 accepted a proton from protonated Lys65 (Fig. 6g, h and S8†).

Fig. 6 Mechanism of the aminolysis reaction of DAP with Ser62-acylated and D-Ala-OEt obtained from QM/MM ABMD simulations. Time-course results are shown from (a)–(h). DAP is represented in grey carton, atoms treated QM (Ser62, Lys65, Tyr153, Asn155, H₂O and D-Ala-OEt) are shown in ball and stick representation (Asp481, also treated as QM, is not shown for clarity). The carbon atoms of the acyl-intermediate are shown in green color, carbon atoms of D-Ala-OEt are displayed in pale green. AI, acyl-intermediate; TS, transition state.

We speculate that Asn155 and the water molecule, which are close during the course of the reaction, help reduce the energy barrier of the proton transfer from the amino group of the nucleophile to Tyr153, stabilizing the tetrahedral intermediate dipeptide species that retains the two protons on its amino group. In addition, Asn155 seems to stabilize the charge on Lys65, allowing it to accept a proton from Tyr153. However, in the reactions with L- and D-Ala-OEt, the distance between the water molecule and the implicated residues (Asn155 and the amino group of the nucleophile) are not significantly different.

There are some differences between the two systems, one being the orientation of the side-chain amino group in the acyl intermediate (in both systems, the methyl group of the acyl intermediate is located in the hydrophobic pocket formed by the residues Ala226 and Trp220). In the reaction of D-Ala-OEt, the amino group of the acyl intermediate is located above the methyl group, allowing it to interact with Asp479 and the side chains of Ala430 and Ala482. On the other hand, the amino group of the acyl intermediate in the aminolysis reaction of L-Ala-OEt is at the same level as the methyl group. This allows the amino group to interact with Asp479 and form an additional interaction with Asn155, stabilizing this species. Another difference between these two reactions is the distance between His285 and Tyr153 in the active site. The methyl group of the L-Ala-OEt nucleophile is oriented towards His285, and due to steric hindrance, it does not allow His285 to be close to Tyr153 (Fig. S9†). On the other hand, in the system with D-Ala-OEt, the methyl group of the Ala-OEt nucleophile is not oriented in that direction, and His285 can be closer to Tyr153, forming a stabilizing interaction with Tyr153 (Fig. S9†). Another difference that can be observed between the reactions is the dihedral angle, psi (ψ), between the N atom of the nucleophile and the C1 atom in the formed dipeptide. The newly formed peptides displayed different phi (φ) and psi (ψ) dihedral angles due to the position of their methyl groups, which changed the angle of the nucleophile attack. D-Ala-OEt presents a φ dihedral angle of ∼−160 degrees, whereas L-Ala-OEt presents an angle of approximately −110 degrees (Fig. S10†). The ψ angle displayed by the L-alanine dipeptide was ∼−20 degrees, and the dihedral angle of the D-alanine dipeptide was ∼158 degrees. In the case of L-Ala-OEt, these dihedrals angles are within the allowed alpha region. In contrast, the D-alanine dipeptide dihedrals are not in the allowed regions for the majority of the reaction, as the methyl of the second alanine and the amino N-terminal group are trans to each other (Fig. 6h).⁷³ This could partially explain the higher energy barrier for the conversion of the AI of the D-alanine dipeptide to the product (Fig. 5), but it does not account for the higher energy barrier of the aminolysis reaction involving L-Ala-OEt.

The energy obtained for the aminolysis reaction with DAP for D-Ala-OEt (28.8 kcal mol⁻¹) is similar to that previously calculated for papain with the same substrate.³⁰ However, in the case of L-Ala-OEt, the energy calculated here is significantly higher (44.5 kcal mol⁻¹) than that found with papain (12 kcal mol⁻¹).³⁰ In fact, we note that the barrier of the papain-mediated aminolysis reaction for D-Ala-OEt was previously calculated to be 28 kcal mol⁻¹, but papain cannot polymerize D-Ala-OEt. Here, we found that DAP is able to polymerize both substrates, but it is more efficient for D-Ala-OEt, which is consistent with the QM/MM results, but the energies presented here are higher than those presented for the papain-mediated polymerization.³⁰ Recently, the hydrolysis of the desulfated inhibitor avivactam by a class-C β-lactamase was reported with an energy barrier of 40 kcal mol⁻¹ and a t_1/2 of 10²–10⁴ h.⁷¹ Based on these previous reports, we consider that the high-energy barriers we encountered in the aminolysis reaction, especially for L-Ala-OEt, could be overestimations. This overestimation could be due to several factors; the main reason is our inability to reliably replicate the reaction, as all the reactions modeled using different CVs failed or displayed higher energetic barriers. The other reasons are diverse but could include the selection of the QM atoms, which is critical for accurately modeling the reactions. The active site of DAP is buried inside the protein and forms an extensive network of contacts and H-bonds that are difficult to model at the QM level. As stated before, incorporating additional atoms into our QM initial systems reduced the energy barrier, as shown by our first attempts to model the aminolysis reaction with DAP. In general, as the QM region size increases, more accurate results can be obtained,⁷⁴ but at the expense of increasing the computational calculation time. Thus, it is plausible that including more atoms in the QM region of our system would afford a more accurate model of the reaction. In addition, the selection of CVs to appropriately describe the reaction is known to be a critical step and a major limitation of metadynamics and similar techniques.⁶⁷ Different reactions have been modeled using a variety of CVs with no positive results. However, despite our best effort, the selected CVs used to finally model the aminolysis reaction may not be accurate enough to fully describe this complex reaction.

Previous works by Asano et al. showed that DAP displays higher stereospecificity towards D-amino acid derivatives, amides and esters as well as peptides with a D-alanine at their N terminus than towards the corresponding L-amino acid derivatives.^31,32 However, additional studies have shown that DAP can effectively polymerize D-alanine oligopeptides, although those reactions were performed in nonaqueous media.³⁸ Using immobilized DAP in water-saturated toluene and in the presence of an organic tertiary amine, DAP was able to catalyze the polymerization of D-alanine oligomers to give the dimer and trimer in 56% and 6% yield, respectively, from D-alanine methyl ester as the monomer. Here, we used DAP in an aqueous buffer without an organic solvent or tertiary amine, showing its ability to polymerize amino acid monomers into short polypeptides under mild conditions. Despite the fact that DAP is able to polymerize L- and D-Ala-OEt, the experimental analysis confirmed that it clearly has a preference for dextrorotatory substrates. Although DAP can catalyze the formation of L-alanine dipeptide, its activity is ∼26 times lower than its activity with D-Ala-OEt as the substrate. Indeed, the free energy barriers obtained from the QM/MM simulations for the reactions are in accordance with the experimental evidence. The barriers for the acylation (37.2 kcal mol⁻¹) and aminolysis (44.5 kcal mol⁻¹) with L-Ala-OEt are, in both cases, higher than the corresponding barriers for D-Ala-OEt (Table 1). D-Ala-OEt seems to be able to form more interactions with residues in the active site of the enzyme, such as with the conserved residues Asn155 and Asp479, and these interactions do not form in the acylation with L-Ala-OEt, and the methyl group of D-Ala-OEt is located closer to a hydrophobic pocket. The water molecules present in the active site are also likely involved in the aminolysis reaction, stabilizing the acyl intermediate and facilitating proton transfer between residues. In the aminolysis, the acyl intermediate seems more stable in the reaction with L-Ala-OEt, as the amino group interacts with Asp479 and Asn155; this last interaction is absent in the case of the AI with the dextrorotatory isomer. Thus, DAP preferring D-substrates can be attributed to chirality.

Conclusions

In this study, we demonstrated the DAP-mediated polymerization of D-Ala-OEt in an aqueous solution. The present results showed that DAP has higher activity with D-Ala-OEt than with L-Ala-OEt. The QM/MM ABMD results were in agreement with the experimental results. The calculated energy landscapes indicate that both the D- and L-substrates underwent concerted reactions, with aminolysis being the limiting step in the catalytic step because it has higher activation energies for both substrates. The reaction of L-Ala-OEt with DAP displayed higher energy barriers for both acylation and aminolysis relative to those of the D-Ala-OEt substrate. In the acylation reaction, the interaction between the amino group of the substrate with Asn155 reduced the energetic barrier faced by D-Ala-OEt, whereas in the aminolysis reaction, the orientation of the amino group in L-Ala-OEt towards Asp479 and Asn155 stabilized this species, and the orientation of the attacking nucleophile could increase the energetic barrier relative to that with D-Ala-OEt. The combined experimental and computational results could provide a platform for synthesizing polypeptides that incorporate D-amino acids and help clarify the stereospecificity of proteases. The possibility of using DAP in aqueous buffers without organic solvents broadens the applications chemoenzymatic synthesis and allows the incorporate D-amino acids as substrates. With their polymerization and incorporation into polypeptides, more biostable peptides could be achieved, and such species may offer benefits for biotechnological and pharmaceutical industries for use as inhibitors or delivery systems.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This study was financially supported by JST ERATO (grant number JPMJER1602). We acknowledge RIKEN ACCC for providing access to the Hokusai/GreatWave and BigWaterfall supercomputer resources.

References

1. A. D. Radkov and L. A. Moe, Appl. Microbiol. Biotechnol., 2014, 98, 5363.
2. I. Güell, J. Cabrefiga, E. Badosa, R. Ferre, M. Talleda, E. Bardají, M. Planas, L. Feliu and E. Montesinos, Appl. Environ. Microbiol., 2011, 77, 2667.
3. S. Fernandez-Lopez, H.-S. Kim, E. C. Choi, M. Delgado, J. R. Granja, A. Khasanov, K. Kraehenbuehl, G. Long, D. A. Weinberger, K. M. Wilcoxen and M. R. Ghadiri, Nature, 2001, 412, 452.
4. G. Genchi, Amino Acids, 2017, 49, 1521.
5. M. Kuwada, T. Teramoto, K. Y. Kumagaye, K. Nakajima, T. Watanabe, T. Kawai, Y. Kawakami, T. Niidome, K. Sawada and Y. Nishizawa, Mol. Pharmacol., 1994, 46, 587.
6. A. M. Torres, C. Tsampazi, D. P. Geraghty, P. S. Bansal, P. F. Alewood and P. W. Kuchel, Biochem. J., 2005, 391, 215.
7. G. Kreil, Annu. Rev. Biochem., 1997, 66, 337.
8. V. Erspamer, P. Melchiorri, G. Falconieri-Erspamer, L. Negri, R. Corsi, C. Severini, D. Barra, M. Simmaco and G. Kreil, Proc. Natl. Acad. Sci. U. S. A., 1989, 86, 5188.
9. A. Vater and S. Klussmann, Drug Discovery Today, 2015, 20, 147.
10. P. Vlieghe, V. Lisowski, J. Martinez and M. Khrestchatisky, Drug Discovery Today, 2010, 15, 40.
11. Z. Feng and B. Xu, Biomol. Concepts, 2016, 7, 179.
12. K. Fosgerau and T. Hoffmann, Drug Discovery Today, 2015, 20, 122.
13. L. Zhao and W. Lu, Curr. Opin. Chem. Biol., 2014, 22, 56.
14. A. Rodriguez-Granillo, S. Annavarapu, L. Zhang, R. L. Koder and V. Nanda, J. Am. Chem. Soc., 2011, 133, 18750.
15. H. E. Stanger and S. H. Gellman, J. Am. Chem. Soc., 1998, 120, 4236.
16. K. M. Makwana and R. Mahalakshmi, Pept. Sci., 2016, 106, 260.
17. L. E. Zawadzke and J. M. Berg, Proteins: Struct., Funct., Bioinf., 1993, 16, 301.
18. L.-W. Hung, M. Kohmura, Y. Ariyoshi and S.-H. Kim, Acta Crystallogr., Sect. D: Biol. Crystallogr., 1998, 54, 494.
19. M. Novotny and G. J. Kleywegt, J. Mol. Biol., 2005, 347, 231.
20. A. Aliashkevich, L. Alvarez and F. Cava, Front. Microbiol., 2018, 9, 683.
21. A. E. Rabideau and B. L. Pentelute, ACS Cent. Sci., 2015, 1, 423.
22. R. Milton, S. Milton and S. Kent, Science, 1992, 256, 1445.
23. D. F. Veber, F. W. Holly, R. F. Nutt, S. J. Bergstrand, S. F. Brady, R. Hisrschmann, M. S. Glitzer and R. Saperstein, Nature, 1979, 280, 512.
24. M. Uppalapati, D. J. Lee, K. Mandal, H. Li, L. P. Miranda, J. Lowitz, J. Kenney, J. J. Adams, D. Ault-Riché, S. B. H. Kent and S. S. Sidhu, ACS Chem. Biol., 2016, 11, 1058.
25. J. R. Pappenheimer, C. E. Dahl, M. L. Karnovsky and J. E. Maggio, Proc. Natl. Acad. Sci. U. S. A., 1994, 91, 1942.
26. K. Yazawa and K. Numata, Molecules, 2014, 19, 13755.
27. K. G. Varnava and V. Sarojini, Chem.–Asian J., 2019, 14, 1088.
28. F. Bordusa, Chem. Rev., 2002, 102, 4817.
29. P. J. Baker and K. Numata, Biomacromolecules, 2012, 13, 947.
30. J. Gimenez-Dejoz, K. Tsuchiya and K. Numata, ACS Chem. Biol., 2019, 14, 1280.
31. Y. Asano, A. Nakazawa, Y. Kato and K. Kondo, J. Biol. Chem., 1989, 264, 14233.
32. Y. Asano, Y. Kato, A. Yamada and K. Kondo, Biochemistry, 1992, 31, 2316.
33. C. Bompard-Gilles, H. Remaut, V. Villeret, T. Prangé, L. Fanuel, M. Delmarcelle, B. Joris, J. M. Frère and J. Van Beeumen, Structure, 2000, 8, 971.
34. J. Arima, K. Shimone, K. Miyatani, Y. Tsunehara, Y. Isoda, T. Hino and S. Nagano, FEBS J., 2016, 283, 337.
35. I. G. Khaliullin, D. a. Suplatov, D. N. Shalaeva, M. Otsuka, Y. Asano and V. K. Svedas, Acta Naturae, 2010, 2, 66.
36. J. Arima, H. Usuki, T. Hatanaka and N. Mori, Appl. Environ. Microbiol., 2011, 77, 8209.
37. W. Lee, M. A. McDonough, L. P. Kotra, Z.-H. Li, N. R. Silvaggi, Y. Takeda, J. A. Kelly and S. Mobashery, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 1427.
38. Y. Kato, Y. Asano, A. Nakazawa and K. Kondo, Biocatalysis, 1990, 3, 207.
39. B. D. Welch, A. P. VanDemark, A. Heroux, C. P. Hill and M. S. Kay, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 16828.
40. J. A. Chuah, T. Yoshizumi, Y. Kodama and K. Numata, Sci. Rep., 2015, 5, 7751.
41. J.-A. Chuah and K. Numata, Biomacromolecules, 2018, 19, 1154.
42. F. Milletti, Drug Discovery Today, 2012, 17, 850.
43. R. E. W. Hancock and D. S. Chapple, Antimicrob. Agents Chemother., 1999, 43, 1317.
44. F. Cava, H. Lam, M. A. de Pedro and M. K. Waldor, Cell. Mol. Life Sci., 2011, 68, 817.
45. T. Kigawa, T. Yabuki, N. Matsuda, T. Matsuda, R. Nakajima, A. Tanaka and S. Yokoyama, J. Struct. Funct. Genomics, 2004, 5, 63.
46. K. Numata, Y. Motoda, S. Watanabe, N. Tochio, T. Kigawa and Y. Doi, Biomacromolecules, 2012, 13, 3450.
47. A. Spirin, V. Baranov, L. Ryabova, S. Ovodov and Y. Alakhov, Science, 1988, 242, 1162.
48. J. M. Ageitos, K. Yazawa, A. Tateishi, K. Tsuchiya and K. Numata, Biomacromolecules, 2016, 17, 314.
49. O. Trott and A. Olson, J. Comput. Chem., 2010, 31, 455.
50. G. Morris and R. Huey, J. Comput. Chem., 2009, 30, 2785.
51. J. C. Hargis, J. K. White, Y. Chen and H. L. Woodcock, J. Chem. Inf. Model., 2014, 54, 1412.
52. D. A. Case, D. S. Cerutti, T. E. I. Cheatham, T. A. Darden, R. E. Duke, T. J. Giese, H. Gohlke, A. W. Goetz, D. Greene, N. Homeyer, S. Izadi, A. Kovalenko, T. S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, D. Mermelstein, K. M. Merz, G. Monard, H. Nguyen, I. Omelyan, A. Onufriev, F. Pan, R. Qi, D. R. Roe, A. Roitberg, C. Sagui, C. L. Simmerling, W. M. Botello-Smith, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu, L. Xiao, D. M. York and P. A. Kollman, AMBER 2016, University of California, San Francisco, 2016.
53. J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, C. Simmerling, S. Biology, S. Brook, S. Brook, Q. Biology and S. Brook, J. Chem. Theory Comput., 2015, 11, 3696.
54. J. Wang, R. M. Wolf, J. W. Caldwell and P. A. Kollman, J. Comput. Chem., 2004, 25, 1157.
55. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926.
56. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327.
57. T. E. I. Cheatham, J. L. Miller, T. Fox, T. A. Darden and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 4193.
58. D. R. Roe and T. E. Cheatham III, J. Chem. Theory Comput., 2013, 9, 3084.
59. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33.
60. G. A. Petersson, A. Bennett, T. G. Tensfeldt, M. A. Allaham, W. A. Shirley and J. Mantzaris, J. Chem. Phys., 1988, 89, 2193.
61. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 09, Revision A.02, 2016.
62. A. W. Götz, M. A. Clark and R. C. Walker, J. Comput. Chem., 2014, 35, 95.
63. R. C. Walker, M. F. Crowley and D. A. Case, J. Comput. Chem., 2008, 29, 1019.
64. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1986, 7, 718.
65. V. Babin, C. Roland and C. Sagui, J. Chem. Phys., 2008, 128, 134101.
66. V. Babin and C. Sagui, J. Chem. Phys., 2010, 132, 104108.
67. A. Barducci, M. Bonomi and M. Parrinello, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 826.
68. A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562.
69. B. Ensing, A. Laio, M. Parrinello and M. L. Klein, J. Phys. Chem. B, 2005, 109, 6676.
70. B. F. Gherman, S. D. Goldberg, V. W. Cornish and R. A. Friesner, J. Am. Chem. Soc., 2004, 126, 7652.
71. C. K. Das and N. N. Nair, Phys. Chem. Chem. Phys., 2018, 20, 14482.
72. R. Tripathi and N. N. Nair, J. Am. Chem. Soc., 2013, 135, 14679.
73. C.-L. Towse, G. Hopping, I. Vulovic and V. Daggett, Protein Eng., Des. Sel., 2014, 27, 447.
74. H. J. Kulik, J. Zhang, J. P. Klinman and T. J. Martínez, J. Phys. Chem. B, 2016, 120, 11381.

---

Footnote

† Electronic supplementary information (ESI) available: Fig. S1 to S9. See DOI: 10.1039/d0ra01138j

---