Extreme Catalytic Power of Ketosteroid Isomerase Related to the Reversal of Proton Dislocations in Hydrogen-Bond Network.

Dynamic electrostatic catalytic field (DECF) vectors derived from transition state and reactant wavefunctions for the two-step reaction occurring within ketosteroid isomerase (KSI) have been calculated using MP2/aug-cc-pVTZ and lower theory levels to determine the magnitude of the catalytic effect and the optimal directions of proton transfers in the KSI hydrogen-bond network. The most surprising and meaningful finding is that the KSI catalytic activity is enhanced by proton dislocations proceeding in opposite directions for each of the two consecutive reaction steps in the same hydrogen network. Such a mechanism allows an ultrafast switching of the catalytic proton wire environment, possibly related to the exceptionally high KSI catalytic power.

Introduction

Proton transfer between hydrogen-bonded biomolecules is one of the most important phenomena in biology. In particular, recent studies indicate that the catalytic activities of numerous enzymes are closely related to the proton transfer in hydrogen-bond chains extending from the active site.[1−8] For example, recent experimental evidence indicates that change in the proton position in the second shell residues forming the hydrogen-bond chain enhances the rate of catalysis in aspartate aminotransferase.[1] Due to the large size of the involved molecular systems and usually insufficient information regarding the positions of hydrogen atoms, the corresponding molecular mechanisms could not be easily explored by conventional experimental or theoretical methods. Therefore, for the purpose of this work, we use a more cautious term “proton dislocation”, which is more adequate for short and strong hydrogen bonds, where the degree of partial proton transfer is usually not known precisely.

The aim of this contribution is to apply for the first time a method allowing the identification of the magnitude of dynamic catalytic effects and the preferred proton dislocation directions in hydrogen-bond chains extending from reactants and leading to a reduction in the activation barrier. For this purpose, we decided to test a new application of dynamic electrostatic catalytic fields (called thereafter as DECFs), defined here as the difference of the transition-state and substrate electric fields at the midpoint of every hydrogen bond, which constitutes catalytic environment. Although the theoretical concept of DECF had been originally introduced earlier,[9,10] it has not been practically applied to any enzyme reaction yet. This was due to considerable difference between the typical enzyme reaction timescale and much longer time required to dislocate active site residues or their side chains containing heavy atoms. However, this limitation disappears in the case of ketosteroid isomerase (KSI), where rapidly oscillating very light protons located in a relatively rigid hydrogen-bond network surrounding an active site have a chance to be coupled with the enzyme reaction coordinate. Due to the existing extensive experimental and theoretical data indicating an important role of the HB network,[11−13] KSI has been selected to test the DECF approach for the first time for an enzyme system. Whenever a strong DECF is detected, it could be regarded as the driving force for catalytic proton dislocation, possibly propagating in the neighboring HBs due to the domino effect. This way, our approach extends and supplements previous attempts of mapping proton wires in proteins.[2] However, catalytic fields are derived using a bottom-up approach, where the analysis of the charge redistribution from the reactants to the transition state yields properties of an optimal catalytic environment in a compact form,[9,10,14−18] constituting an inverse solution to the optimal catalyst problem. The only input data required to predict favorable proton-transfer directions for different reaction stages are the reaction pathway stationary points structures and the positions of electronegative atoms of the molecular environment only, which could serve as hydrogen bond proton donors or acceptors. No information about the usually doubtful positions of the hydrogen atom is necessary.

Methods

The catalytic activity could be expressed in terms of perturbation theory[8,14] as the difference of the transition state TS and substrate RS interaction energies ΔE in a catalytic environment C, constituting differential transition state stabilization energy EDTSS approximating activation barrier loweringA detailed analysis of the electrostatic, exchange, delocalization (induction), and correlation (dispersion) EDTSS constituents for chorismate mutase,[10] cAMP-dependent protein kinase,[15] ribonuclease,[16] several KSI mutants,[17] and Si → Al substitutions in zeolites[18] indicated the dominant role of the electrostatic term EEL, which yields the results significantly correlated with higher theoretical (RHF, MP2) levelsEDTSS can be further approximated (1) by the sum of the products of catalyst atomic charges qi and the difference in the molecular electrostatic potentials of the transition state VTS and substrate VRS. This difference defines the static catalytic field ΔS(i), which represents the activation energy lowering resulting from the presence of a unit point charge +1 at any point “i”The highest catalytic activity is achieved when charges representing a static molecular environment coincide with extremal −ΔS values. Due to the additive nature of electrostatic interactions, −ΔS represents the solution of the optimal catalyst problem, allowing an inverse design of the new catalysts and the rapid activity estimates of mutated enzymes.[17,19] This constitutes a bottom-up approach of catalyst design based on the knowledge of reactant charge redistribution during activation at a given reaction step.[19] It is fundamentally different from the conventional top-down methods requiring considering entire protein and numerous assumptions related to quantum mechanical/molecular mechanical (QM/MM) boundaries, selection of protonation states, and arbitrary empirical force fields.

As the rigid structure of the catalyst C is assumed to be derived from the static catalytic fields ΔS, the coupling of reaction coordinates with much slower protein motions cannot be then considered. However, it is well known that rate-promoting residue vibrations could be detected in some enzymes[20] and certainly rapid oscillations of protons in hydrogen bonds could couple directly with reaction coordinates. Let us define, for a given hydrogen bond between atoms X and Y, the quantity as a negative gradient vector field ΔS calculated along the line connecting X and Y at its middle point (eq ). Such a one-dimensional dynamic electrostatic catalytic field (DECF) would therefore indicate the direction of the favorable unit probe charge movement with respect to reactants, resulting in a further decrease of the activation barrier (this idea has been illustrated in Figure 1 in ref (9) and expressed in eqs 4 and 5 in ref (9) and eq 16 in ref[10]). In general, the three-dimensional −grad ΔS could be oriented in any direction and the DECF magnitude is the fraction of its length. In other words, vector indicates the optimal direction of proton movement in a linear hydrogen bond, leading to an increase in catalytic activity. In case of the same coordinates of atoms X and Y in both RS and TS structures, ΔD can be rigorously calculated as the difference of the electric field vectors and located in the middle of hydrogen bond X–H···Y, projected on the X···Y bond axis described by and position vectorswhere the dot denotes a scalar product. In the following text, is the vector quantity of the dynamic catalytic field and ΔD, as the signed scalar quantity, is the magnitude of with a positive sign if the vector is facing from X to Y and negative otherwise (i.e., the first factor of the product in eq ). Vectors and indicate the positions of X and Y atoms forming the X–H···Y hydrogen bond. Electric field vectors and are calculated at the midpoint of the hydrogen bond and then projected on the X–H···Y axis.

Equation can be directly applied when the coordinates and are the same in TS and RS structures, i.e., when the proton movement is the only catalytic effect. In practice, to make independent of the orientation of the TS and RS coordinates and to allow the application of the method to an ensemble of microstates of the catalytic complex, it can be calculated as the difference of the electric field vector projections on the X···Y axes, calculated independently for TS and RS geometriesIt is worth noticing that the change in X and Y coordinates may be an artifact of the computations; even when the proton dislocation is the only catalytic effect, the orientation and position of X and Y atom pair may change due to the independent optimization of each stationary structure.

Scalar values of ΔS and vectors represent in compact form the electrostatic effects of outside reactants due to their charge redistribution on the reaction pathway from substrates to transition state.

Calculations include the following steps:

The transition state and substrate stationary reaction structures are obtained using the standard QM/MM or ONIOM technique.[21] The QM part should include any active site residue forming covalent bonds with the reactants in the reaction course (ASP38 in KSI case). In this work, stationary points obtained in the B3LYP/6-31+G*/CHARMM27 QM/MM approach[17,22,23] have been used.

Coordinates of the entire examined protein including reaction pathway stationary points in PDB format are used as an input by the locally implemented Python script to locate all electronegative atoms located within the assumed distance from reactants (12 Å from O3 in this work) capable to form hydrogen-bonded linear and branched chains separated by a specified distance range (D = 0–3.1 Å in this work). In the case of a multistep reaction, like in KSI, it is recommended to use stationary point coordinates obtained in the same reaction coordinate scan obtained for a given reaction step and use MM minimized enzyme structure with a substrate of every step.

For every hydrogen bond included in the HB chains determined earlier, expectation values of electric field vectors and are calculated at every HB midpoint, directly from the corresponding wavefunctions of reaction stationary points. Any standard quantum chemical program, like Gaussian[24] or GAMESS,[25] could be used for this purpose.

The electric field vectors and are individually projected on each HB axis and then subtracted to obtain the DECF vector and its one-dimensional projection on the bond ΔD, according to eq .

A list of all hydrogen bonds forming linear or branched chains and characterized by the same ΔD sign and DECF value above assumed threshold (0.001 au in this work) is calculated.

To visualize the results in the form of vectors, the program VMD (visual molecular dynamics)[26] was used.

Results and Discussion

To validate our approach, we selected ketosteroid isomerase enzyme, in which proton transfer between TYR14, TYR55, and ASP99 plays a crucial role in their exceptionally high catalytic activity and has been extensively studied experimentally and theoretically.[11−13,22,23,27,28] Isomerization of 5-androstene-3,7-dione (RS) to 4-androstene-3,7-dione (PS) shown in Figure proceeds in two steps via the first transition state TS1, the intermediate complex IS, and the second transition state TS2.[17,22,23] In this case, we could examine the possible catalytic role of various KSI HB chains for each reaction step, i.e., RS-TS1 and IS-TS2 separately. All stationary points have been obtained within the B3LYP/6-31+G*/CHARMM27 QM/MM setup.[17,22,23] In an analogy to our previous study related to KSI mutants,[17] the QM part included, besides reactants, the side chain of ASP38 residue directly involved in the isomerization reaction. Using MP2/aug-cc-pVTZ results presented in Table as the reference, the maximum error of DECF values estimated at MP2/aug-cc-pVDZ theory level does not exceed 0.001 au, whereas same qualitative results could be still obtained from RHF/aug-cc-pVTZ or RHF/aug-cc-pVDZ (0.0015 au max error), RHF/6-311G** (0.0035 au max error), and B3LYP/6-31+G* (0.003 au max error) with a much smaller computational effort. Complete results obtained using all of the above-mentioned theory levels are presented in Table S1 in the Supporting Information along with hydrogen bond lengths. These results could be useful in selecting appropriate theory level in the case of large size reactants requiring more computational resources.

Figure 1

Mechanism of isomerization of 5-androstene-3,7-dione (RS) to 4-androstene-3,7-dione (PS) catalyzed by ketosteroid isomerase. QM part adapted from refs (17) and (21) included two first rings of the steroid and the side chain of Asp38 residue, which is involved in reactant hydrogen transfer from C4 to C6 position. This scheme constitutes a modified version of Figure 1 from our previous work (ref (17)https://pubs.acs.org/doi/10.1021/acs.jctc.6b01131) and further permissions related to the material excerpted should be directed to the ACS.

Table 1

Components of Dynamic Electrostatic Catalytic Field ΔD Values in [au] along Hydrogen Bonds in Branched HB Chains Involving TYR14, TYR55, ASP38 Residues, and Water W1, W2, W3 Molecules for RS-TS1 and IS-TS2 Reaction Steps in Ketosteroid Isomerase Enzymea

	ΔD/∑ΔD [%]		ΔD [·10–3 au]		cumulated ∑ΔD [·10–3 au]
hydrogen bond	RS-TS1	IS-TS2	RS-TS1	IS-TS2	RS-TS1	IS-TS2
(TYR14)O···H→···O3(AND)	36.1	29.4	10.1	–8.52	27.9	–29.0
(ASP99)O···H→···O3(AND)	21.6	26.0	6.03	–7.53
(TYR55)O···H→···O(TYR14)	10.1	9.5	2.81	–2.78
(ASP38)O···H→···O(W1)	14.5	19.3	4.00	–5.62
(W1)O···H→···O(W2)	9.6	8.0	2.70	–2.31
(W1)O···H→···O(W3)	8.3	7.8	2.31	–2.27

AND stands for substrate molecule (5-androstene-3,7-dione). All results have been obtained at MP2/aug-cc-pVTZ theory level.

Initially, 43 hydrogen-bonded chains (D < 3.1 Å) have been detected in KSI without considering water and 81 when water at the distance of 12 Å from O3 is included. The chains of hydrogen bonds extending from substrate and active site Asp38, including water up to 12 Å from substrate O3 oxygen, were selected for the calculation of . The first reaction step RS-TS1 related to the abstraction of hydrogen from the substrate C4 atom is strongly enhanced by the unidirectional proton movement from TYR14 and ASP99 to substrate oxygen O3 and from ASP38 to water W1, with additional unidirectional proton dislocations in TYR55, water W2, and W3 forming branched HB chains (Figure a). The second critical reaction step IS-TS2 characterized by the higher activation barrier[17] (Figure b) is related to attaching hydrogen to the substrate carbon C6. The second reaction step IS-TS2 is strongly catalyzed by proton dislocations in the reverse directions in the same HB chains, with the dominant role of TYR14 and somewhat smaller role of water W1 and ASP38. The corresponding ΔD and ∑ΔD values are presented in Table . These results are in agreement with the recent experimental and theoretical results indicating that 98% changes of electric fields in the KSI active site originate from TYR14 and ASP99 residues.[11−13,22,23,27] Catalytic field directions obtained for the second reaction step IS → TS2 are opposite to those of the first step RS → TS1; however, the magnitudes of ΔD and their relative ΔD/∑ΔD contributions remain very similar. This indicates that proton dislocation assisting in consecutive reaction stages RS → TS1 and IS → TS2 have to be only reversed within the same rigid hydrogen bond scaffold. Hydrogen bond lengths given in Table S1 in Supporting Information for RS, TS1, IS, and TS2 stationary points change less than 2% (0.076 A). This particular feature could be very meaningful and possibly responsible for very high KSI activity. After the completion of the first reaction step RS → TS1, minimal, possibly femtosecond, time would be required to switch the direction of proton dislocation assisting the second reaction step IS → TS2. To the best of our knowledge, such a finding has not been reported in literature yet. Similar magnitudes of ΔD values for both reaction steps are remarkable but until more other systems are examined, we cannot definitely prove that this is not accidental. Possibly, this reflects a similar orientation of particular HB with respect to both local reaction centers at C4 and C6, as ΔD values are more symmetric for remote HBs.

Figure 2

Directions of proton dislocations along HB chains extending from the KSI active site with relative magnitudes of catalytic effects indicated by the corresponding dynamic catalytic field vectors in [10–3 au] obtained using the MP2/aug-cc-pVTZ approach for (a) the first reaction step RS → TS1 (green arrows), (b) the second reaction step IS → TS2 (violet arrows), and (c) KSI reaction profile.[17] The part of the system included in QM calculations[17,22] presented using thicker ball-and-sticks representation.

In general, computational analysis of strongly correlated multiple proton-transfer processes in hydrogen-bond chains would be extremely demanding due to its multidimensional character. On the other hand, a single proton-transfer process in the hydrogen bond characterized by the largest dynamic catalytic field ΔD value will certainly induce additional proton dislocation in the neighboring hydrogen bonds and could constitute the driving force for a collective PT. Therefore, we consider the sum of individual DECF projections along HB axes ∑ΔD as the quantitative measure of the catalytic activity of the entire chain and possibly connected branches. The largest relative DECF contributions ΔD/∑ΔD obtained for TYR14···AND (36.1 and 29.4%), ASP99···AND (21.6, 26.0%), and ASP38···W1 (14%, 19.3%) are in agreement with the recent path integral ab initio molecular dynamics results indicating the possibility of partial ionization of KSI ASP38, TYR14, and ASP99 residues and their dominant contribution to electric fields in the active site.[11−13,28] These conclusions are supplemented by the recent computational study of the KSI product release indicating the essential role of the ASP99 side chain.[28]

The importance of extended HB networks in KSI has been confirmed by extensive ab initio path integral molecular dynamics studies.[13] Recent independent analysis using quantum theory of atoms in molecules (QTAIM) reported a significant correlation between charge density, external electric field, and activation energy changes in KSI.[29] These results obtained by entirely different methodologies are complementary with ours reinforcing conclusions in both studies. This is in line with the well-known proton-coupled electron-transfer effects.[30]

Whereas earlier top-down computational simulations[13,28,31] focused on electric fields exerted by KSI residues on reactants, in our bottom-up approach we investigated the electric field changes, resulting from the reaction progress, acting on protons in highly polarizable HBs outside reactants. Such an inverse approach allows obtaining more detailed insight into the dynamics of the catalytic environment, which cannot be examined by the existing conventional theoretical or experimental methods.

As most of the HBs considered in this work are extremely short, possibly with a single flat minimum, the actual proton dislocation may have a partial character observed in highly polarizable HBs,[32] and in KSI in particular.[11,12]

Conclusions

Dynamic electrostatic catalytic fields could provide new insight into the mechanisms of enzyme reactions involving extensive hydrogen-bonded networks in relatively rigid proteins with high catalytic activity, like carbonic anhydrase or green fluorescent protein.[31] The ability to reverse proton-transfer directions could also constitute an important principle in the design of theozymes aimed to catalyze multistep reactions. This approach could be also useful in scanning enzyme hydrogen-bond networks for more accurate theoretical or experimental studies related to a controversial idea of low barrier hydrogen bonds (LBHBs)[5,7] or postulated covalent nature of transition state interactions with most active enzymes, usually via hydrogen bonds.[33]