Consecutive Marcus Electron and Proton Transfer in Heme Peroxidase Compound II-Catalysed Oxidation Revealed by Arrhenius Plots

Electron and proton transfer reactions in enzymes are enigmatic and have attracted a great deal of theoretical, experimental, and practical attention. The oxidoreductases provide model systems for testing theoretical predictions, applying experimental techniques to gain insight into catalytic mechanisms, and creating industrially important bio(electro)conversion processes. Most previous and ongoing research on enzymatic electron transfer has exploited a theoretically and practically sound but limited approach that uses a series of structurally similar (“homologous”) substrates, measures reaction rate constants and Gibbs free energies of reactions, and analyses trends predicted by electron transfer theory. This approach, proposed half a century ago, is based on a hitherto unproved hypothesis that pre-exponential factors of rate constants are similar for homologous substrates. Here, we propose a novel approach to investigating electron and proton transfer catalysed by oxidoreductases. We demonstrate the validity of this new approach for elucidating the kinetics of oxidation of “non-homologous” substrates catalysed by compound II of Coprinopsis cinerea and Armoracia rusticana peroxidases. This study – using the Marcus theory – demonstrates that reactions are not only limited by electron transfer, but a proton is transferred after the electron transfer event and thus both events control the reaction rate of peroxidase-catalysed oxidation of substrates.

Oxidation and reduction reactions are the most common biochemical transformations in biology. They are catalysed by oxidoreductases, a diverse class of enzymes evolved to facilitate electron transfer (ET). Therefore, an improved qualitative and quantitative understanding of ET catalysed by these enzymes is pivotal to advance practical applications, such as the selection of mediators for bioconversion systems or the rational engineering of enzymes for bioelectrochemical sensors, biofuel cells, and other technological processes. For example, peroxidases have applications in feedstock delignification 1 , dye decolorization 2 , xenobiotics degradation 3 , and biosensing 4 . As a result, ET is a common subject of scientific research, especially for enzymes of current practical importance, such as heme peroxidase (HP).
HP is a diverse class of oxidoreductases found in fungi, plants, and animals. These enzymes catalyse oxidation reactions via the following reaction scheme 5  Here, the reduced form of enzyme E red is oxidized by hydrogen peroxide. The resultant intermediate form compound I (CpdI) is reduced back to E red via CpdII with two reduction steps (each consisting of one electron and one proton transfer). Possibly the most remarkable feature of HP is its wide range of possible substrates, S, that can be oxidised. These include phenols, indoleacetic acids, phenylendiamines, phenoxazines,and phenothiazines [6][7][8][9][10] .
Notably, the most comprehensive data are available for the rate constant k 2 and a range of substrates, because it represents a much slower step than k 1 [6][7][8][9][10] . On one hand, this lack of specificity allows for a wide range of applications of HP enzymes, e.g. horseradish peroxidase (HRP) or with the lesser known, but more active peroxidase from Coprinopsis cinerea (CIP). On the other hand, it may appear puzzling how these enzymes are able to oxidise such a variety of substrates. There have been various attempts resolve this apparent paradox [6][7][8][9][10] , all using Marcus theory, created in the late 1960s and laid out in seminal papers relating the bimolecular ET reaction rate constant with properties of reactants via the famous equation [11][12][13][14][15][16][17][18] : where V if is the electronic coupling between the initial and final state electronic wave functions, T the temperature, k B the Boltzmann constant, ħ the reduced Planck constant, λ the solvent and inner reorganisation energy of a reaction, ΔG 0 Gibbs energy of the reaction, and K eq the equilibrium constant of transitional complex formation from reactants. Previous studies attempted to explain the observed oxidation rate constants,k 1 and k 2 from Eq. 1, based on the physicochemical properties of substrates 6-10 , measuring rate constants using different series of homologous compounds. This led to the axiom that homologous compounds should have similar reaction reorganisation energies and pre-exponential factors in Eq. 2. Therefore, the only variable that changes for different sets of homologous compounds is ΔG 0 and there should exist a parabolic relation between ln(k ET ) and ΔG 0 . However, this relation can only exist if ET limits the reaction rate. This type of reasoning has also been applied to other oxidoreductases [19][20][21][22][23][24][25][26][27] . In fact, semi-parabolic dependencies between both k 1 and k 2 and ΔG 0 were found for both CpdI and CpdII of HRP and CIP and a variety of substrates, which led to the conclusion that both k 1 and k 2 are ET [6][7][8][9][10] . This conclusion appears surprising because the substrate oxidation steps II and III in Eq. 1 involve both ET and proton transfer (PT). Due to the greater weight of a proton, the PT is generally much slower than an ET, suggesting that k 1 and k 2 are PT rather than ET limited, which calls the previous research on HP catalysis into question.
We believe that this contradiction arises from the use of homologous series of compounds, which limits the range of ΔG 0 and makes it impossible to investigate the entire range of catalytic properties of HPs and other oxidoreductases. Therefore, we suggest that the investigation of non-homologous series of compounds can yield an explanation that is in better agreement with the points raised above. To test this, we selected CpdII from HRP and CIP and a set of substrates that had been used in previous publications, which differ in reduction potential, structure, and the number of electrons and protons donated during the oxidation process (12 substrates, 13 oxidation reactions in total). To provide a picture of the CpdII reduction mechanism that is as complete as possible, we also measured the temperature dependencies of k 2 for each reaction, calculated the quantum chemical self-exchange solvent and inner reorganisation energies (λ s ), estimated the change of λ s for compounds in enzyme substrate complexes, and measured the kinetic isotope effects for the oxidation rate of the chosen compounds.

Results and Discussion
Kinetics and thermodynamics of substrates catalysed by cpdii. We measured the rate constants (k 2 ) and activation free energies (ΔG ‡ ) with pre-exponential factors (A) for all 12 substrates at pH 7.00 (experimental details, structures, and names of substrates are available in the Supporting Information (SI)). At 25 °C, the rate constants ranged from about 2 × 10 6 to 5 × 10 8 M −1 s −1 for CIP and from about 1 × 10 5 to 3 × 10 8 M −1 s −1 for HRP (concerns about possible diffusion-limited kinetics are addressed in the SI). These values are comparable with those published previously 9,10 . Measured ΔG ‡ values varied from 0 ± 1 to 11.6 ± 2 kcal mol −1 for CIP and from 4.7 ± 0.9 to 9.7 ± 0.8 kcal mol −1 for HRP. To the best of our knowledge, there are no published values to which we could compare these results.
The reaction free energies were calculated from CpdII 1 1 , using the relevant one-electron reduction potentials (RP) of CpdII (E CpdII 1 = 0.982 V for CIP 28 and 0.93 V for HRP 29 relevance of these reduction potentials to the ET is discussed in the SI (Page S68). One can expect the reduction potentials of enzymes and substrates to be different in the enzyme-substrate complex. However, since the reduction potentials were measured in previous work using redox titration with different substrates, we assume that these effects have already been included in the published redox potentials of both CpdIIs. The measurements and calculations of E s 1 RP and λ s for substrates are described in the SI. The results for CIP are summarised in Table 1 (HRP results are presented in the SI, Table S3) and illustrated in Fig. 1. As expected for non-homologous substrates, k 2 depends on ΔG 0 .
The reduction potentials, E s 1 , of the substrates vary from 0.162 to 0.790 V -a broad range that is unattainable by using any kind of series of homologous compounds. The self-exchange (solvent and inner) reorganisation energies, λ s , were estimated quantum chemically by using optimised structures from reduction potential calculations (cf., Table 1 with more details in the SI) [30][31][32][33][34][35][36] . The calculated values of λ s range from 3 to 12 kcal mol −1 (i.e., from 0.13 to 0.52 eV).The λ s values are comparable for structurally similar compounds; e.g., the self-exchange (solvent and inner) reorganisation energies of phenylenediamines and phenoxazines are 11 and 3 kcal mol −1 , respectively. In contrast to homologous series of substrates, our non-homologous set allows for probing a broad range of ΔG 0 and λ s .
For comparison, a previous publication, also using the framework of Marcus theory but assuming a limiting ET step, reported the reorganization energy for the oxidation of phenothiazines and phenoxazines as 7 kcal mol −1 (0.3 eV) 10 . The reorganisation energy of the overall reaction (assuming additivity) can be approximated using 18 : s s heme where λ s are the inner and solvent self-exchange reorganisation energies calculated in water and listed in Table 1.
In order to account for reorganisation energy in an enzyme-substrate complex, an additional factor, θ s , is used (see Table 1 for CIP values;values for HRP are in Table S3). θ s is calculated as the ratio of solvent-accessible surface (SAS) of the substrate in docked enzyme-substrate complexes to the SAS of free substrate (docking methodology, structures, and calculation details are provided in the SI). This relation can be justified from the large difference between static and optical dielectric constants in water and protein 37,38 . Possible uncertainties in θ s are due to uncertainties in the docking procedures. Using existing in the literature values of λ s , θ s and λ heme (0.5 eV) for cytochrome c and HRP, we obtained similar values (~0.3 eV) from Eq. 3 to those published for the overall reaction reorganisation energies for phenothiazines and pheoxazines 39,40 . Therefore, we conclude that the assumed additivity of reorganisation energies (Eq. 3) as well as the scaling of the reorganisation energy with θ s and the quantum-mechanically calculated values of λ s , accurately summarises observations made by other authors for the overall reorganization energies of oxidation reactions catalysed by CIP. In Eq. 3, there is no dependence of the reorganisation energy on the distance between electron donor and electron acceptor. The analysis of docked complexes revealed that the centres of mass of all susbtrates are at similar distances from the centres of mass ofheme. We therefore conclude that the distance-related correction factor will be similar for all compounds and will reside in theheme reorganisation energy, λ heme , which is calculated from measurements of k 2 . www.nature.com/scientificreports www.nature.com/scientificreports/

calculation of activation free energies ΔG ‡ and λ s does not support et-limited cpdii catalysis.
It has been hypothesised that the CpdII-catalysed oxidation of substrates is ET-limited [6][7][8][9]41 . We tested this hypothesis by calculating the reorganisation energies, λ, from the following form of the Marcus equation: , which should apply if the ET were the rate-limiting step in the reaction. According to Eq. 3, there should be a correlation between the reorganisation energies calculated from this equation and theoretical computations of λ s . However, we found no such correlation ( Fig. 2A), which, in the light of the similarities between the calculated and experimentally found reorganisation energies for phenoxazines and phenothiazines (vide supra), appears even more worrying.
Therefore, we hypothesise that the observed value fork 2 is actually a synthetic rate constant that accounts for both the ET and PT rates. This is supported by the linear relation between the measured free energies of activation and the logarithm of pre-exponential factors of the rate constants (Fig. 2B). It seems unlikely that such a relation would exist merely by chance (the probability of such an event for rate parameters of CIP is 1.4 × 10 −5 ). A similar correlation could be plotted from the values for formation kinetics of salts of organic amines, although the authors do not claim a linear relation of ΔG° and A 42 . The existing correlation could be explained by assuming that k 2 is a synthetic rate constant composed of both ET and PT steps. k 2 consists of pt and et rates in cpdii catalysis. There are only few alternatives to the hypothesis that ET is the limiting step in CpdII reduction (Scheme 1). If PT is the limiting step, a proton can be transferred either from the substrate or from a nearby protonated amino acid (i.e., internal PT, from HIS 55 for CIP) to the CpdII active centre (heme). However, there can be no PT from the substrate as radical cations are formed in the 9 observed reactions and substrates do not donate protons to the overall reaction (pK a values and additional discussion can be found in the SI, page S67). Instead, water self-ionisation is the proton source for HIS 55. Internal PT (from HIS 55 to heme) can be ruled out as the limiting step because the rate constants for all substrates would have to be the same. However, experimental values of k 2 vary by several orders of magnitude for different substrates; Table 1 and Fig. 1).
In the recent literature illustrated pathways in Scheme 1 are collectively termed as proton coupled ET (PCET) 43 . The term, however, was intended to be used only for the reactions where a proton and electron are transferred simultaneously, without intermediate states 44 . Nowadays, it is commonly termed as concerted proton electron transfer (CPET) 43 . The CPET mechanism can be ruled out as a mechanism due to the small kinetic isotope effect (KIE), which has previously been measured for indoleacetic acids as k h /k d ≈ 1.3 -typical KIEs are much larger in CPET reactions [45][46][47] . However, the CPET mechanism sometimes exhibits KIE values close to unity 48 ,yet theoretical predictions of the dependence of KIE on ΔG° indicate that KIE should be maximal for a zero driving force 45 . We observed an entirely different behaviour, where KIE was minimal for a near zero driving force and became maximal for a large driving force (cf., the SI section on KIE) 45 . In addition, most substrates lack a labile proton that could accompany ET in the CPET mechanism 48,49 . Therefore, the only remaining mechanism www.nature.com/scientificreports www.nature.com/scientificreports/ is a consecutive ET-PT. The ET-limited mechanism should be independent on the driving force and exhibits KIE close to 1. The PT mechanism can exhibit nearly any KIE value above 1, depending on many factors, such as barrier height, PT distance, vibronic coupling, etc. Yet, for the internal PT it should be independent on the driving force of the reaction, because previously mentioned factors should be constant. Our finding above and KIE measurements (cf. SI section on KIE) dismiss these cases. Therefore, the only alternative is ET-PT mechanism, where oxidation rate of substrates is governed both by ET and by PT rates. This mechanism can be modelled using apparent bimolecular ET and consecutive internal PT schemes that are indistinguishable under quasi-steady state conditions: The apparent bimolecular rate constant, k 2, for the total process under these quasi-steady state conditions is dependent on the concentration of substrate (see the SI for a derivation):  www.nature.com/scientificreports www.nature.com/scientificreports/ where the rate constant, k ET , is defined as in Eq. 2 (using the average initial substrate concentrations [S] for the analysis; Table S2). The pre-exponential multipliers from Eq. 2, K eq and V if , contain the dependence of electronic coupling on the distance between donor and acceptor. The equilibrium constants of the transitional complexes are combined to apparent electronic coupling ′ V if (therefore, its new dimension is cal mol −1 M −1/2 ), which is calculated for every substrate type 11,17,18 . One could argue that these factors should be estimated in order to apply ET theory. However, any uncertainty or error in estimating the equilibrium constants of the formation of transitional complexes, the dependence of electronic coupling on distance, etc., will ultimately reside in the electronic coupling factor, V if , giving a false impression of certainty in the pre-exponential factor, which is actually composed of uncertain multipliers.
The rate constant k PT is based on Marcus theory as well. The PT step in Eqs (6) and (7) is a separate internal charge-transfer reaction and depends on pK a of the proton-donor and the proton-acceptor in the CpdII active centre. Therefore, the PT driving force is constant and the free energy of activation of the PT step can be approximated with a constant, E a , as suggested by Hammes-Schiffer 49 : The substitution of Eqs (2) and (7) into Eq. (6) results in an explicit expression of k 2 . This expression in its linearized form was used to simultaneously fit the pre-exponential factors,A, and the free energies of activation, ΔG ‡ (Eqs S10 to S15 in the SI) for all substrates. In the calculations, λ s , θ s , and ΔG° are used to calculate the parameters E a , V′, λ heme , and ′ V if ( Table 2). The correlations between the calculated and measured values of A and ΔG ‡ (Fig. 3) show that Eq. 6 (involving consecutive ET and PT) is in good agreement with experimental results for both CIP and HRP. Therefore, we conclude that CpdII reduction is a consecutive ET and PT process, where PT also influences the rate constant k 2 . This conclusion is further supported by the fitted reorganisation energy of heme ( Table 2, 39 and in HRP by Andreu et al. 40 . The undefined ′ V if value of TMPD (undefined means that it was fitted as an unbound value) indicates that the rate of CpdII-catalyzed oxidation of TMPD is exclusively limited by a PT step.
In contrast to the widely held assumption regarding similarity of V if for homologous compounds (which, to the best of our knowledge, has never been discussed, yet is frequently used in applying Marcus ET theory to series of homologous substrates [6][7][8][9][10][19][20][21][22][23][24][25]50 Table 2) varies by up to 5 times between homologous compounds (e.g., TH and MB, or CPZ and PZ). Therefore, traditional investigations of the dependence of enzymatic ET rate constants on potentials of homologous substrates only partially reveal the trends predicted by Marcus theory. We measured the rate constants at different temperatures and the results suggested a new way to estimate ′ V if . This approach is general and does not depend on the assumption of similar ′ V if values for a set of homologous substrates. The proposed model (Eq. 6) splits the consecutive ET-PT process into separate steps and allows estimating of the Marcus free energy of activation of ET Δ ‡ G ET for each compound. Δ ‡ G ET is small and ranges from 0 to 5 kcal/mol. The fitted E a of the PT step is 11.7 ± 0.3 kcal/mol and the rate constant k PT = 520 ± 40 s −1 at 298 K for CIP. For HRP, this becomes 10.7 ± 0.1 kcal/mol with a rate constant of k PT = 200 ± 30 s −1 at 298 K. Thus, the PT step is the main factor controlling ΔG ‡ . In this context the Marcus quadratic dependence of activation free energy on the driving force is only a first approximation of a much more complicated ET. Based on the driving forces, four of the www.nature.com/scientificreports www.nature.com/scientificreports/ investigated compounds appear to fall into the inverted region, where the dependence of activation energy on driving force becomes asymmetric. In such a case, the original Marcus law greatly overestimates the activation energy of the reaction. This asymmetry is more pronounced when the range of driving forces is greater than 2 eV and reactants change their structures appreciably upon ET (which is not applicable here) [51][52][53][54] . Then, the vibrational modes start to participate in ET and lower the activation energies. This would affect the compounds VB, TMPD, MB, and TH and could lower their activation free energies of ET by up to several kilocalories per mole (which is within the margins of measurement error). In this context, the experimental free energy of activation for TMPD (11.6 ± 2 kcal/mol) appears high and can only be explained by invoking ET-PT or PT-ET mechanism rather than by pure ET step. physical constraints of pt reveal that et precedes pt step in cpdii catalysed oxidation. The calculated values of PT rate constants and E a can be directly related to known structures of peroxidases. It is commonly accepted that in CpdII reduction to ferric peroxidase, the proton is transferred from a protonated HIS 55 (in case of CIP) to a heme-bound oxygen atom 5,[55][56][57][58] . This step can precede or follow the ET. We assume that the main difference between the ET-PT and PT-ET mechanisms comes from the distances involved in PT. These distances are controlled by the oxidation state of the heme iron and the position of bound oxygen, either in the form of a hydroxyferryl or hydroxide complex. The geometry of ferric (reduced) CIP is known (PDBID:1h3j), and its active centre is illustrated in Fig. 4 59 . If ET precedes PT, the positions of heavy atoms before and after the PT step should not change significantly. Therefore, the distance of PT (1.0 ± 0.2 Å) can be estimated by subtracting the bond lengths of N-H and O-H from the distance between the HIS 55 nitrogen and the oxygen atom (3 Å) (Fig. 4). Similarly, by assuming an iron -oxygen distance of about 1.7 Å for CpdII 60 , we estimate the PT distance in the PT-ET mechanism as about 1.7 Å. The PT rate constants are calculated using the multi-state non-adiabatic proton theory by Borgis and Hynes (seethe SI) 61 .
E a of the PT step (11.2 ± 0.3 kcal/mol) and k PT (520 ± 40 s −1 at 298 K) over a distance of 1.7 Å (PT-ET) can be reproduced with λ = . For PT-ET, pK a of the proton acceptor in heme should be about 42, which would be too large. Published pK a values for CpdII are about 7 for the histidine side-chain and about 9 for heme-bound water 64,65 . The reduction of CpdII should therefore proceed via the ET-PT mechanism for both CIP and HRP (as both have similar catalysis centres and rate parameters). CpdI and CpdII compounds are quite similar in nature, therefore, histidine side-chain pK a values of CpdI can partially reflect pK a of CpdII. In the literature these values cover a range from 4 to 8 [66][67][68] . However, even taken this information into account, ET-PT mechanism is still much more reasonable compared to PT-ET.
Kie estimations for cip lend support to the et-pt hypothesis. If we assume that ΔG 0 , λ p, and the deuteron transfer (DT) distance are the same as in the ET-PT mechanism, we obtain a DT rate constant of 295 s −1 and a KIE of 1.8 ± 0.2 for the PT step for CIP. The KIE of the full ET-PT process differs for each substrate (Table S2). The lowest KIE value (about 1) is predicted for CPZ, while the largest value (about 1.8) is calculated for phenylenediamines TMPD and AMB. We found the average KIE value to be about 1.4 which appears in good agreement with previously published values for indoleacetic acids 6 . The KIE values for all three CpdII-catalysed www.nature.com/scientificreports www.nature.com/scientificreports/ oxidation reactions and substrates TMPD, CPZ, and ABTS for CIP were measured to be 1.6 ± 0.2, 1.0 ± 0.1, and 1.3 ± 0.1, respectively (for HRP-CpdII we found a similar behaviour; see the SI). They are in good agreement with calculated values (cf. values in the Supporting Table S2) and further support our hypothesis that ET is followed by the PT mechanism during the oxidation of substrates catalysed by CpdII.

conclusions
We demonstrated a general approach for the investigation of electron and proton transfers in oxidoreductase-catalysed reactions. The measured dependence of the reaction rate constants on temperature (i.e., Arrhenius plots) yielded the reaction-free energies of activation as well as the pre-exponential factors needed for the rigorous comparison of theoretical predictions and experimental observations. The approach was applied to the oxidation of non-homologous substrates catalysed by Coprinopsis cinerea and Armoracia rusticana (horseradish) peroxidases' compound II. By combining the above approach with quantum chemical calculations and docking studies we were able to compare the experiments with theory which led to new insights about the catalysis mechanism of CpdII. Contrary to commonly held assumptions 8,9,19,[22][23][24][25][26][27] , the apparent electronic coupling constant ′ V if for different series of homologous substrates differ considerably. Moreover, ET is not the rate-limiting step for CpdII-catalysed oxidation. It was found that a PT step follows the initial ET step described by Marcus theory. We could explain the value of the PT rate constant based on molecular geometry and solvent reorganisation energies of PT using the harmonic oscillator approximation and non-adiabatic proton transfer theory. The proposed model correctly predicts the observed kinetic isotope effects and further supports that a PT follows ET in CpdII catalysis.

Methods
Reaction rate constants and free energies of activation. All experiments were carried out in 50 mM phosphate buffer solution, with a pH of 7.00, using a thermostated stopped-flow spectrophotometer (Otsuka RA-401) cell at temperatures ranging from 10 to 30 °C. Kinetic curves were measured for five different initial concentrations for each substrate reaction, and experiments were repeated 5-10 times. Further details on the materials used and the experiments are available in the SI.

electrochemical measurements. Electrochemical experiments were performed with a Series G 750
Potentiostat/Galvanostat/ZRA from Gamry Instruments, Inc. (Warminster, PA, USA), controlled by dedicated PHE200 software. All measurements were carried out in a three-electrode glass cell using an Ag|AgCl|3 M KCl reference electrode (210 mV vs. NHE) and a platinum wire counter electrode from BASi (West Lafayette, IN, USA). All electrochemical measurements were performed in 50 mM phosphate buffer solution, with a pH of 7.00, using glassy carbon (3 mm diameter) or gold (1.6 mm diameter) working electrodes, obtained from BASi. Prior to any measurement, the electrodes were polished with diamond paste (particle size of 1 µM, BASi), alumina slurry (particle size of 0.3 µM, BASi) and then sonicated. Gold electrodes were also cleaned electrochemically by running 40 potential sweep cycles from 0 to 1.9 V vs. NHE in 0.5 M H 2 SO 4 solution. Additional experimental details are presented in the SI.
Structural and standard reduction potential calculations. All quantum chemical calculations were performed using GAMESS 32 . Specifically, the molecular structures were first optimized by using the MINI basis set and then used as a starting structure for HF/6-31 + G(d,p) energy optimization. After optimization, the Hessians were evaluated and all structures were determined to be at their energetic minima with no imaginary frequencies. Final optimizations were carried out using DFT 33 with the hybrid B3LYP functional 34 and the 6-31 + G(d,p) basis set. Energies of solvated species were evaluated at the same level. Restricted open-shell Hartree-Fock wave function (ROHF) was used for radicals; for other structures, the restricted Hartree-Fock wave functions (RHF) were utilized. All calculations were performed using the solvation model density solvent model method 35 as implemented in GAMESS. Zero point energies were calculated at the HF/6-31 + G(d,p) level. Further details about the calculation are available in the SI.
Docking studies and calculations of solvent-accessible surface areas. For docking studies, the rCip structure (PDBID:1h3j) 69 and the HRP (PDBID:1h55) 5 were used. Crystallized water molecules were removed from the peroxidase structure and the oxygen-iron distance was reduced in order to reflect the existence of the compound II state rather than the reduced state of the peroxidase active centre containing complexed water. The docking of substrates was performed using AutodockVina 70 . Solvent accessible substrate surface areas in the substrate-enzyme complexes were estimated from the docked structures that exhibited the best energy score when using the MSMS method 71 , implemented into USCF Chimera 72 . The docking results are available in the SI.