The Role of Hydrogen Bond in Catalytic Triad of Serine Proteases†

In order to investigate the origin of catalytic power for serine proteases, the role of the hydrogen bond in the catalytic triad was studied in the proteolysis process of the peptides chymotrypsin inhibitor 2 (CI2), MCTI-A, and a hexapeptide (SUB), respectively. We first calculated the free energy profile of the proton transfer between His and Asp residues of the catalytic triad in the enzyme-substrate state and transition state by employing QM/MMmolecular dynamics simulations. The results show that a low-barrier hydrogen bond (LBHB) only forms in the transition state of the acylation of CI2, while it is a normal hydrogen bond in the acylation of MCTI-A or SUB. In addition, the change of the hydrogen bond strength is much larger in CI2 and SUB systems than in MCTI-A system, which decreases the acylation energy barrier significantly for CI2 and SUB. Clearly, a LBHB formed in the transition state region helps accelerate the acylation reaction. But to our surprise, a normal hydrogen bond can also help to decrease the energy barrier. The key to reducing the reaction barrier is the increment of hydrogen bond strength in the transition state state, whether it is a LBHB or not. Our studies cast new light on the role of the hydrogen bond in the catalytic triad, and help to understand the catalytic triad of serine proteases.


I. INTRODUCTION
Serine proteases are one of the largest groups of enzymes, which are ubiquitous throughout the living system, and responsible for many different biologic functions [1]. The essential functional unit of conventional serine proteases family is usually a catalytic triad (Ser/His/Asp) [1,2], in which serine serves as a nucle- † Part of Special Issue "John Z.H. Zhang Festschrift for celebrating his 60th birthday".
* Authors to whom correspondence should be addressed. E-mail: zhouyz@nju.edu.cn, dqxie@nju.edu.cn ophilic attacking group to initiate the whole hydrolysis process, histidine acts as a general acid-base element, and aspartic acid is assumed to be responsible for stabilizing the positive charge on histidine imidazole ring in transition states mainly through electrostatic interaction. However, the origin of catalytic power for serine proteases is still unclear and has long been an important issue. The catalytic role of the low-barrier hydrogen bond (LBHB) in enzymatic reactions was proposed nearly thirty years ago [3][4][5][6], which can supply extra stabilization to transition states or intermediates and enhance the reaction rate. In the enzyme-substrate complex, the hydrogen bond formed between His and Asp of the catalytic triad is a normal weak hydrogen bond, while in the transition state, a LBHB formed. In LBHB, the proton is delocalized between the two heteroatoms, which makes it stronger than a normal hydrogen bond. The increased strength of the LBHB can provide extra energy to overcome the reaction barrier [7]. But it is elusive and debatable about the existence and function of LBHB in enzyme catalysis [8,9].
Due to the high efficiency of enzymes, transition states and intermediates formed in the catalytic process are very short-lived, therefore the LBHB is inherently formidable to prove. In experiments, crystal structure, nuclear magnetic resonance spectrum, neutron diffraction, infrared spectrum, isotope separation experiments are usually used to examine the property and strength of hydrogen bonds in enzymes [10][11][12]. Hydrogen bonding interactions have also been extensively studied by theoretical methods [10,13,14]. However, the results are very controversial. Numerous studies have shown that a LBHB can be formed in the transition state and promote the forward reaction [4,6,[15][16][17]. While some studies suggested that the LBHB does not exist in the enzyme at all [13,14,[18][19][20]. In addition, other experimental evidences and theoretical calculations found that the LBHB does exist in NMR spectrum or X-ray structure, but it is much weaker in protein than in the gas phase and unlikely to have a significant catalytic effect [21][22][23]. A recent work pointed out the existence of LBHB requires certain conditions [24], which suggests that there is no universal answer for this topic.
Our group has been working on serine proteases for many years, which provides us a good opportunity to study on this special topic of LBHB in serine proteases. In our previous studies on the mechanism of peptides hydrolyzed by serine proteases, we only recognized a LBHB for subtilisin-CI2 system. In this work, we will perform further theoretical investigations for the properties and functions of His· · · Asp hydrogen bond in the acylation reaction of peptides CI2, MCTI-A and a hexapeptide SUB hydrolyzed by serine proteases employing Born-Oppenheimer ab initio QM/MM molecular dynamics simulations and umbrella sampling. The hydrogen bond in the enzymatic environment is of essential importance for understanding protein structure and predicting the activity of enzymes. Our studies would throw new insights on the role of hydrogen bonds in enzymatic reactions, and help understand the mechanism of the catalytic triad in serine proteases.

II. COMPUTATIONAL DETAILS
The initial structures of the enzyme-peptide complexes were obtained from the X-ray crystal structures of subtilisin BPN ′ in complex with chymotrypsin inhibitor 2 (CI2) (PDB ID: 1LW6) [25], and porcine βtrypsin with MCTI-A (PDB ID: 1MCT) [26], respectively. The sequence of the hexapeptide SUB was taken from the MCTI-A at the active site (Cys3-Pro4-Arg5-Ile6-Trp7-Met8). The protonation states of charged residues were determined at constant pH=7 based on pK a calculations via the PROPKA program [27] and the consideration of the local hydrogen bonding network. Next, each system was neutralized, solved, and equilibrated with a series of minimizations interspersed by short molecular dynamics simulations with periodic boundary condition. Finally, an extensive molecular dynamics simulation of 100 ns at constant temperature and pressure was carried out and a snapshot was randomly chosen for the subsequent QM/MM simulations. The pressure was maintained at 1 atm, and the temperature was controlled at 298 K for CI2 and 310 K for MCTI-A and SUB.
All ab initio QM/MM calculations were performed with modified Q-Chem [28] and Amber12 programs [29] for CI2 system. The QM sub-system includes the side chains of the catalytic triad of the target enzyme and the scissile peptide portion of the peptides, and was described at the B3LYP/6-31++G * * level. All other atoms were described by amber99SB force field [30,31]. The QM/MM boundary was treated by the pseudobond approach [32,33] with the improved parameters [34]. There are 47 QM atoms and 6 pseudo atoms for CI2 system, and 38 QM atoms and 7 pseudo atoms for MCTI-A system, respectively. For more details, please refer to our previous work [35][36][37].
The prepared QM/MM system was first minimized and then employed to map out a reaction path with QM/MM minimizations. The reaction coordinates were chosen as d N−H −d O−H for proton transfer and the distance between N and O atoms was chosen for the hydrogen bond strength. With the chosen reaction coordinate, first a minimum energy path was mapped out with ab initio QM/MM restrained minimizations using the reaction coordinate driving (RCD) method [38]. For each determined structure along the path, a 500 ps classical molecular dynamics (MD) simulation was performed to equilibrate the MM subsystem, with the QM subsystem being frozen. Since the length of ab ini-FIG. 1 (a) Illustration of acylation reaction mechanism for hydrolysis of peptides by serine proteases. The peptide is in blue, while the serine protease is in black. For MCTI-A or SUB system, the catalytic triad is "Ser195/His57/Asp102". While for CI2 system, it is "Ser221/His64/Asp32". (b) Illustration of the proton transfer between His and Asp in serine protease-peptide complex in transition state.
tio QM/MM MD simulations is quite short, this step is important that the MM sub-system can be relaxed given the change of QM sub-system. The final snapshot from the MM sub-system equilibration was used as the starting structure for Born-Oppenheimer QM/MM MD simulations with umbrella sampling for that window. The window size was chosen as 0.2Å along the reaction path, and 0.1Å around the transition state. Along the reaction path, about 10 windows were chosen for each step. For every window, at least 40 ps simulations have been carried out and a harmonic potential force constant of 30−50 kcal·mol −1 ·Å −2 was employed. The last thirty ps QM/MM MD simulations were used to calculate the free energy profile with the weighted histogram analysis method (WHAM) [39,40].

A. The mechanism for acylation reaction of serine proteases
Serine proteases hydrolyze peptides or proteins in two successive stages: acylation and deacylation. In this work, we focus on the acylation stage of CI2 hydrolyzed by subtilisin BPN ′ , and MCTI-A or SUB hydrolyzed by porcine beta-trypsin on the QM/MM level. CI2 [41] is a member of the potato inhibitor 1 family of serine proteinase inhibitors, which is a well-defined peptide without disulfide bonds. The scissile peptide bond of CI2 lies between Met59 and Glu60, and a catalytic triad (Ser221, His64, and Asp32) is employed as essential catalytic functional unit by the enzyme subtilisin. MCTI-A [26] is a potent trypsin peptide inhibitor of the squash family with 28 amino acid residues and a cysteine-knot structure. The sequence of the SUB was derived from MCTI-A at the active site (Cys3-Pro4-Arg5-Ile6-Trp7-Met8). For MCTI-A and SUB, the scissile bond lies between Arg5-Ile6 and the catalytic triad of porcine βtrypsin is Ser195, His57 or Asp102. For convenience, we will use CI2, MCTI-A or SUB system to represent the corresponding serine protease-peptide systems.
From our previous QM/MM-MD simulations, the acylation reactions of the three peptides share a similar two-step catalytic mechanism as shown in FIG. 1(a) [35][36][37]. First, in the enzyme-substrate state (ES), the OG atom of Ser nucleophilic attacks the carbonyl carbon atom of the scissile peptide bond to form a C−O bond. Meanwhile the proton transfers to the N atom of His, and then the first tetrahedral intermediate TI1 is formed. In the next step, His serves as a general acid to protonate the leaving group for the breaking of the TI1, the C−N bond of peptide breaks leading to a stable product EA1. During the reaction, aspartic acid serves to orient His57 and to neutralize the charge developed on it in transition states responsible for stabilizing the transition states. The free energy profiles in FIG. 2 show that the reaction barrier of MCTI-A is significantly high in comparison with CI2 and SUB. One interesting finding is that for CI2 system, spontaneous transfer of the proton between the N atom of imidazole ring of His64 and the O atom of the carboxyl group of Asp32 occurs from time to time in the transition state region as shown in FIG. 1(b), which suggests that the hydrogen bond between His64 and Asp32 might be a low-barrier hydrogen bond (LBHB). While for MCTI-A and SUB systems, however, this proton transfer phenomenon never happens. This indicates that the properties of the hydrogen bond between His and Asp of the serine proteases are different in different serine proteases-peptide complexes. Under which conditions, the LBHB would form in the catalytic triad? And will this LBHB help to decrease the energy barrier? What's the function of a hydrogen bond in the catalytic triad? To answer those questions, we would investigate this hydrogen bond in the catalytic triad thoroughly.

B. The properties of the hydrogen bond in the catalytic triad
The His· · · Asp hydrogen bond is inherently formidable to prove in biomolecules experimentally due to the short-lived nature of the transition state. The calculated energy of hydrogen bond proton transfer is considered to be a straightforward way to determine the property of hydrogen bond theoretically  [14,23]. Therefore, we determined the free energy profile for the proton transfer from His to Asp in the enzyme-substrate state ES and transition state TS2 * with QM/MM MD simulations and umbrella sampling. The reaction coordinate is chosen as d N−H −d O−H , along which the H atom transfers from His to Asp. In FIG. 3(a), the free energy profile shows that in the TS2 * state of CI2, the free energy barrier of the proton transfer from His to Asp is only about 0.3 kcal/mol, which indicates a LBHB formed, in consistence with the proton transferring phenomenon observed in the transition state region of CI2. On the other hand, it can be clearly seen from FIG. 3(b) that the free energy curves of proton transfer only have a single well for the three systems in the ES state and for MCTI-A or SUB system in the TS2 * state, which indicates that the hydrogen bonding proton is located only on the side of His57 as a normal hydrogen bond. In a word, the His· · · Asp hydrogen bond is a normal hydrogen bond in the reactant state for the three systems we   Table I. Clearly, this distance is significantly shorter in the LBHB (2.565Å) than in other normal hydrogen bonds (longer than 2.7Å), which indicates that distance is a key factor to form LBHB. From ES to TS2 * , the distance becomes shorter, because the interaction between His and Asp is stronger in TS2 * due to the accumulated positive charge on His. The distance at ES for CI2 is much shorter than that of MCTI-A or SUB, so it is easier for CI2 to form a LBHB in TS2 * state.
Because the distance between the heteroatoms is always correlated with the strength of hydrogen bond [42,43], the change in the distance gives us a hint about how the strength of hydrogen bond changes along the reaction path. For SUB or CI2, the change of distance is more than 0.1Å, while for MCTI-A, the distance is nearly unchanged possibly due to the rigidity of MCTI-A caused by the cysteine-knot. Therefore, we expected that the hydrogen bond would be much stronger in TS2 * state for CI2 and SUB, but with little change for MCTI-A.

C. The strength and function of the hydrogen bond
The strength of hydrogen bond is important for understanding protein interactions and estimating the activation energy of enzyme reactions, so it would be of great significance to have a direct method to estimate the hydrogen bond strength. The strength of HB can be calculated as the energy difference between the total energy of the whole system and the sum of the individual energies of the donor and acceptor fragments [44][45][46]. But it is impossible to divide the enzyme into two parts under its complicated bonding and non-bonding interaction pattern. Hydrogen bond energy acts as a function of separation between donor and acceptor, and reaches equilibration when the distance is enough long [44]. So we can calculate the hydrogen bond energy by pulling the two fragments far away and the energy required to pull is the bond strength. But in a realistic enzyme system, it would be impossible to calculate the absolute strength of a hydrogen bond by pulling the donor and acceptor too far apart [47,48]. Typical hydrogen bonds in proteins are 2.8−3.0Å long between N and O atoms, extending from the natural bond length of the hydrogen bond to 3.5Å which would be enough to break a hydrogen bond within enzyme, so the energy required for this process can be approximately used as the strength of the hydrogen bond [47,48]. According to this we calculated the strength of His· · · Asp hydrogen bond in the transition state TS2 * and reactant state ES of the acylation process with QM/MM-MD simulations and umbrella sampling by elongating the N−O distance between His and Asp to 3.5Å.
The structure and strength of hydrogen bonds are usually strongly correlated with the distance between the donor and acceptor atoms [42,43,49]. The relationship between hydrogen bond strength and the distance between heavy atoms has been widely discussed in the gas and liquid phases. But, the quantification of the hydrogen bond strength and the distance between heavy atoms in enzymes has been rarely reported, which may be significantly changed by specific enzyme environments. Such enzymatic environment in a hydrogen bond strength is important to the understanding of protein interactions, including drug design and predicting enzymatic reactions barrier. To fill the gap of understanding this problem, we elucidated the relationship between the heavy atom distance of the hydrogen bond and the corresponding hydrogen bond strength as shown in FIG. 4. The results showed that there was also a linear correlation between N−O distance and the hydrogen bond strength. This means that the hydrogen bond strength within the enzyme can also be predicted by the distance between the donor and acceptor atoms.
The contribution of the hydrogen bond to decreasing the acylation barrier is estimated by ∆G HB , which is the difference between the strength of the hydrogen bond in TS2 * and ES states: FIG. 5 illustrates the strength of His· · · Asp hydrogen bond in the ES and TS2 * states for the three systems. These energy profiles show that LBHB is stronger than the normal hydrogen bonds. As shown in FIG. 5 and Table I, the LBHB strength of His64· · · Asp32 in TS2 * state of the CI2 system is 10.4±0.1 kcal/mol, which is 4.6 kcal/mol higher than that in ES state. This means that LBHB can provide stabilization energy of 4.6 kcal/mol to decrease the energy barrier of acylation. These results provide a direct theoretical evidence for the proposal that LBHB promotes the enzyme catalytic efficiency. The HB provides 4.5 kcal/mol extra stabilization energy for the acylation reaction of peptide SUB, while only provides about 1.0 kcal/mol for MCTI-A. This finding is consistent with the acylation energy barrier of the three systems as listed in Table I. These results mainly suggest two critical conclusions: first, the strength of a normal hydrogen bond is weaker than a LBHB; secondly, the His· · · Asp hydrogen bond in TS2 * state is stronger than at ES, reflecting the stability role of the hydrogen bond in the transition state. By comparing the strength of hydrogen bonds in the three systems, we find that apart from LBHB, normal hydrogen bonds can also reduce the energy barrier of the reaction. How much the reaction energy barrier can be reduced by the hydrogen bond in the catalytic triad is determined by the increment of hydrogen bond strength in the transition state, rather than the property of the hydrogen bond itself, which provides a new insight into the role of the hydrogen bond in enzyme reaction.

IV. CONCLUSION
In summary, we have calculated the free energy profile of the proton transfer within the His· · · Asp hydrogen bond of the catalytic triad of serine protease and also estimated its contribution to decreasing the acylation energy barrier by employing Born-Oppenheimer ab initio QM/MM MD simulation and umbrella sampling. The results show that the His· · · Asp hydrogen bond is a normal HB in ES state, and becomes LBHB in TS2 * state only for CI2 system. One explanation of the LBHB hypothesis has been proved that it plays an important role in reducing the barrier of enzyme reaction. Different from previous understanding in this field, we find that the property of the hydrogen bond is not the original source to decrease the reaction barrier, but the increment of hydrogen bond strength in transition state is the key factor. Apart from LBHB, normal hydrogen bonds may also reduce the energy barrier of the acylation reaction. In addition, we constructed a quantitative relationship between the strength of the hydrogen bond and the distance of N−O atoms between His and Asp in the enzyme, intuitively showing a linear correlation between the distance and strength of the hydrogen bond. Our studies on the above systems throw new insights and direct theoretical basis for understanding the role of the hydrogen bond in enzyme reaction, and emphasize the importance of hydrogen bond in regulating the reaction rate and understanding the reaction mechanism.

V. ACKNOWLEDGMENTS
This work is supported by the National Key Research and Development Program of China (2017YFA0206500) and the National Natural Science Foundation of China (No.22073040). We thank the High Performance Computing Center of Nanjing University for providing computational resources.