Electron transfer activation of a second water channel for proton transport in [FeFe]-hydrogenase

There has long been interest in understanding the family of biological redox enzymes known as hydrogenases.^1–6 One of the main reasons for this interest is the potential biomimetic functionality the enzymes present through their ability to reversibly catalyze the reduction of protons and the oxidation of molecular hydrogen, i.e.,

for use as an energy source.⁷ The current alternatives for these processes are expensive and relatively rare. The uncatalyzed cleavage of the dihydrogen bond is very difficult, requiring on the order of +100 kcal mol⁻¹ and the best synthetic catalyst for this redox reaction is the scarce platinum metal. The hydrogenase proteins, with base metals of iron and nickel, provide a possible framework to reduce the operational cost of the production and consumption of H₂.

The hydrogenase family is categorized by its transition metal components in the enzyme active site. [NiFe]-, [FeFe]-, and [Fe]-only hydrogenases are the three main classes^2,3,8,9 and only the [Fe]-only, or cluster free hydrogenase, does not catalyze the reduction of protons or the oxidation of molecular hydrogen. [NiFe] is used typically for hydrogen oxidation, while [FeFe] participates mostly in proton reduction, although bidirectional enzymes are also known to exist. The typical functionality of each class of hydrogenase is evidenced by their different catalytic rates.

The crystal structure of the [FeFe]-hydrogenase Clostridium pasteurianum (CpI) is shown in Fig. 1(a).^10–12 The protein is composed of four subdomains, of which the active site subdomain is the largest. Within the active site domain is located the eponymous metal center of the protein surrounded by over 300 amino acids and a Fe₄S₄ cuboid.¹¹ The di-iron cluster active site and the cuboid make up what is referred to as the H-cluster (shown in Fig. 1(b)). They are joined together by a cysteine ligand bridge. The Fe atom in the di-iron cluster closest to the cuboid is termed the proximal iron and the other is termed the distal iron. In the active site, each Fe atom is coordinated by both a CN and CO ligand. The di-iron atoms are bridged together by a dithiolate ligand typically assigned to be the di(thiomethyl)amine (DTMA),¹² di(thiomethyl)ether,¹³ or propanedithiolate¹⁴ groups. Recently, it was shown through controlled metalloenzyme activation, using an active site synthetic analogue, that the DTMA group is the bridging group in the native [FeFe]-hydrogenase active site.^15,16 These results were also supported by electron paramagnetic resonance spectroscopy and Fourier transform infrared spectroscopy. Another CO ligand also bridges the two Fe atoms, but, depending on the redox state, it is also found to be bound to only the distal iron atom.¹⁷ In the remaining subdomains, three more iron-sulfur cuboids, denoted as clusters A, B, and C, are situated increasingly farther away from the active site. Finally, a Fe₂S₂ cluster is also found in the protein.

Much research has been devoted to the hydrogenase catalytic cycle in both the [FeFe] and [NiFe] enzymes.^2,4,5,17,18 The structure and spin state of the active site have been extensively studied along with the reaction mechanism of the catalysis, especially in the nickel complex.^17,19–26 Since the functionality of the hydrogenase enzymes are predicated on charge separation and/or association, many recent studies have also focused on the different pathways protons and electrons take.^27–35 In particular, in a study by Shaw and co-workers,³⁰ the proton transport pathways were explored, and it was found that the amino acids play a large role in transferring the excess proton through the protein. A recent study by McCullagh and Voth³⁴ looked at the electron transport pathway in the [FeFe]-hydrogenase as one electron hopped from the exterior of the fully oxidized protein to cluster C (electronic state C^2‑), then to cluster B (electronic state B²⁻), then to cluster A (electronic state A²⁻), and finally to the H-cluster. As the electron moved towards the enzyme active site, the pathway was shown to be thermodynamically favorable. These authors also explored the electrostatic influence of these electronic states and various coarse-grained domains of the protein in possible proton transport pathways. The present work builds upon these past studies in investigating the proton transport pathways of the [FeFe]-hydrogenase enzyme, especially via a second, electronic state-dependent water channel predicted to exist in Ref. 34, called WC2.

Nearly all of the prior studies which deal with proton transfer through the [FeFe]-hydrogenase protein into the enzyme active site focus on the residue-based pathways that we have called “WC1,” consisting of three different proton transfer pathways (shown in Fig. 2 colored in blue) originating from the protein surface near the Glu282, Lys571, and Glu368 residues, respectively.³⁰ All of these pathways include protonatable amino acids and converge near the Cys299 residue before reaching the active site. The importance of one of these residue pathways on the production and oxidation of molecular hydrogen is clear. Mutations made to Cys299, Glu279, Glu282, and Ser319, all residues that most likely participate in the pathway, were experimentally shown to block catalytic activity in the enzyme.²⁹ These results were later corroborated by a theoretical study.³⁰ Recently, a water channel (WC2) was found that was activated by electron transport to the active site and that did not involve the Cys299 residue.³⁴ Specifically, when cluster A is reduced (state A²⁻), a loop containing residues 231 to 241 was shown to separate from another loop containing residues 521 to 536. These two loops are shown in Fig. 2, colored in orange and purple, respectively. In this region, the figure also shows the presence of water molecules, which fill the space left between the separated loops near the enzyme active site. Because the production of molecular hydrogen requires the presence of two protons, we have hypothesized that this second water channel could work in a controlled fashion with the electron transfer process in bringing a second proton to the active site. This would likely occur after the initial proton was passed through the experimentally determined residue pathway. Understanding this newly predicted water channel, and the role it possibly plays in the concerted proton and electron transfer in the [FeFe]-hydrogenase may be important to fully understanding the biomimetic functionality and potential of these enzymes. This is the focus of the present paper.

In the present study, we use multistate empirical valence bond (MS-EVB) reactive MD simulations^36–40 to investigate the proton transport process in WC2 of the [FeFe]-hydrogenase. Since the charge defect associated with the excess proton prefers to be delocalized across several molecular species, in particular water molecules, the protons are known to “hop” through water wires via the intermediate Zundel-like (H₂O—H⁺—OH₂) structure. This is known as the Grotthuss mechanism.³⁸ The MS-EVB MD methodology allows for such dynamic alteration of the bonding topology during MD simulations, which can be used in general for any chemical reaction, but here is restricted to only proton transfer between water molecules. Proton transport pathways were found using MS-EVB combined with umbrella sampling, and the free energy profiles were also characterized. Mutations in state A²⁻ were also engineered in silico to better understand the water wire network and how the protein structure responds. The double mutation Gly418 to phenylalanine (G418F) and Gly502 to valine (G502V) was thus studied for state A²⁻. These two residues were chosen because of their proximity to the water wire and proximal Fe atom at the active site in WC2. Initially, the mutations were studied independently, but only when they were combined as a double mutant was the water wire disrupted.

In Sec. II, the details of the atomistic MD simulation for electronic states A²⁻, B²⁻, and C²⁻ are presented. These states refer to the reduced Fe₄S₄ cuboids moving away from the active site. We have modeled the enzyme active site as the Fe^IFe^II oxidation state. The double mutated state A²⁻ simulation setup is described as well. Also in this section the method used to find proton hopping channels and to characterize their free energy surfaces is detailed. The subsequent discussion section is organized to first look at the proton transport pathways of the electronic states and then to elaborate on their free energy profiles. The double mutant pathway and free energy profile are also explained here.

The [FeFe]-hydrogenase protein, taken from the X-ray crystal structure of Clostridium pasteurianum (PDB code 1FEH),¹¹ was used as the starting geometry for all simulations. The protein contains a total of 574 amino acid residues and six metal centers. The central metal center, the active site of the protein, consists of a [2Fe]-subcluster to which carbonyl and cyanide groups are ligated. In this work, the bridging groups were taken to be a di(thiomethyl)amine (DTMA) and a carbonyl group. The active site is covalently linked to a Fe₄S₄ cuboid metal center by a cysteine residue, while three other cuboid metal centers are present throughout the protein. The remaining metal center is of the chemical formula Fe₂S₂. The force field parameters for all standard amino acids were taken to be the CHARMM27 configurations with CMAP corrections.⁴¹ For the Fe clusters and the corresponding amino acid ligands, the force field parameters were taken from the published results of Chang and Kim.⁴² The details of the parameter modifications done in the current work can be found in the supplementary material of the published work by McCullagh and Voth.³⁴ 

All of the initial configurations of the hydrogenases at varying electronic states were taken from snapshots of trajectories from the published results of McCullagh and Voth,³⁴ where the proteins were simulated for at least 50 ns. In these previous simulations, the fully oxidized protein was selectively reduced, by adding electrons to the iron cluster sites shown in Figs. 1(a) and 2. The proteins were completely solvated by waters forming a rectangular simulation box, an additional proton was added to the bulk water near the active site of the hydrogenase, and counter ions were added to neutralize the system charge. The systems were equilibrated again for 2 ns using a Langevin thermostat and barostat set to 300 K and 1.0 bar using the NAMD 2.8 software.⁴³ The same protocol was used for the double mutant simulation except that after the mutant residues were added to the snapshot along with the waters and ions of the A²⁻ electronic state, the protein was minimized using the conjugate gradient for 1000 steps. The system was then equilibrated for 50 ns with the NVT ensemble.

The details of the MS-EVB methodology are described quite extensively in many places^{36–38,40,44} and to provide a general introduction to the method, a brief description is included here. For more information about the force field, the reader is encouraged to look at the other publications on the topic. In the present work, the MS-EVB3.1 model, described in the supplementary material of Ref. 45, was used.

The MS-EVB model is a molecular mechanics scheme in which bond topologies are not fixed, thus allowing for chemical reactions, such as in this case of proton transport, throughout the MD simulation. To do this, the bonding topologies, labeled by |i⟩, represent diabatic states in the MS-EVB Hamiltonian, given by

where |$h_{ij} ( {\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r} } )$|$hij(r⇀)$ and |$\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r} $|$r⇀$ represent the matrix elements and the nuclear degrees of freedom, respectively. The matrix is dynamically constructed from the dynamically selected diabatic states and diagonalized at each MD timestep. The off-diagonal elements of the MS-EVB Hamiltonian describe the coupling between the diabatic states and are parameterized to reproduce ab initio potential surfaces and experimental properties.

Also, in order to better describe the protonic charge defect of the excess proton, the center of excess charge (CEC) variable is used, given by

where |$c_i^2 ( {\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r} } )$|$ci2(r⇀)$ is the contribution of diabatic state i to the ground state potential energy surface, N_EVB is the total number of diabatic states, and |$\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r} _i^{{\rm COC}} $|$r⇀i COC$ is the center of charge vector, given by

where q_k is the partial charge on the kth atom. In short, the center of charge, |$r_i^{{\rm COC}} $|$ri COC$⁠, can be thought to represent the position of the hydronium ion for the diabatic state i, while the CEC coordinate tracks the full excess proton charge defect, which can be delocalized by the Grotthuss mechanism. The MS-EVB simulations were performed using the Rapid Approach for Proton Transport and Other Reactions (RAPTOR)^38,40 add-on package of the LAMMPS⁴⁶ MD code.

The free energy simulations were carried out using the umbrella sampling (US)^47–49 method. The umbrella windows were positioned at equidistant points along a straight line connecting the bulk water to the protein active site. For each electronic state and the double mutant, a US simulation was performed using cylindrical coordinates, where a force constant of 10 kcal mol⁻¹ Å⁻² was applied in the axial direction and a smaller potential of 0.1 kcal mol⁻¹ Å⁻² was applied in the perpendicular directions. As such, the excess proton CEC was mostly unhindered in the perpendicular direction and was able to explore a large radial area. The initial US simulations were “path-finding” and even an arbitrary, curvilinear path was able to be identified. A separate US simulation was performed again for the calculations of the free energy also with cylindrical coordinates. The weighted histogram analysis method (WHAM)⁴⁷ was used to construct the potential of mean force (PMF) from the umbrella windows. The PMF is the free energy profile for the proton transport via both standard diffusion and Grotthuss shuttling.

In order to better understand the nature of the proton transport pathways in WC2, water density maps were constructed from equilibration simulations at each electronic state of the [FeFe]-hydrogenase protein. The maps were generated as Gaussian densities of the frequency that the oxygen atoms in water molecules of WC2 occupy points on a 0.5 ų three-dimensional grid throughout the equilibration trajectory. The figures represent density isovalues of 3.0. These are shown in Fig. 3. The density in state A²⁻ (top panel in Fig. 3(a)) is the most pronounced, with large pockets of water in several locations in WC2. A water wire leading to the proximal Fe of the active site is also present. In fact, a large pocket of water molecules is observed between Ser232 and Gly502 (Fig. 1(b)). A clear contrast is observed in electronic state B²⁻ (Fig. 3(b)), where there is almost no water density within the water channel. A small pocket of water is located near Val249, Thr250, and Val529. In state C²⁻, shown in Fig. 3(c), this pocket is shifted closer to the active site. Somewhat surprisingly, there exist two other small regions of water density near the proximal Fe and FeS cuboid in state C²⁻, which are obstructed from the active site by Cys503 and Ser232. Finally, the bottom panel of Fig. 3(a) shows the water density map of the double mutant G418F/G502V in state A²⁻. Similar to the wild type protein, many water molecules fill WC2; however, at any given snapshot of the double mutant trajectory, there was no longer a water wire leading from the bulk to the enzyme active site.

For each electronic state considered, an initial “path-finding” US simulation was performed to generate a pathway for the excess proton CEC connecting the active site to the bulk. Because comparisons were only made between pathways in different electronic states, the path-finding US method was deemed sufficient, as opposed to a more elaborate metadynamics-based approach.³¹ To this point, the pathways were not completely unbiased since some aspects of the paths were known a priori. In each simulation the hydrated proton (hydronium-like species) was originated in bulk water. In each electronic state A²⁻, B²⁻, and C²⁻, the excess proton CEC was initially positioned near residues Thr250 and Gln533, about 20 Å away from the active site. Here, the systems were equilibrated once again for an additional 50 ps at 300 K. The umbrella windows were then evenly spaced in each of the different states (40 windows for electronic state A²⁻, B²⁻, and C²⁻) along a straight line starting from the active site to the initial hydronium position. Each of the path-finding US simulations consisted of 50 ps of sampling at each umbrella window.

In electronic state C²⁻, a pathway was found as shown in Fig. 4(c). The blue spheres in Fig. 4 indicate the instantaneous excess proton CEC position during the simulation. Here, the hydronium is initially located near residues Asp248, Val249, and Thr250. After equilibration the excess proton CEC is sampled through a pocket of amino acids, including Arg234, Ala415, and Leu530, that hinder access to the active site. The CEC changes direction once it passes Leu530, as a direct path to the H-cluster is obscured by the Pro231 and the large Phe417 residues near the distal Fe (left side in Fig. 4(c)) of the active site. A very curvilinear path is taken to the distal side of the active site by the CEC, where a small opening allows the CEC to approach the metal center. No water wire was found for state C²⁻, even after the initial US simulation. Thus, at times the hydrated proton structure was only solvated by one other water molecule as it maneuvered along the pathway. It also should be noted that the apparent disruption in the CEC pathway in Fig. 4(c) near residues Ala415 and Gln533 is an artifact of the large interval between successive snapshots of the CEC trajectory and the perspective of the image. In fact, this pathway is continuous, although less sampled than other regions.

The pathway exhibited by electronic state B²⁻ (Fig. 4(b)) is similar to that of state C²⁻. Again the hydrated proton moves through a water pocket created by residues Arg234, Ala415, and Leu530 at about 10.5 Å away from the active site. Yet, unlike in state C²⁻, there is more space in this region, and the excess proton CEC can more easily move between these residues. This is also the case closer to the active site. Whereas in state C²⁻ the large phenylalanine and proline residues obscure the access of the CEC to the proximal Fe of the subcluster, here there is just enough space opened up for the CEC to directly access the active site. As Fig. 4(b) shows, the pathway found is not nearly as curvilinear as it makes its way through the protein. As was previously the case, because of the restrictive nature of the path, no water wire existed along the determined path. Instead, only a few water molecules and the excess proton are observed in the obstructing residues.

In electronic state A²⁻, shown in the top panel of Fig. 4(a), a curvilinear path was found, wrapping mainly around Pro231, positioned directly under the [2Fe]-subcluster of the active site as is shown in Fig. 1(b). The excess proton CEC followed the water wire that was initially observed throughout the entire 50 ns production run in the previous work by McCullagh and Voth³⁴ and remained during the equilibration runs of the present study. The Ser232, Arg234, Gly418, Gly502, Cys503, and Val504 residues contour the water wire in WC2. At about 6 Å from the active site and near the Ser232 residue, the water pathway becomes something of a branch point (shown as the purple sphere in the top panel of Fig. 4(a)), where it appears to connect to other nearby water pathways, much in the same way the pathways in WC1 converge near the Cys299 residue. This underlines the importance of the residues in this region in allowing or disallowing proton transport to the active site. The ability of any mutation to obstruct access to the active site is directly related to this branch point. Residues near the mouth of WC2 would have little effect on disrupting the proton transport pathway, and thus mutations were made to the glycine residue near the active site. The efficacy of any mutation to block CEC access to the active site is directly related to its ability in disrupting this branch point near the active site.

In electronic state A²⁻ with mutations to residues G418F and G502V, a proton transport path was also found using the path-finding US simulations, and the instantaneous CEC positions are shown in the bottom panel of Fig. 4(a). The hydrated proton was initially positioned about 18 Å away from the active site, near Arg234. The mutations of the glycine residues near the [2Fe]-subcluster were shown to obstruct the water wire from approaching the proximal Fe atom as was the case in the wild type state A²⁻. This is mostly due to the presence of the large hydrophobic phenylalanine group, but the valine residue also aids in preventing the water molecules from approaching the active site near the bridging Cys503 ligand. Preliminary equilibrations were done with only a single site mutation G418F, and a water wire pathway similar to state A²⁻ was still present. The CEC approaches the active site, instead, on the opposite side of the Pro231 residue towards the distal Fe atom. There is no water wire seen in this case, but there remains enough space for the CEC to maneuver through to the active site between the Asn269, Phe417, Pro231, and G418F residues.

Because in the current simulation the MS-EVB model did not allow for the protonation of amino acids, it could be the case that the CEC might be more stabilized by these neighboring residues along WC2. However, because of the general restricted nature of the pathways in both state B²⁻ and C²⁻, this does not seem very likely. The only protonatable amino acid along the proton pathways is Asp248, and this is located near the mouth of the water channel in reduced cluster states B and C. In electronic state A²⁻, the water wire is easily discernible and no protonatable residues are present near the active site.

Because the loops in electronic states B²⁻ and C²⁻ work to envelop the active site, thus preventing water from approaching the enzyme center, the proton transport pathways in the different electronic states are noticeably different and this is reflected in the PMFs, all shown in Fig. 5. The PMFs of state B²⁻ and C²⁻, shown in Figs. 5(b) and 5(c), respectively, are similar to each other, yet have some distinguishable features. The main characteristic of both of the PMFs is their nearly continuous uphill nature from the bulk situated near Gln533 and Asp248 (located at X = 16 Å) to the center of mass of the metal center (near X = 0 Å). However, the two surfaces are not identical. To start, the PMF in electronic state C²⁻ experiences a much steeper rises as the hydrated CEC approaches the enzyme center. Also there are almost no plateau regions except when the CEC is positioned about X = 13 Å away from the active site. This corresponds with the water pocket shown in Fig. 3(c), formed by Leu530, Arg234, and Val249. In state A²⁻, this pocket opens up and includes water molecules from the bulk. The same plateau region is observed in state B²⁻ but is more pronounced. This agrees with what we would expect given the presence of more water molecules in the Arg234, Ala415, and Leu530 pocket in state B²⁻. Also, unlike in state C²⁻ where the CEC experiences relatively easy access to the pocket near Arg234, in the reduced cluster B state access to this pocket is very much constrained by Asp248, Val249, Thr250, and Gln533, resulting in a steep rise at the right edge of the PMF. This rise could be slightly mitigated by the protonatable aspartic acid residue, but the general trend still stands. As the CEC moves closer to the H-cluster from the pocket, the PMF surface again rises quickly in both states. In state B²⁻ this free energy rise is the result of the CEC moving through the arginine and proline residues by the active site until it levels off in the water pocket near Arg234 at about 9 Å from the active site. The PMF of state C²⁻ rises continuously only as it passes through Thr250 and Gln533. It is at this point in the excess proton CEC pathways of states B²⁻ and C²⁻ that the pathways diverge as was discussed in Sec. III A. As the CEC in the reduced cluster B state follows Pro231 and Ser232 towards the proximal Fe of the active site, the PMF increases steadily. The free energy of state C²⁻ also increases as the CEC approaches the active site, this time from the distal Fe atom. Again, the increase is steady until about 6 Å, after which the free energy rises much more. However, what seems clear is that the pathway taken by the hydrated proton CEC in electronic states B²⁻ and C²⁻ is unfavorable in either case, especially compared to the free energy profile in electronic state A²⁻.

The black line in Fig. 5(a) shows the PMF (free energy profile) for proton transport in electronic state A²⁻. The initial portion of the PMF on the right side of Fig. 5(a) has the excess proton CEC at the mouth of WC2 from the bulk near the pocket created by Arg234, Leu530. This is shown as an initial small uphill climb in the free energy of nearly 1 kcal mol⁻¹ with a peak at 10 Å from the active site. As the CEC bends its way around Pro231, the free energy slowly decreases until it reaches 7 Å from the active site. In this region a small pocket is created by Ser232, Val504, after which the free energy experiences a much more steep decrease until the minimum is reached. The main feature of the PMF of state A²⁻ is its minimum about 3–4 Å away from the center of mass of the active site (X = 0 Å). It is important to note that the radius of the di-iron active site with its carbonyl and cyanide ligands is around 3 Å. Since no proton transfer was allowed in the present MS-EVB model between the excess proton on the water molecules and the di-iron complex as happens in the catalytic cycle, the CEC encounters the wall of the active site. This is shown in the PMF of state A²⁻ as a sharp rise beginning at closer than 3 Å.

The reason for the relatively small free energy barrier in the PMF of electronic state A²⁻ is a direct result of the unbroken water wire leading from the bulk to the active site. It was shown previously³⁴ that the electrostatic potential (ESP) along WC2 became more negative as a function of the inverse distance to the active site. This attractive ESP not only allows for the presence of the existing water wire through the separation of two loops, but it also adds an attractive force which drives the excess proton CEC towards the active site of the enzyme. Once the CEC is past the pocket near the positively charged Arg234, its free energy is completely exergonic until it reaches the H-cluster.

The above results suggest a possible special role for WC2 when the enzyme is in the A²⁻ state, i.e., an electron transfer activated proton transport pathway. In order to explore this behavior and to make an experimentally testable prediction, the proton transport PMF was calculated for electronic state A²⁻ with the double mutation G418F and G502V. The result is shown as the red line in Fig. 5(a). The most noticeable feature is that the transporting of the CEC towards the active site is no longer exergonic after passing by the Arg234 residue. Instead, the free energy becomes quite unfavorable. The free energy rises about 20–25 kcal mol⁻¹ above what was observed for the A²⁻ state of the wild type enzyme. This result suggests that the G418F/G502V double mutant could significantly affect the proton transport behavior of the enzyme and likely its catalytic production of molecular hydrogen. (Unless, of course, the enzyme has robust “self-rescue” capabilities for proton transport pathways that we have not tested here.)

The MS-EVB reactive molecular dynamics simulations were used along with free energy sampling to investigate the proton transport process in a second water channel (WC2) of [FeFe]-hydrogenase. The proton transport in WC2 was explored in three different electronic states, mimicking the electron transport pathway from cluster C to cluster B and then to the cluster neighboring the active site domain, cluster A. The reduced cluster state A was also studied with the double mutation G418F/G502V, to show how the excess proton CEC is able to access the enzyme active site (or not). The path shown for electronic state A²⁻ without mutations, which curves around the Pro231 residue leading to the proximal Fe of the active site, was the only path where a water wire was shown to exist throughout the simulations. A stable pocket was formed by Arg234, Ala235, Ala419, Ala415, and Ile416, which allowed for water molecules to fill this area, below the Pro231 residue. In the other electronic states, this cavity was no longer present. The Arg234 residue was shown to obstruct any waters from filling this area due to the convergence of two loops, one containing Asp248, Val249, and Thr250 and another loop containing Asn532, Gln533, and Asp534. The excess proton CEC traveled through two separate pathways in electronic state B²⁻ and C²⁻. The double mutant proton transport pathway was similar to the one taken in electronic state C²⁻ in which the CEC approached the distal Fe atom from the water pocket located near Arg234.

The free energy surfaces (PMFs) were also explored qualitatively with the US method for the three electronic states. It seems clear that the electronic state of the enzyme plays a large role in the ability of protons to move to and from the active site in WC2. The pathway taken in state A²⁻ was the energetically favorable path, which agrees well with previous assumptions about WC2.³⁴ The electronic states B²⁻ and C²⁻ show much more steep PMFs, which are strongly uphill in nature. When the two glycine residues in WC2 were mutated (G418F/G502V), the proton transport PMF of state A²⁻ more closely resembled that of the other two states with an uphill free energy profile to the active site.

The free energy profile in electronic state A²⁻ lends support to the role WC2 could play in bringing a second proton to the active site for the production of molecular hydrogen. This water channel would necessarily work in concert with the residue pathway and the electron transport pathway. Although, the current work expands our understanding of the possible electronic state-dependent (activated) proton transport through WC2, there still exist many possibilities of further study in this area. For instance, a more thorough examination of the different proton pathways taken in WC2 could be useful, including the relative energy differences associated with the excess proton CEC approaching the proximal or distal Fe of the active site. Such results, in turn, would be tied to how the different catalytic cycles of the active site either facilitate or hinder proton transport in WC2.

The authors acknowledge the University of Chicago Research Computing Center and the Extreme Science and Engineering Discovery Environment (XSEDE) for computation resources. This research was supported by the (U.S.) Department of Energy (DOE) under Contract No. DE-AC02-06CH11357.