Coupling Protein Dynamics with Proton Transport in Human Carbonic Anhydrase II.

The role of protein dynamics in enzyme catalysis is one of the most highly debated topics in enzymology. The main controversy centers around what may be defined as functionally significant conformational fluctuations and how, if at all, these fluctuations couple to enzyme catalyzed events. To shed light on this debate, the conformational dynamics along the transition path surmounting the highest free energy barrier have been herein investigated for the rate limiting proton transport event in human carbonic anhydrase (HCA) II. Special attention has been placed on whether the motion of an excess proton is correlated with fluctuations in the surrounding protein and solvent matrix, which may be rare on the picosecond and subpicosecond time scales of molecular motions. It is found that several active site residues, which do not directly participate in the proton transport event, have a significant impact on the dynamics of the excess proton. These secondary participants are shown to strongly influence the active site environment, resulting in the creation of water clusters that are conducive to fast, moderately slow, or slow proton transport events. The identification and characterization of these secondary participants illuminates the role of protein dynamics in the catalytic efficiency of HCA II.

Introduction

Although proton transfer/transport (PT) in simple solvents has been studied for several decades, investigations into the molecular mechanism of enzyme facilitated PT events have only witnessed an intense surge of activity in the past few years.[1−14] It is now well understood that an excess proton can be transferred through an enzyme along dynamic pathways comprised of hydrogen bonded networks of amino acid residues and water molecules.[15,16] The inherent multiscale complexity of enzymes[17−28] can render fluctuations of such networks to be nontrivially correlated, over both space and time, to the motion of the excess proton. Thus, modeling of key PT steps in enzymes becomes prohibitively difficult and requires input from both theory and experiments at various levels of detail (temporal and spatial). This is well demonstrated by human carbonic anhydrase II (HCA II), which has long served as a prototype for enzymes containing a rate determining PT step.[29−49] HCA II is not only one of the fastest known enzymes, it is also one of the most extensively studied systems with a wide range of structural, mutation, and kinetic experiments, supplemented by a large number of theoretical and computer simulation studies. Based on these investigations, detailed insight into the underlying mechanism of PT in HCA II has been presented.[39,50−53] The focus of this article is on understanding if the motion of an excess proton is correlated to fluctuations in the surrounding protein and solvent matrix that may be rare on the pico- and subpicosecond time scales of molecular motions.

HCA II is a zinc metalloenzyme that catalyzes the reversible hydration of carbon dioxide with a maximal rate constant value of 0.8 μs–1.[31,39] The HCA II mediated catalysis is widely believed to proceed through two steps:[54,55] viz., (1) formation of bicarbonate resulting from a nucleophilic attack by a zinc-bound hydroxide on carbon dioxide, and (2) intramolecular transport of a proton from the zinc-bound water molecule through a hydrogen-bonded water cluster to a histidine residue (His64) that is 8 to 10 Å away, and near the mouth of the active site cavity. Under physiological conditions, when the exogenous buffer concentration is in excess, the second intramolecular PT step is widely believed to be rate determining.[56] However, when the exogenous buffer concentration is limited, an intermolecular PT between His64 and the buffer becomes rate determining.[39,57,58] Mutation of His64 to alanine, a nonproton conducting residue, has been reported to lower the rate of intramolecular PT by ∼20-fold.[59] This observation suggests that His64 acts as the dominant proton donor/acceptor during catalysis, but that alternate pathway(s) are also capable of rapid proton release. Moreover, the kcat/KM (i.e., catalytic efficiency) can be significantly impacted by various mutations of Trp5, Tyr7, Asn62, and Gln92, suggesting that these amino acids also play a role in the creation and viability of PT pathways.[34,60,61]

Two mechanistic components clearly contribute to the efficient PT, and hence fast catalytic rate, in HCA II: the orientation of His64[62,63] and the formation of water clusters between the zinc-bound water and His64 side chain.[31,36,39,40,63−65] Detailed analysis of the high resolution X-ray crystal structure of the wild type enzyme indicates the presence of two conformers of His64 with different side chain orientations.[37,66,67] The inward conformer corresponds to an rZn–Nδ1 separation of about 8 Å with the imidazole ring of the His64 side chain buried inside the active site cavity. The outward conformer, directed away from the active site, is characterized by an rZn–Nδ1 separation of about 10 Å.[67] As reported by molecular dynamics (MD) simulations,[62,63] the relative populations of these conformers and free energies of their interconversion depend on the protonation state of the imidazole ring (see Table S1).

Due to the lack of other protonatable residues between the donor and acceptor, water molecules are hypothetically essential for PT (Figure ). This hypothesis is supported by neutron structure data,[47] which reveals the presence of a hydrogen bonded water cluster in the active site of HCA II that bridges the zinc-bound water to the inward conformer of His64.[43,63] Extensive MD simulations have been performed with the proton residing either on the zinc-bound water or on His64 and averaged over long-time dynamic fluctuations.[43] When His 64 is oriented inward, an efficient PT path is created via 3–4 water molecules that are stable up to a picosecond. In comparison, when His64 is oriented outward, larger clusters containing 5–6 water molecules may be accommodated, with only transient formation of a complete bridging path and having lifetimes of about 0.2 ps. Computational[44,51,65,68] and experimental[57,59,69−71] studies on a large number of single and double mutants of HCA II also suggest that the catalytic efficiency depends on the formation of less branched clusters involving 2–3 water molecules and on the presence and orientation of His64. It has also been found that fluctuations in the size of active site water clusters are significantly correlated to changes in the orientation of His64.[43] In particular, the entry and exit of the side chain from the transition state region delineating the inward and outward conformations has been reported to be associated with a change in the number of water molecules in the active site within about 5–10 fs.[40] Thus, the rate determining PT step appears to be optimized by the conformational switching of His64 and associated dynamic reorganization of the bridging water clusters between the zinc-bound water/hydroxide and His64.

Figure 1

Hydrogen bonded water network in the active site of HCA II depicting the conserved water cluster connecting the catalytic zinc and the purported proton acceptor sink, His64.

The orientational preference of His64 and the infrequent transitions between stable conformations are hypothesized to be coupled also to the dynamics of other residues present in the enzyme. As shown with Monte Carlo sampling,[63] the outward to inward rotation along a minimum-energy path requires His64 to pass through a narrow channel-like region bordered by Asn62 and Trp5. Experimental studies[72] have also shown that Asn62 maximizes the PT efficiency by creating a favorable environment for the two conformers of His64. A subsequent transition path sampling simulation study[40] showed that residues Asn62, Trp5, and Tyr7 rearrange as the His64 side chain passes through the transition state between its outward and inward conformers. This study further demonstrated how the coupled fluctuations of Asn62 and Asn67 instigate the inward rotation of His64 from its outward to inward conformer.[40]

Although interesting, these classical simulations were not designed to address the question of what protein motions enable the actual PT event in and out of the active site cavity. In contrast, the multistate empirical valence bond (MS-EVB) method[39] was utilized to examine the mechanism of PT in HCAII, including the correlation of His64’s orientation with the distribution of proton pathways through the active site. The following three models were investigated (see Table S2 for relative PT free energy barriers):

System 1: His64 in its inward orientation and included in the MS-EVB basis set.

System 2: His64 in its outward orientation and excluded from the MS-EVB basis set

System 3: His64 in its outward orientation and included in the MS-EVB basis set

System 2 represented an important baseline model for which the estimated free energy barrier was 14.6 ± 0.4 kcal/mol, corresponding to a rate constant of 0.0004 μs–1, which is substantially slower than the experimentally observed rate constants of 0.8 and 0.02 μs–1 for the wild type and His64Ala mutant, respectively.[39] The mechanism of PT was found to involve the following: (i) a weakly stabilized contact ion pair between the zinc-bound hydroxide and the hydronium cation around the position of water W1 (Figure ) where the character of the hydrated excess proton is predominantly Zundel-like; (ii) Eigen-like localization of the center of excess charge (CEC) on W1; (iii) transfer of the excess proton to W2, where the Zundel-like character persists; and (iv) a high free energy Eigen-like cation around the location of W3b (marking the free energy barrier for PT to bulk), which is in close proximity of the residues Asn62 and Asn67. Destabilization and increased mobility of the active site waters occupying locations W2, W3a, and W3b, as observed in the mutant His64Ala of HCA II, are associated with a further increase in the free energy barrier in the absence of any ionizable residues or chemical rescue agents.[48]

Building on this baseline study, Maupin et al. showed that the free energy barrier remains around the location of W3 (a or b) when His64 is included in the MS-EVB basis state (i.e., participates explicitly in PT). The presence of His64 as the proton acceptor was found to lower the barrier for PT to 10.0 ± 0.4 kcal·mol–1 when His64 was in its inward conformation (System 1), and to 11.4 ± 0.4 kcal·mol–1 when His64 was in its outward conformation (System 3).[39] Therefore, it is hypothesized that in wild type HCA II the most favorable pathway corresponds to His64 accepting the excess proton from the active site while it occupies the inward conformation. His64 then rotates to the outward orientation where it can release the excess proton to the surrounding bulk environment. A secondary pathway corresponding to His64 accepting the excess proton from the active site while it occupies the outward conformation is also probable, but this pathway is associated with a higher free energy barrier.

It is evident from the above discussion that the free energy barriers encountered by the excess proton after deprotonation of the zinc-bound water, peak during the transition from a Zundel-like to an Eigen-like cation near the location of W3. More specifically, the free energy maximum corresponds to an Eigen-like cation centered primarily on W3b when His64 is in the outward conformation and on W3a when His64 is in the inward conformation. Interestingly, the localization of the excess proton near the sites of W3a or W3b brings it in close proximity to residues such as Asn62, Asn67, and Tyr7 (Figure ), which are intrinsically coupled to the motions of His64 and, therefore, hypothesized to influence the rate limiting PT event.[40,46,63] In this article, our main objective is to probe the dynamics across the PT free energy barrier, and hence the influence of protein motions on the mechanism of PT.

We have focused our present investigations on these three systems. First, the fastest uptake of the excess proton by His64 side chain in its inward orientation (system 1) is studied using a combination of MS-EVB[2,3,39,48,73] and forward flux transition path sampling simulations[74−78] to reveal the active site structures that facilitate the fastest PT events. Second, alternative PT pathways are investigated with His64 in its outward orientation (systems 2 and 3). Schematic representations of these model systems are shown in Figure . Both Systems 2 and 3 give much slower PT events, allowing us to examine the role of slow protein motions on PT. In system 2, by additionally eliminating His64 from the MS-EVB state space (i.e., not allowing it to participate in PT), we investigate how the excess proton uses alternate water-mediated paths to enter/exit the active site, without utilizing His64. In System 3, these alternative paths relay the excess proton to the edge of the active site for its subsequent uptake by His64 in its outward orientation. Systems 2 and 3 further probe the dynamical importance of His64 in the wild type enzyme, in which His64 does participate in the proton relay. These systems are also helpful to understanding mutant systems, such as His64Ala. We have carried out an extensive transition path sampling study of System 2 and 3 in combination with MS-EVB simulations and analyzed the mechanism(s) of PT across the corresponding free energy barriers. A comparison of the dynamical trajectories of these models is found to elucidate the role of several active site residues on the mechanism of the rate determining PT step in HCA II, and illuminate the coupling of protein dynamics and catalysis. A nontrivial coupling is detected between the rate of PT to His64 and the side chain orientations of residues such as Asn62 and Gln92 that may slow down the transfer of the excess proton by nearly 3 orders of magnitude.

Figure 2

Schematic representation of PT pathways investigated by transition path sampling simulations in model systems 1 (left) and 2,3 (right). Black circles indicate typical locations of the excess charge when the respective system resides at the free energy barrier region. Schematic transition paths originating from these barrier regions are shown using curved arrows delivering the excess charge to His64 in its in conformation (purple, system 1), to the zinc-bound hydroxide ion (pink, system 2), and to His64 in its out conformation (blue, system 3). In systems 2 and 3, the entry into the active site is primarily observed through the region between Asn62 and Asn67 (red arrow). In system 3, His64 located in a channel (highlighted in green) bordered by Asn62 and Trp5.

Methods

MS-EVB Methodology

A full description of the MS-EVB methodology can be found in the literature and, therefore, will only briefly be covered here.[2,73,79,80] In the MS-EVB formalism the lowest energy solution of the MS-EVB Hamiltonian matrix defines the potential energy surface upon which the system propagates, thereby dictating the evolution of the system. This matrix is expressed in terms of a dynamically adaptive basis set of empirical valence bond states (EVB),|i⟩, which form the MS-EVB Hamiltonian matrix elements from the operator where the vector represents the complete set of nuclear degrees of freedom. The diagonal elements of the EVB Hamiltonian, h, represent a specific bonding topology of the system with the corresponding potential energy as determined by the underlying classical parm99 force field for each MS-EVB state |i⟩. The off-diagonal elements, h, represent the coupling between MS-EVB states |i⟩ and |j⟩. The functional form of the off-diagonal elements is represented bywhere q is the PT coordinate, R is the distance between donor and acceptor atoms, and V is an adjustable parameter. The functional form of f(R) for the ionizable residue-hydronium pair iswhere g(q) = e( and C, α, β, γ, ε, a, b, and c are adjustable parameters. The PT reaction coordinate, q, is a geometric reaction coordinate described bywhere r = r0 – λ(R – R0) and r, λ, and R0 are additional adjustable parameters. An additional modification made to the underlying force field is the substitution of the standard harmonic approximation by a Morse potential, which allows for a more accurate description of the dissociable bond between the excess proton and the donating/accepting atom.[80] Due to the delocalized nature of the excess proton, the center of excess charge (CEC) is used to follow the charge defect. The CEC is defined aswhere c are the coefficients obtained after the diagonalization of the MS-EVB matrix at each time step and the center of charge (COC) is given bywhere the summation is over all the atoms contained in the particular MS-EVB state i with their respective partial charges q. The PT reaction coordinate ξ defines the distance between the donor molecule of interest and the CEC, as determined by ξ = | – |, where D is the donor molecule, such as the catalytic zinc-bound water. Using the MS-EVB formalism, one may generate dynamical trajectories on the potential surface obtained under varying conditions of temperature and volume/pressure. In addition, due to the high computational efficiency, it is possible to carry out accurate free energy simulations. In combination with umbrella sampling and the weighted histogram analysis method (WHAM), MS-EVB simulations have been used to study free energy profiles projected along suitably chosen variables such as ξ = | – |.

HCA II Systems

To adequately evaluate all aspects of the rate limiting PT event in HCA II, three systems (see Introduction) were simulated by means of the coupled transition path sampling (TPS) algorithm and the MS-EVB methodology (TPS-MS-EVB). The starting configurations for the three systems were obtained from previous long time MD and MS-EVB simulations as outlined in Maupin et al.,[39,62] and therefore will only be briefly described here. The HCA II structure was obtained from X-ray crystal data (PDB 2CBA)[67] and was fully solvated in a cubic simulation box of modified TIP3P[80] water molecules (L = 60 Å), where the protein was described by the AMBER parm99[81] force field and the parameters for the catalytic zinc active site were taken from previously published work.[62] The systems were then subjected to over 2.5 ns of equilibration in classical MD simulations followed by 0.5 ns of equilibration in MS-EVB. The resulting structures were subsequently simulated in biased MS-EVB simulations with the position of the CEC restricted to regions of space defining the movement of the CEC through the rate limiting PT event. Each of the restrained windows was simulated for ∼1 ns and resulted in the PMF describing the rate limiting PT event for the three systems outlined above. It was from these biased simulations that the coordinates were obtained for the TPS-MS-EVB simulations presented in this manuscript.

Transition Path Sampling

The theoretical background and numerical implementation of the TPS algorithm and its variants have been extensively documented in the literature.[74−78] In the present work, we have developed and implemented an interface between the existing methodologies of MS-EVB and TPS in order to explore the dynamics of an excess proton traversing the transition state for the rate limiting PT event. This is a multistep procedure requiring a suitable choice for the order parameter that is capable of representing the transition from a reactant to a product state and vice versa. Although a priori knowledge of the underlying free energy profile is not required, in the present case, a detailed profile of the PMF (see section and ref (39)) was available for the system under consideration, providing a sound basis for introducing the definition of the reactant and product states. As mentioned earlier, the utilized transition path sampling method harvested only forward flux for system 1, and the full dynamical trajectory, consisting of both forward (blue arrow, system 2; pink arrow, system 3 in Figure ) and backward (red arrow) propagation segments from the shooting point, for systems 2 and 3. These details are summarized in Table .

Table 1

Systems, Order Parameters (OP), Location of the Free Energy Maxima (Fmax), and Barrier Regions Employed in the Transition Path Sampling (TPS) Simulations

						Definition of states
System	Description	Method of simulation	OP	Fmax	Definition of the barrier region for shooting of trajectories	Reactant (R)	Product (P)
1	His64 in its inward orientation and included in the MS-EVB basis set	Forward flux TPS	rZn-CEC	rZn-CEC = 5.2 Å	5.1 Å ≤ rZn-CEC≤ 5.3 Å	CEC already at the barrier region	CEC taken up by inward conformer of His64 side chain
2	His64 in its outward orientation and excluded from the MS-EVB basis set	Full dynamical TPS	rZn-CEC	rZn-CEC = 6.5 Å	6.4 Å ≤ rZn-CEC≤ 6.6 Å	Entry of CEC from the protein surface in to the active site at the barrier region	CEC taken up by zinc-bound hydroxide ion
3	His64 in its outward orientation and included the MS-EVB basis set	Full dynamical TPS	rZn-CEC	rZn-CEC = 6.5 Å	6.4 Å ≤ rZn-CEC≤ 6.6 Å	Entry of CEC from the protein surface in to the active site at the barrier region	CEC taken up by outward conformer of His64 side chain

Order parameter and definition of reactant/product states

Forward flux transition path sampling

For system 1, where His64 is present in its inward orientation, the deprotonation of the zinc-bound water molecule and subsequent protonation of the His64 side chain requires the CEC to pass through the free energy barrier at rZn-CEC = 5.2 Å[39] located around W3a. The observed PT pathways across this barrier region fall under the category of a fast proton pathway where the CEC traverses the hydrogen bonded pathway containing W1–W2–W3a, where a major stabilizing contribution comes from the hydrogen bonding between the W3a and Tyr7 side chains. In the context of this order parameter, the barrier region was defined as a narrow range from 5.1 ≤ rZn-CEC ≤ 5.3 Å. Accordingly, the product region state P was defined as corresponding to rZn-CEC > 5.3 Å. All the trajectories harvested in this system were generated from the barrier region or from the edge of the reactant state R (rZn-CEC < 5.1 Å). In addition to the rZN-CEC based order parameter, an additional order parameter defined as rNδ1-CEC was chosen to monitor the fate of the forward flux system. A dynamical trajectory originating from the shooting zone is designated as reactive if rNδ1-CEC (distance between Nδ1 of His64 and CEC) falls below 1.0 Å on forward propagation. This order parameter was selected for system 1 due to the relatively small distance separating states R and P, and the rapid movement of the CEC between these states.

Full Dynamical Transition Path Sampling

The distance, rZn-CEC, between the zinc atom and the CEC has been chosen as the order parameter in setting up the transition path sampling simulation for systems 2 and 3. It may be noted that rZn-CEC by definition is a collective variable, dependent on the instantaneous location of the CEC, which in turn depends on the position of the excess proton delocalized across several solvation shells.[2] A close inspection of the free energy profile of deprotonation of the zinc-bound water molecule[39] reveals that the free energy maximum appears at rZn-CEC = 6.5 Å with a slow decrease of the free energy for larger rZn-CEC distances. We define the barrier region as a narrow range centered at the barrier top defined by 6.4 ≤ rZn-CEC ≤ 6.6 Å. Although this choice is arbitrary, it is expected to be reasonable when probing the dynamics of the system with high free energy. In addition to the barrier region, the reactant state R was defined as corresponding to rZn-CEC > 6.6 Å, and the product state P was defined as corresponding to rZn-CEC < 6.4 Å. As shown in Table , the state P corresponds to the uptake of CEC by the zinc-bound hydroxide in system 2 and by the His64 side chain in its outward orientation for system 3. In both cases, the state P encompasses not only the excess proton being taken up, but also the intermediate Zundel-like state where the CEC is primarily centered on the W1 water site (Figures and 2).[39] State R, as defined, includes the state points representing diffusion of the CEC along the protein surface at larger distance, and its entry into the active site cavity to reach the barrier region.[39] Therefore, a full dynamical transition path within systems 2 and 3 includes segments corresponding to entry of the CEC into the active site, fluctuations in the barrier region, and eventual passage to the respective product states.

Generation of the initial reactive trajectory

A reactive trajectory is generally defined as one that connects the reactant region (defined as state R) and the product region (defined as state P) after passing through the barrier region one or several times. In order to conduct transition path sampling, at least one trajectory that traverses from state R to state P or vice versa is required as a seed for the Monte Carlo sampling of trajectory space.

While generating initial trajectories for system 1, the dynamical migration of W3a away from Tyr7 was found to push the CEC toward W3b, thereby resulting in a relative slowing down of His64 protonation. Preliminary searches for initial transition paths within this model have, however, been focused on those originating from the region between Gln92 and Asn62. These transition paths were generated by restricting the CEC to reside within 5.2 Å of Zn2+ using a biasing harmonic potential with a force constant of 40 kcal mol–1 Å–2. In addition, we have also used, in a few simulation runs for system 3, the same harmonic potential to maintain a mean distance of 11–12 Å between Asn62 and Trp5. The latter ensures that the exit channel, which is bordered by Asn62 and Trp5, remains in an open state. The open state of the channel was found to facilitate the searching of the MS-EVB state space as the CEC moves toward the His64 side chain and for the protonated His64 (H+His64), carrying the excess charge away from the active site toward the exit channel. Under such conditions, the excess proton falls on the inward His64 side chain rapidly within 500 fs. These biased trajectories were then subjected to a further unbiased propagation for ∼1 ps. The shooting points for system 1 were chosen to lie in a narrow barrier region spanning 4.8 ≤ rZn-CEC ≤ 5.4 Å with a reactive trajectory satisfying the criteria that the rNδ1-CEC distance falls below 1.0 Å on forward propagation. The results of this process resulted in the generation of four appropriate seed trajectories.

Full dynamical transition path sampling

The search for a suitable seed trajectory was conducted by collecting a series of configurations from biased samplings of the region centered between 6.2 < rZn-CEC < 6.8 Å. Each of the selected configurations in this region was then submitted to an additional MS-EVB simulation[73] in the presence of a biasing umbrella potential for 1.5 ps in the constant volume and constant temperature (NVT) ensemble with a time step of 0.5 fs. After the initial 1.5 ps, the bias was removed and the system was allowed to equilibrate for another 10 ps, with all other conditions remaining unaltered. The use of this two-step procedure ensures enhanced sampling of the region near the barrier top, thereby facilitating the search for an initial trajectory. In the event that a trajectory developed complications, due to large fluctuations in the system energy upon removal of the biasing potential, an intermediate biasing trajectory was introduced before the unbiased simulations. The intermediate biasing simulations utilized a moderate force constant of 40 kcal mol–1 Å–2 and were conducted for 1.5 ps. Systems utilizing this intermediate biasing protocol were then checked for stability after at least 1 ps of simulation in the unbiased simulation. Stable unbiased trajectories were then simulated for 10 ps.

Sampling of the Transition Paths

Forward Flux Transition Path Sampling

Within this modified framework, we have harvested a transition path ensemble starting with a set of 4 seed trajectories each of length 1 ps, saved at an interval of 0.5 fs. Subjecting each of these seed trajectories to 50 shooting trials followed by forward propagation for 1 ps is found to yield a total of 130 reactive trajectories. The same methodological details of sampling were followed in the forward direction as those mentioned earlier. The high acceptance ratio and short length of each of these reactive trajectories (1 ps) indeed show no or little recrossing into the barrier region and a more facile PT to the His64 side chain in comparison to the baseline model discussed earlier.

Full Dynamical Transition Path Sampling

Using each of the initial reactive trajectories, a conventional TPS algorithm[75] was utilized to carry out a Monte Carlo sampling of the trajectory space. For this purpose, a microcanonical ensemble was generated for transition paths using forward and backward shooting moves. First, a time slice of a trajectory was randomly chosen from the barrier region of the seed trajectory and the details of the coordinates, velocities, and forces at that point were retrieved. A small random displacement from a zero mean Gaussian distribution was then added to the original momentum, generating a new momentum, which was then rescaled in order to conserve the total linear and angular momenta as well as the total energy of the system.[75] For an acceptable move to this new phase space point, the trajectory was propagated in both forward and backward directions using the TPS-MS-EVB simulation method.[73] If the total trajectory satisfied the criteria of being a reactive trajectory, it was saved and used as the starting point for the subsequent Monte Carlo move. However, if the new trajectory was a nonreactive trajectory, a new point was chosen from the barrier region and the steps outlined above repeated until another reactive trajectory was generated. For each trajectory that was used in the current seed, a maximum of 50 trial shootings were carried out in the TPS-MS-EVB algorithm with an acceptance ratio of ∼18%. Following this process, a transition path ensemble was generated consisting of 50 reactive trajectories of length 10 ps each, saved at an interval of 50 fs.

Occupancy analysis of the mechanism of PT

In order to probe the dominant mechanism of PT while the system resides at the free energy barrier top, the active site water molecules that are stable throughout the reactive trajectory were investigated. For this purpose, occupancy plots for the water oxygen atoms inside the barrier region, region P, and region R were created using the volmap tool in VMD.[82] In addition, the occupancy plots for amino acid residues Asn62, Asn67, Trp5, and Tyr7 were also created due to their close proximity to the stable water molecule locations inside the active site cavity.

It is well documented that a judicious choice of collective variable(s) to represent the order parameter is crucial for the construction of an accurate reaction coordinate.[60,78,83] This article presents several such order parameters. To understand the coupling between PT and protein dynamics, we focus on the side chain dihedral angles of active site residues already found to be kinetically important from experiments. Alternatively, several equivalent distance based order parameters could be chosen. However, the use of side chain dihedral angles is found to be better in discriminating the fast PT mechanism from the slow ones. The relative importance of these additional variables to the dynamics along an accurate reaction coordinate should ideally be assessed in terms of the commitment probability of the transition paths. Obtaining the distribution of commitment probability from the barrier top (or, more precisely, the separatrix) is computationally demanding in spite of recent developments to improve the efficiency of such a task.[40,78,84−88] Alternative machine learning based methods within transition path sampling may be used in the future to extend the current observations and arrive at an accurate reaction coordinate.[60]

Results and Discussion

System 1: Fastest PT Paths Mediated by His64 in Its Inward Orientation

From the four seed trajectories and subsequent 200 shooting trials, an ensemble of 130 reactive trajectories were generated that propagated the system from the free energy barrier region around rZn-CEC = 5.2 Å to protonate the His64 side chain. Analysis of these trajectories shows that, by far, the most efficient and fastest PT event occurred when His64 was in its inward orientation. Temporal variations of the reactivity order parameter, d2≡ rNδ1-CEC, along selected seed trajectories as well as some of the transition paths harvested are depicted in Figure . The generated reactive trajectories in system 1 could be further classified into two categories based on the values of τrxn:

Figure 3

Variation of the order parameter, d2 = rNδ1-CEC, along four initial dynamical trajectories (top), and four harvested transition paths when His64 is in its inward orientation and participating in the PT (bottom).

Ultrafast PT with τrxn ≤ 50 fs.

Fast PT with 51 ≤ τrxn ≤ 250 fs.

Out of the 130 transition paths harvested, 52 were found to be ultrafast, while 76 were classified as being fast. Only two were found to have τrxn equal to 400 and 551 fs, and were not considered for further statistical analyses in system 1. The associated distributions of 7 side chain dihedral angles (χ, i = 1,7) for residues in close proximity to the PT event were evaluated (Table ). The systems presented here were monitored along the entire 1 ps length of each transition path. In addition, the distributions were examined along the truncated segment when the His64 side chain remains protonated. Differences, if any, between the distributions of χ along the full and truncated transition paths are expected to reflect coupling of a preferred orientation of the associated side chain to the PT. Our results indicate coupling of χ1, χ2, and χ3 to PT in system 1. All other dihedral angles had nearly superimposable distributions along the full and truncated trajectories (Supporting Information) and are not discussed for system 1.

Table 2

Side Chain Dihedral Angles Investigated along the Fastest PT Paths Starting with His64 in Its Inward Orientation

Side chain dihedral angle	Residue	Constituent atoms
χ1	His64	N-Cα-Cβ-Cγ
χ2	His64	Cα-Cβ-Cγ-Nδ1
χ3	Asn62	Cα-Cβ-Cγ-Nδ2
χ4	Asn67	Cα-Cβ-Cγ-Nδ2
χ5	Trp5	Cβ-Cγ-Cδ1-Nε1
χ6	Tyr7	C*-Cα-Cβ-Cγ
χ7	Gln92	Cβ-Cγ-Cδ-Nε2

The side chain of His64 was found to remain mostly directed toward the active site cavity in a position to facilitate fast proton uptake. As shown in Figure , along the ultrafast PT paths, the His64 side chain adopted χ1 values resembling both in (positive χ1) and out (negative χ1) conformations. The ultrafast transfer appeared to be facilitated by those χ2 values that correspond to the imidazole ring remaining nearly parallel to the catalytic zinc.[62] Along fast PT paths (Figure ), the His64 side chain started from a wider range of conformations. But, upon PT, the distribution peaked between 90°–180° for χ1 (typical of in orientation). The distribution of χ2 shifted to higher positive values upon PT. The latter corresponds to ring flipping of the imidazole ring from parallel to a more perpendicular orientation relative to the zinc.[62] Such transitions along the fast PT paths would make τrxn marginally longer compared to the ultrafast PT path.

Figure 4

Distribution of key side chain dihedral angles, χ1, χ2(His64, upper panel), and χ3(Asn-62, lower panel) in system 1 for ultrafast proton uptake by His64 in its inward orientation with τrxn ≤ 50 fs. Each distribution is sampled along full dynamical trajectories (black) and along the truncated segment when the excess proton resides on the His64 side chain (red). The black and red lines for χ2 are depicted by dotted lines for clarity. Asn-62 shows a clear conformational preference when His64 is protonated (red, lower panel).

Figure 5

Distribution of key side chain dihedral angles, χ1, χ2 (His64, upper panel), and χ3 (Asn62, lower panel) in system 1. Fast proton uptake by His64 is in its inward orientation for transition paths with 50 < τrxn ≤ 250 (fs) (black). Each distribution is sampled along full dynamical trajectories (black) and along the truncated segment when the excess proton resides on the His64 side chain (red). The black and red lines for χ2 are depicted by dotted lines for clarity.

Apart from χ1 and χ2, the side chain dihedral angle of Asn62, χ3, also showed a distinct change in its distribution upon PT (Figures and 5). Along both ultrafast and fast PT paths, the Asn62 side chain oriented away from the active site (χ3 below −90°) side chain when the excess proton is attached to the His64 side chain in its in orientation. Such reorientation was encountered in spite of starting with a wide range of values of χ3 (Figure ) and appeared to trigger the eventual transfer of the excess proton from the barrier region to the His64 side chain. The correlation between the speed of the PT event and χ3 indicates that Asn62 pointing away from the active site creates an environment conducive to the fast protonation of His64 (i.e., within 250 fs). Though statistically insignificant within system 1, the two trajectories with τrxn = 400 and 551 fs also exhibit fast proton uptake upon localization of the Asn62 side chain away from the active site.

A preliminary analysis of dynamical changes of the order parameters listed in Table has also been carried out by focusing on a narrow (∼10 fs) interval around τrxn along both ultrafast and fast trajectories (Supporting Information). These data suggest, in addition to Asn62, transient but sizable reorganizations of χ4-χ7 may also trigger the PT step.

The observation of ultrafast and fast PT events clearly reveals the highly optimized nature of HCA II in moving a proton between the active site zinc water/hydroxide and His64 via W1–W2–W3 with the His64 side chain in the inward orientation and the Asn62 side chain pointing away from the active site. A minimum recrossing of the barrier region, coupled to the underlying PMF for PT, reveals that system 1 should possess a significantly large rate coefficient due to favorable activation energy and a relatively larger pre-exponential factor for the rate coefficient. As shown later, this rate is expected to be larger as compared to systems 2 and 3.

System 2: PT Paths Excluding His64

To evaluate the specific role of His64 on the rate and efficiency of the rate limiting PT event, it is beneficial to evaluate the control system consisting of the HCA II enzyme with His64 in the outward conformation and not participating in the proton relay (i.e., excluded from the MS-EVB state space). This hypothetical system allows for the direct evaluation of His64’s impact on stabilizing the water cluster inside the active site and modulating the PT barrier. As discussed below, the pathways excluding His64 have been classified in five groups. For each of these groups, distributions of active site water molecules and amino acid side chains are examined and their effects on the formation of PT paths analyzed. The observations reported in this section are derived by sampling all the trajectories within a given group and hence provide a unique, statistically meaningful perspective of the enzyme function in the absence of its primary reaction pathway.

Variation of the Order Parameter along Initial Trajectories and Sampled Transition Paths

Following a combination of biased and unbiased MS-EVB simulation runs as described in the previous section, four initial trajectories were obtained and used as seeds in the subsequent Monte Carlo sampling of the trajectory space connecting states P and R. The temporal variation of the order parameter, rZn-CEC, along these four initial trajectories is presented in Figure (top), which depicts various transitions that originated in state R and that cross and recross the free energy barrier region (6.4 Å ≤ rZn-CEC ≥ 6.6 Å, Table ) and ultimately end in state P. It is evident that the seed trajectories (represented by different colors) differ from each other in terms of both the residence time, τ, as well as the number of recrossings, N, into the barrier region. The variation of rZn-CEC along the transition paths in the ensemble are also highlighted in Figure (bottom) in terms of four representative trajectories with different values of τ and N. The current procedure, as seen by the various trajectories, is capable of sampling different trajectory paths such as early/late entry/exit to/from the barrier region, leading to the changes in τ and N. The sampled trajectories recrossed into the barrier region on average 20 ± 14 times before exiting the transition region for the final time.

Figure 6

Variation of order parameter during (top) seed trajectories and (bottom) representative paths from the transition path ensemble. The region rZn-CEC < 6.4 Å is P, and rZn-CEC > 6.6 Å is R, while the narrow region between 6.4 and 6.6 Å represents the top of the free energy barrier. This figure highlights the diversity of the trajectories used in the present study in the ensemble of transition paths.

Classification of fast and slow transition paths

From the observed wide distribution of values for both τ and N, the transition path ensemble was classified into five broad categories (Table ). This classification was primarily based on increasing values of τ. As expected, more time spent at the barrier region resulted in a higher number of recrossings, N. The transition paths belonging to group A are designated as fast paths, those belonging to groups B and C are defined as moderate paths, and the transition paths belonging to groups D and E are labeled as slow paths. The ensemble of transition paths generated for system 2 consisted of predominantly slower paths with the relative populations of the fast, moderate, and slow categories mentioned above being approximately 14, 34, and 52%, respectively.

Table 3

Classification of the Transition Paths in System 2 in Terms of the Residence Time, τb, and the Number of Recrossing Events, Nb, at the Free Energy Barrier (i.e. transition state)

Group	Residence time,τb (ps)	Number of recrossings, Nb	Path type
A	τb ≤ 0.5	2–5	Fast
B	0.5 < τb ≤ 2	3–10	Moderate
C	2 < τb ≤ 5	5–30	Moderate
D	5 < τb ≤ 8	10–40	Slow
E	8 < τb ≤ 10	20–65	Slow

In the following discussion, the mechanistic relevance of such a classification is demonstrated by highlighting the factors that distinguish the fast paths from the slower ones. For a given set of paths, we first examined the entire length of each transition path carrying the excess proton from region R through the barrier top and terminating in region P. In addition, we separately investigated the barrier states in a given group, defined in the present study as the time slices of each transition path in that group residing in the barrier region (6.4 Å ≤ rZn-CEC ≤ 6.6 Å, Table ).

PT along Fast Paths

The fast paths, schematically shown in Figure , are characterized by τb ≤ 0.5 ps and an average value of ⟨τb⟩ = 0.26 ps. The fast paths were found to involve a short, hydrogen bonded pathway formed by 3 water molecules (W1′-W2′-W3c) that were pre-existing in nearly optimum orientations at the time of entry into the barrier region. The water molecules nearest to the zinc-bound hydroxide, W1′ and W2′, along the fast paths mimicked the role played by wild type W1 and W2. W1′ remained hydrogen bonded to the zinc-bound hydroxide and Oγ1 of Thr200 while W2′ was stabilized by hydrogen bonding to W1′. When the CEC resided in the barrier region, a third water molecule, labeled in this study as W3c, was found to form the fast path by sharing the excess proton with W2′. No analogue of W3c was observed in the proton pathways detected in wild type HCA II. Only a few barrier states were encountered along these fast paths. Interestingly, the side chain of Gln92 (highlighted as a blue van der Waals surface in Figure A) was found to partition the active site roughly into two regions. The fast proton paths passed along the left side of Gln92, along the active site boundary lined predominantly by hydrophobic residues such as Phe93, Phe95, and Leu118. Other active site water molecules were located mostly on the opposite side of the Gln92 side chain, close to the active site wall bordered by hydrophilic residues such as Asn62 and Asn67 (highlighted by a pink surface in Figure A). A couple of water molecules were detected around the crystallographic location W2. However, they were not part of the proton relay. No localized water molecule was found at the site of W3b (i.e., active site water that is hydrogen bonded to both Asn62 and Asn67). The water molecule nearest to the spatial position typically occupied by W3b was found to contain an equal distribution of waters in two distinct locations, both of which were within hydrogen bonding distance with Asn62. Examination of the water structure around Tyr7 for all the barrier states revealed two water molecules that remained hydrogen bonded to the hydroxyl side chain of Tyr7. However, none of them connect, even transiently, to W2′. Therefore, PT along the fast paths is concluded to contain only W3c-W2′-W1′. No analogue of this path was observed in system 1 when His64 is in the inward conformation.

Figure 7

Fast PT pathway. (A) Gln92 (shown in blue) divides the active site waters into two chains. Fast PT occurs through the left chain along the active site wall with hydrophobic residues (not shown) and opposite to the wall with hydrophilic residues such as Asn62 and Asn67 (shown in pink). (B) The three waters that mediated the fast PT path, W1′, W2′, and W3c, lead to the zinc-bound hydroxide (licorice). The CEC is represented by a yellow sphere while the gray sphere indicates the van der Waals surface of the zinc ion.

Water and Proton Path Distribution along Slower Paths

In contrast to the fast paths, the slow paths (Figure ) typically entered the free energy barrier region near the residues Asn62 and Asn67. Inspection of these paths revealed that the side chain of Gln92 hindered the motion of the CEC toward the region where the fast path occurred (i.e., W3c-W2′-W1′). Under these circumstances, the CEC is observed to explore the water clusters extending in the region bordered by Asn62 and Asn67, Trp5 and Tyr7. As shown in Figure , there are several water molecules in this region that can potentially form a hydrogen bonded water cluster to transfer the proton to region P. To gain insight into the dynamical nature of the water cluster and the average location of the water oxygen atoms, the distributions of water molecules within a distance of 6.5 Å from the catalytic Zn were averaged while the CEC resided in the barrier region for groups A, C, and E (Supporting Information Figure S3). While the locations of a water molecule analogous to W1 were not difficult to identify in the water density plots, some disordering of its position was clearly visible as the values of τ increased. In contrast, the hydration site corresponding to W2 was also found to be more disordered for the barrier states of slower paths. These observations were further substantiated by comparing the barrier state distributions to those of long-lived water sites in the active site cavity, sampled from the entire lengths of fast and slow transition paths. The active site water isosurfaces (Figure ) along fast paths (Group A) revealed the dominant location of water molecules, facilitating the persistence of an isolated and stable three-water cluster as discussed above (Figure ). However, for the slower paths of groups C and E, the respective isosurfaces gradually moved toward the Asn62 side chain, away from the location preferred by the fast paths. This dislocation of higher populated hydration sites from one side of the active site cavity to another for slower paths is highlighted in Figure B and C. Consequently, the moderately slow paths were detected in the region between W2′ site of the fast path (Figure ) and the Asn62 side chain. Along the slow paths (e.g., Group E), most probable hydration locations (Figure C) are found closest to Asn62 side chain. Accordingly, compared to the fast paths, the slower PT pathway appeared as a highly branched network of water molecules that was shifted to the opposite side of the active site cavity. Finally, as shown in Figure D, the excess charge was also stabilized closer to the Asn62 side chain with increasing residence time, τ. Consequently, the dynamics along moderately slow and slow transition paths are also expected to differ substantially from that presented in connection with the fast paths. This difference is examined in detail next.

Figure 8

Snapshots of a slow transition path entering the barrier region near Asn67 (A) and a typical barrier state with the CEC (yellow) residing near Asn62 (B).

Figure 9

Isosurfaces representing active site water molecules with 90% occupancy. (A) Isosurfaces along transition paths (red) and barrier states (yellow) for group A; (B) Isosurfaces along transition paths for groups A (red, transparent), C (orange), and E (green); (C) Isosurfaces for the barrier states of groups A (red), C (orange), and E (green). (D) Isosurfaces showing sites with 90% occupancy of the CEC when in the barrier region for groups A (red), C (orange), and E (green). The transparent surface around the amino acid residues corresponds to the van der Waals surface of Gln92. The zinc ion and the hydroxide bound to it are shown using their corresponding van der Waals surfaces.

For the moderately slow transition paths (trajectories in group C) with ⟨τ⟩ = 3.9 ps, a pair of water molecules (shown as green molecules in Figure ) was located in 28% of all the barrier states and treated as the W1–W2 pair, conducive to the formation of a path resembling the fast path structure. At these barrier states, a water molecule disconnected from this pair (shown in purple) was found to be hydrogen bonded to the Tyr7 side chain in addition to an analogue of W3a. However, less than 10% of the barrier states were found to be associated with such configurations. The barrier states were found to be dominated by a different pair of water molecules, shown in purple in Figures B, C, and D, at the locations of W1 and W2, with the W1 site now populated by the water molecule, mentioned above, that was hydrogen bonded to the Tyr7 side chain (shown in purple, Figure A). The water pair of W1 and W2 (shown in green) was now found in the region between Gln92 and Asn62. Interestingly, for either of the two pairs (green or purple), irrespective of their location at the active site, the pairing water molecules always remained hydrogen bonded to each other at all barrier states. It is also noted that the side chains of Asn62 and Gln92 both contributed to the variation in the relative positions of the green and purple water pairs by stabilizing one or two water molecules through hydrogen bonding. A closer inspection of the CEC′s location reveals that it spent 37.5% of the time within 3.5 Å of the Asn62 side chain. During this period, no direct path into region P could be detected. Finding its path disrupted, the CEC continued to fluctuate across a disordered analogue of W2 and two other water molecules located in the intervening space between Asn62 and W1, thus leading to a much longer residence time at the barrier region. The CEC devoted less than 2% of its time exploring the region that brought it within a distance of 4 Å of Gln92. Within this short excursion near Gln92, no water cluster resembling the fast path and contributing to the solvation shells of the CEC were encountered. These observations are supportive of the research that indicates small nonbranched water clusters are more conducive to PT than larger, branched water clusters.[43]

Figure 10

Transient hydrogen bonded clusters formed in the active site obtained by sampling of the barrier states of moderately slow and slow transition paths. Panels A, B, and C correspond to barrier states of group C and highlight the three-water mediated path in addition to the stabilization of the branched water network by Asn62 and Gln92. A transient barrier state formed in Group E, D shows the Asn62 side chain in its minor side chain orientation hydrogen bonded to His64 and relocation of the excess charge (yellow) closer to the putative pair of W1 and W2 (purple). The two special pairs of water molecules which populate these barrier states are shown in green and purple in each case.

An interesting scenario emerged in our examination of the mechanism of PT along the slow transition paths where the trajectories belonging to group E, characterized by ⟨τ⟩ = 9.1 ps, were analyzed. As in the case of moderately slow paths, the CEC entered the barrier region from state R near Asn62 and Asn67, close to the location of W2. Once again, a special pair of water molecules (shown in green, Figure D) was detected. These water molecules, due to their location and orientation, are believed to play the role of W1 and W2 along the path. However, unlike the moderately slow paths, these two water molecules remained hydrogen bonded in about 85% of the barrier states. A third water molecule (Figure D) did establish a hydrogen bond to the W1 site, facilitated by its stabilization in that region by Gln92. In 15% of the barrier states, the water molecule previously belonging to the green pair and occupying the W2 site was found to be hydrogen bonded to the side chain of Asn62. Following the dynamical movement of the CEC along the entire length of the slow reactive trajectories, it was found to fluctuate in the same region as observed in the case of moderately slow trajectories with repeated recrossings into the barrier region. During such nonreactive fluctuations, the CEC spent 82.1% of the time within 3.5 Å of the side chain Nδ2 atom of Asn62. The CEC also spent 6.2% of its time at the region within 4 Å of Gln92. However, less than 1% of the states explored during this excursion corresponded to the barrier region. Consequently, the loss of connectivity between W1 and W2 coupled to the predominant localization of the CEC near Asn62 made the transition to state P much slower in this group compared to the instances considered earlier. It is worth noting that the structure and dynamics observed in the slower paths is most likely the dominant pathway in the His64Ala mutant, and this helps explain the presence of enzyme activity, although it is greatly reduced (∼20-fold slower) as compared to the WT enzyme. In addition, evaluating the fast and slow pathways adds insight into the chemical rescue phenomena involving His64Ala and chemical rescue agents such as imidazole. In recent work[48] it has been hypothesized that the chemical rescue agent needs to bind deep in the active site, thereby setting up small nonbranched water cluster pathways, instead of binding on the rim of the active site, which would favor larger branched water cluster pathways.

In view of the differences in average and dynamical occupancies of water molecules along the transition paths, the microscopic nature of the CEC is expected to assume significantly different attributes as we go from fast to slow paths. As shown in Figure , the slow paths are dominated by the formation of Eigen-like cations, c2 ≥ 0.65, while the fast paths possess substantial contributions from the Zundel-like cations, c2 ≈ 0.5 as well. Therefore, it is concluded that the transitions between regions

Figure 11

Distribution of the square of the maximum amplitude of the CEC, c2, showing the variation of the microscopic nature of the excess charge along trajectories belonging to (A) group A, (B) group C, (C) Group D, and (D) Group E with increasing residence time, τ, at the barrier region.

Role of Active Site Residues

We next probe the role of some selected amino acid residues along the reactive trajectories and at the barrier states for groups A, C, and E. It is important to note that residues such as Asn62 and Gln92 bordering the active site cavity have been repeatedly implicated in the preceding discussion on the variation of the underlying mechanism of fast and slow paths. In the present work, the distributions of side chain orientations for Gln92, Asn67, Asn62, Trp5, and Tyr7 at the barrier states of groups A, C, and E are evaluated. For the sake of comparison, we have also considered the analogous distributions obtained by averaging over the complete transition paths in each group. As the residence times at the barrier states increased, the distributions for the complete slower paths were found to resemble the behavior at the barrier states more than those associated with the faster paths. The most important changes were observed for Asn62 and Gln92 side chains in going from the fast to the slow paths. These are shown in Figure . The preferential orientations of these side chains clearly demonstrate how they contributed to the localization of the CEC between Gln92 and Asn62 through hydrogen bonding with the intervening water clusters. Along the slow paths, the Asn62 side chain existed predominantly in the conformation pointing toward the catalytic zinc. However, during its transient excursions to the minor conformation away from the active site, it formed a hydrogen bond with the His64 side chain, as shown in Figure D. It was with these minor side chain conformations of Asn62 that the CEC moved toward the center of the active site and a rare PT path was formed, such as the one shown in Figure D. Interestingly, the side chain of Gln92 exhibited side chain flexibility that was in direct contrast to the behavior of Asn62 as we go from fast to slow paths. Both side chain dihedral angles of Gln92 were found to span a wider range of values at the barrier states belonging to groups C and E, while they are relatively more localized at the barrier states of group A. Both Trp5 and Tyr7 exhibited transitions between two dominant side chain orientations at the barrier states. However, the persistence of both the conformations at the barrier states with marginal changes in the relative population makes it difficult to infer their exact dynamical role in decoupling the time scales of the transition paths, unlike what was found for Asn62 and Gln92. Although the mechanistic details reported here are arguably dependent on the time span probed at the free energy barrier, extending a set of ten slow paths up to another 10 ps shows that the dual control exerted on the CEC by Asn62 and substantial branching of the active site water network near Gln92 still remained pertinent.

Figure 12

Distribution of side chain dihedral angles of amino acid residues Asn62 (a, b) and Gln92 (c, d) at the active site sampled from barrier states (a, c) and from full transition paths (b, d). In each panel, distributions corresponding to different transition path groups have been shown (red: Group A, fast; green: Group C, moderate; and blue: Group E, slow).

PT Paths Mediated by His64 in Its out Orientation

In section the discussion focused on important dynamical changes associated with the movement of the CEC from the free energy barrier region toward the active site (i.e. toward the reactant, region R). In this section the focus is on the probable dynamical paths that carry the excess charge away from the barrier region into region P. Analysis of the present transition path ensemble clearly showed that in the absence of a strong proton sink, such as His64 in its inward conformation, the CEC utilized the water network extending through the intervening region between Asn62 and Asn67 to exit from the active site, which is hypothesized to be the primary pathway in the His64Ala mutant. We have next investigated the model where the EVB state search included His64 in its outward orientation. Once again, the free energy barrier appears at rZn-CEC = 6.6 Å around the location of W3b. The same protein motions come into play as discussed in section and will not be repeated here. In our sampling of the transition path ensemble, we did not observe the formation of a water network originating at this barrier region and leading to the outward conformer of His64.

The absence of any proton exit path leading to His64 in its outward conformation, even if it is included in the EVB basis sets, warrants further investigations. The distance of closest approach, r, between the side chains of Asn62 and Trp5 is distributed roughly between 3–7 Å, indicating a rather narrow opening of the channel that the excess proton (e.g., Eigen/Zundel cation, His64H+) must transverse to enter/leave the active site cavity. Investigation into the transfer of the excess proton from the free energy barrier top toward His64 in an alternative model where the His64 side chain occupied the out conformation, and directly participated in the PT path due to its inclusion in the MS-EVB basis set has previously been explored.[39] Continuing this investigation, an extensive search for a suitable dynamical trajectory that delivers the excess charge to the side chain Nδ1 atom of His64 using a combination of biased and unbiased MS-EVB simulation runs of lengths 10 ps −1 ns has been conducted. Contrary to expectations, a productive PT path terminating on His64 was not encountered; instead, all stable trajectories allowed the CEC to exit from the active site through the water network extending in the region between Asn62 and Asn67, as mentioned above. In the case where the CEC was restricted to reside within 3 Å of the Nδ1 atom of the His64 side chain and the distance, r, was fixed at about 11 Å, a fast PT was observed in less than 20 fs and has been represented in Figure . These simulation studies indicate the importance of relaxation of the channel housing His64 in its outward conformation on the formation of a PT path involving His64.

Figure 13

Relaxation of the mouth of the channel bordered by Asn62 and Trp5, allowing placement of the excess charge (yellow) within 3 Å of the His64 side chain (A) and the following rapid transfer of the proton to His64 (B). The shortest distance between the side chains of Asn62 and Trp5 at the mouth of the channel is found to be more than 11 Å for any such transfer to be observed.

The results detailed above are a crucial outcome of the present study. The residence time and the number of recrossing events at the free energy barrier region are accordingly expected to be governed by the time scale of relaxation of this channel coupled to a concurrent excursion of the CEC into the barrier region. A close inspection of the associated PMF curve for this model (Figure 5 in ref (39)) reveals a transition between a Zundel-like and an Eigen-like cation as the system exits from the barrier region, an approximately thermoneutral reaction, and an overall barrier to PT that is roughly 1 kcal/mol higher than when His64 occupies the inward conformation. The detailed results obtained from the presented transition path sampling simulations indicate that the PT to His64 in the outward conformation will most likely have a significantly smaller pre-exponential factor for the overall PT rate coefficient as compared to the inward conformation due to the large number of recrossing events and the entropic cost associated with the reorganization/relaxation of the water cluster and surrounding residue configurations. The elevated barrier and reduced pre-exponential factor for His64 in the outward conformation as compared to the inward conformation highlights the preference of HCA II to utilize His64 in the inward conformation over the outward conformation, while also revealing how protein dynamics can modulate the rate of a reaction by impacting factors that are expressed in the pre-exponential factor of the rate coefficient (i.e., recrossing and entropic contributions).

Conclusions

We have presented a detailed analysis of how protein dynamics (specifically the motions of nonreactive residues in the active site and coupled solvent dynamics) influence the rate limiting PT event in HCA II. Within the framework of MS-EVB, the dynamics of the excess proton is followed in terms of the center of excess charge (CEC), thereby probing dynamical reorganization of active site water molecules that are expected to contribute to the fluctuating solvation shells of the CEC. In addition, our present study sampled the dynamics at configurations populating a narrow window on the top of the highest free energy barrier encountered by the CEC during its passage into/out of the active site of the enzyme. Typically the lengths of transition paths in an ensemble range between 0.1 and 1 ps depending on the system being investigated.[40,84,85,89−92] In view of these systems, the present ensemble of transition paths explores the top of the barrier region for a very long duration, i.e. 0.5 ns, thus ensuring sampling of alternative mechanisms. Delineation of the transfer mechanism in terms of fast, moderately slow, and slow proton paths clearly demonstrates how the “reactive motion” of the CEC along the chosen order parameter, rZn-CEC, may be directly controlled by the dynamical preferences of active site residues, such as Asn62 and Gln92, which otherwise do not participate in the proton relay, but are observed experimentally to favorable impact the kcat/KM.

The major outcome of the present study is an understanding of how slower protein and water motions can influence PT through the HCAII enzyme active site. It is clearly demonstrated that an efficient/fast PT pathway typically involves W1, W2, and a third water molecule. However, hydrogen bonding with residues such as Asn62 and Gln92 can substantially affect the reorientation of this third water molecule (such as W3b or W3c), causing a slowing down of the overall transfer. Although the sites of W1 and W2 are populated along most of the transition paths, they are found to undergo dynamical exchange with other active site water molecules, that in turn may cause transient disruption of hydrogen bonded connectivity between W1 and W2, slowing the transfer to/from the zinc-bound hydroxide/water. Additionally, the uptake of the excess charge by His64 in its outward orientation appears to be ineffective due to the large fluctuations required at the mouth of the channel where the His64 resides. The same fluctuations appear to be important in allowing the passage of His64 with an excess charge from its inward conformation toward the solvent.

The results presented support the rare dynamical variations of the complete PT path that has been recently observed in the neutron diffraction studies.[47] A dual control is exerted by Asn62 and Gln92 in deciding if the CEC, trapped at or near the barrier region for longer than 10 ps, would eventually fall back on the zinc-bound hydroxide (provided a 2–3 water-mediated path). It is also found that hydrogen bonding of active site water molecules by Gln92 introduces substantial branching of the active site water network, thereby slowing down the PT step. When Asn62 is in its minor orientation pointing toward His64, the motion of CEC toward zinc is facilitated by the loss of hydrogen bonding between the Asn62 side chain and its nearby active site water molecule (W3b or its analogue). In a couple of cases, Asn62 remained hydrogen bonded to the solvation shell of the CEC, thereby pushing it toward the region between Asn62 and Asn67. In these cases, the CEC did not fall back on the zinc-bound hydroxide even after 22 ps.

To understand the implications of these results in the context of the overall function of HCA II, a few qualitative inferences may be drawn. In the wild type enzyme, at least three different mechanisms may be operative in the exit/entry of the excess proton from/to the active site following the deprotonation of zinc water. The least probable one would encounter the largest free energy barrier and/or smallest pre-exponential term for the rate coefficient, and allow the slowest PT, with important kinetic control being exerted through residues such as Asn62 and Gln92. Although not so crucial in the wild type enzyme, this mechanism may account for the residual activity in the mutant His64Ala, where the His64 mediated route is absent. In the wild type enzyme, the outward conformer of His64 may indeed provide a second and marginally more probable proton relay in view of the lower free energy barrier encountered along these paths. However, the kinetic efficiency of this mechanism appears to depend additionally on the coupled dynamics of the CEC, Asn62, and Trp5, which is in agreement with available experimental mutation studies. The most probable mechanism, with the lowest estimated free energy of activation and favorable pre-exponential coefficient, involves the inward orientation of His64.

The results presented here help illuminate the role of protein dynamics in enzyme catalysis, which remains one of the most highly debated concepts in enzymology.[93−95] The main controversy centers around what may be defined as functionally significant conformational fluctuations and how they couple to the chemical changes being catalyzed by the enzyme. Recent single molecule experiments using fluorescence resonance energy transfer[96,97] and extensive MD simulation studies[22,24] do indicate the importance of several catalytically important conformers. The results presented here further emphasize a possible need to reanalyze the way kinetics of enzyme activity is conventionally represented,[97] often in terms of only one or two variables.