Fast Proton Transport in FeFe Hydrogenase via a Flexible Channel and a Proton Hole Mechanism.

Di-iron hydrogenases are a class of enzymes that are capable of reducing protons to form molecular hydrogen with high efficiency. In addition to the catalytic site, these enzymes have evolved dedicated pathways to transport protons and electrons to the reaction center. Here, we present a detailed study of the most likely proton transfer pathway in such an enzyme using QM/MM molecular dynamics simulations. The protons are transported through a channel lined out from the protein exterior to the di-iron active site, by a series of hydrogen-bonded, weakly acidic or basic, amino acids and two incorporated water molecules. The channel shows remarkable flexibility, which is an essential feature to quickly reset the hydrogen-bond direction in the channel after each proton passing. Proton transport takes place via a "hole" mechanism, rather than an excess proton mechanism, the free energy landscape of which is remarkably flat, with a highest transition state barrier of only 5 kcal/mol. These results confirm our previous assumptions that proton transport is not rate limiting in the H2 formation activity and that cysteine C299 may be considered protonated at physiological pH conditions. Detailed understanding of this proton transport may aid in the ongoing attempts to design artificial biomimetic hydrogenases for hydrogen fuel production.

Introduction

Enzyme catalyzed hydrogen production is an inspiration for renewable energy research. Electrocatalytic di-iron complexes that mimic the active site of hydrogenase enzymes are able to produce molecular hydrogen fuel by reducing protons, which is an appealing, more sustainable, alternative for the currently prevalent methods based on steam re-forming of fossil fuels and gasification of coal.[1] In particular, the protein family of [FeFe] hydrogenase enzymes is highly efficient in producing H2 and has tremendous potential to function as a precursor for artificial catalysts.[2] Recent rationally designed di-iron compounds show promising activities, but, so far, none of them is remotely as efficient as the natural enzyme.[3]

The active site of the [FeFe] hydrogenase enzyme, commonly referred to as the H-cluster, is located in a hydrophobic cavity inside the enzyme. The transfer of protons and electrons from the protein exterior to the active site follows separate pathways through the protein scaffold, as illustrated in Figure . The electrons transfer via a series of four Fe4S4 cubane cofactors, denoted FS4A to FSAD and one Fe2S2 cofactor, labeled FS2.[4,5] The proton transfer most likely occurs via a Grotthuss mechanism along hydrogen-bonded networks of amino-acid residues and crystallographic water molecules.[6]

Figure 1

(Left) Cartoon structure of the [FeFe] hydrogenase enzyme. The electron pathway, proton pathway, and catalytic center are highlighted. (Right) Zoom-in on the proposed proton transfer channel outlined by residues Glu282, Ser319, Glu279, and Cys299. Residues Lys571, Glu368, and Glu361 may take part in an alternative pathway for proton transfer (see main text).

Three different proton transport pathways have been proposed for [FeFe] hydrogenases.[7] Pathway 1 consists of two glutamic-acid residues (Glu282, Glu279), a serine residue (Ser319), crystallographic water (H2O-615), and a cysteine residue (Cys299) and is considered the primary proton pathway (see right-hand-side panel in Figure ). Site-mutation studies of the residues along this pathway deactivated the enzyme by 95%,[6,8] which is a very strong indication that this channel is indeed the main pathway for the proton transport to the catalytically active H-cluster. Pathway 2 is proposed to start at surface residue Lys571, includes Ser298 and crystallographic water (H20-594, -668, -675), and ends at Cys299.[9] These two pathways are proposed to carry protons to the di-iron bridging di(thiomethyl)amine (DTMA) ligand of the H-cluster, which can act as a proton relay. A third pathway is proposed to protonate the CN– ligand of the H-cluster. This pathway contains surface residues Glu245, Lys237, Ser202, and Glu240 and three crystallographic water molecules.[10]

The residual 5% activity after mutation of key residues in pathway 1 suggests that the other proposed pathways may serve as backup channels to transfer protons to the catalytic site, albeit with a much lower efficiency. However, recent work by Duan et al. suggests that the residual efficiency may also be due to transfer via water molecules in the site mutated channel.[11] Senger et al. showed that the flexible hydrogen-bond network along the primary proton transfer pathway helps in the bidirectional transfer of protons and that the protons follow a discontinuous mechanism in the oxidized state.[12] Lampret and co-workers recently showed that long-range interactions between the electron and proton transfer pathways affect the directionality and efficiency of the catalysis.[13]

Having established pathway 1 as the main proton channel, still many questions concerning the detailed mechanism have remained unanswered. Is the transfer of the proton through the protein also enhanced by concerted protein rearrangements and/or solvent penetration as seen for the electron transfer? Does the charge transfer follow an excess proton or a hole transport mechanism? Does the hydrogen-bonded chain of amino acids and water molecules along the transfer pathway require a directional reorganization after each proton transfer, and would such resetting of the chain be a bottleneck on the transfer kinetics? And more generally, what can we learn from the proton transfer process in the enzyme that could benefit the rational design of improved artificial catalysts for hydrogen fuel production?

In our previous study of the electron transfer pathway, we used classical force field simulations to show that the protein matrix is highly flexible and allows for a significant protein reorganization and solvent penetration near each of the cubane clusters when it is reduced.[5] The protein and solvent reorganization stabilize the negatively charged state of the cubane and help the migration of the electrons toward the catalytic site. Also in other redox active proteins such water penetration has been observed in concert with electron transfer.[14,15]

In this work, we focus our attention on answering the above questions about the proton transfer mechanism in the main pathway in di-iron hydrogenase. In particular, we employ classical force field molecular dynamics (MD) and hybrid quantum chemical/force field (QM/MM) simulations to unravel the mechanistic details of the transfer and compute the free energy landscape of the proton hopping along the pathway. The free energy profile will reveal possible bottlenecks and metastable intermediate states and whether indeed the Cys299 may be expected to be protonated under physiological pH.

Methods

The initial atomistic structure of the [FeFe] hydrogenase enzyme for our simulations is taken from the crystal structure PDB 3C8Y of the protein databank.[16] We use the Gromacs simulation package[17] for carrying out classical molecular dynamics simulations. The protein force field parameters are taken from the CHARMM27 force field with CMAP corrections.[18,19] The TIP3P water model is used for the solvent.[20] For the FeS clusters, we use the parameters devised by Chang and Kim[21] and further modified by McCullagh and Voth.[4] All metal ions are bound to their ligating atoms via harmonic bond potentials. The protein is solvated with 30000 water molecules and sodium and chlorine ions to create a 0.1 M NaCl solution, thereby neutralizing the protein charge. After an initial energy minimization, the system is equilibriated for 10 ns in the NPT ensemble to a temperature of T = 300 K and a pressure of p = 1 atm. A subsequent equilibration run of 20 ns followed by a production run of 100 ns are carried out in the NVT ensemble. These 120 ns NVT simulation runs are repeated five more times, in which every time a different cubane cluster, or the FS2 cluster, is electronically reduced. The total length of these MD simulations of the hydrogenase system is 720 ns, of which 600 ns is used for analysis. All simulations are performed with the default leapfrog integration scheme and a time step of 2 fs, using LINCS[22] to constrain all bonds. The v-rescale thermostat[23] and the Parrinello–Rahman barostat[24] are used for the isothermal–isobaric system.

The mixed QM/MM simulations are performed using the CP2K software,[25] in which the electronic structure in the QM region is described at the density functional theory (DFT) level of theory using the mixed Gaussian and plane wave method as implemented in CP2K. The PBE exchange–correlation functional[26] is used, augmented with Grimme’s D3 dispersion corrections.[27] The valence electronic wave functions are expanded with the DZVP-MOLOPT[28] basis set together with GTH type pseudo potentials[29,30] to describe the interaction with the core electrons and the nuclei. The plane wave expansion is cut off at 300 Ry. The PBE functional is a GGA-type functional that provides a very good compromise between accuracy and efficiency, which is essential for carrying out our DFT-based MD simulations within a QM/MM framework. PBE has been extensively benchmarked by us for our previous studies on hydrogenase catalysis[31,32] and by others, showing for example the good performance on hydrogen-bonded systems.[33−35] For the MM region, the same CHARMM force field is employed as for the Gromacs simulations. The QM/MM simulations[36,37] are carried out in the NVT ensemble with a QM box size of 15 × 15 × 25 Å3. A Poisson solver is used to avoid spurious periodic electrostatic interactions.[38] The v-rescale thermostat[23] with a time constant of 50 fs maintains an average temperature of 300 K.

Constrained molecular dynamics simulations[39] are carried out to compute the free energy profiles along the proton transfer pathway. Here, a collective variable (CV) that represents the reaction coordinate of the proton transfer process is held fixed during the simulation using the method of Lagrange multipliers. The CVs used in this study are all of the type of the difference between two distances, as follows:in which the first right-hand-side term is the distance between the atom that donates the proton, XD, and the proton, H, and the second term is the distance between the accepting atom, XA, and the proton. X is an oxygen, nitrogen, or sulfur atom. The CV, q, is thus guided from a negative number in the reactant state, via q = 0 at which the proton is equidistant from the donor and the acceptor, to a positive number in the product state.

An equilibration run of 3 ps and a production run of 10 ps is performed at each constrained CV value. The first constrained QM/MM run (see also panel A in Figure ) is started from a configuration taken from the final part of the force field MD simulation of the oxidized hydrogenase, in which the Glu279 and Glu282 residues are oriented in a “forward” configuration for proton transfer and most H-bonds in the channel are intact. Each subsequent constrained QM/MM run was then started from the last frame of the 3 ps equilibration run of the former run, by increasing q to the next constrained value and starting the next equilibration run. The free energy profile as a function of the CV, ΔG, is obtained by integrating the accumulated average constraint force (i.e., the Lagrange multiplier) over the CV, ranging from the initial to the final constrained CV value:Note that the derivative of the Hamiltonian with respect to the CV, dH/dq, is the constraint force, and the subscript q′ denotes that the average is obtained from a constrained ensemble.

Figure 4

Top left: free energy profile of the proton transfer from the cysteine sulfur atom to the nitrogen atom of the DTMA residue. The labels A, B, and C denote the initial, transition, and final states, respectively. Center-left: computed average force of constraint. The CV values at which the average force is zero signify the stable states and the transition state. The error bars signify the standard deviations in the force. Bottom-left: average distances between S–H, N–H, and S–H1 atoms. Right: molecular structures of states A, B, and C. The atoms labels are indicated in structure A.

Results and Discussions

The results are organized as follows. We first discuss the analysis of the 100 ns long force field MD simulation to assess the stability of the hydrogen-bond network of the proton transfer pathway outlined by amino-acid residues Glu282, Ser319, Glu279, and Cys299, as shown in Figure . Next, we examine the effect of the protein oxidation state on the protein structure in the vicinity of the proton transfer channel. The second part of the results is dedicated to the proton transfer free energy landscape and mechanism, computed with constrained QM/MM MD simulations. The transport through this long channel is subdivided into four consecutive proton transfer segments: (1) the transfer to the DTMA bridging ligand at the catalytic site from the cysteine residue Cys299, (2) the transfer to Cys299 from Glu279 via encapsulated water molecules, (3) the transfer to Glu279 from Glu282 via Ser319, and finally (4) the protonation of Glu282 from the solvent. The reason that we start the proton transport process from the solvent to the catalytic site from the end point is because we find that the most favorable mechanism follows a “proton hole” transport mechanism. In other words, the process proceeds by transport of a proton deficiency from the catalytic site to the protein exterior and the solvent.

Hydrogen-Bond Network Stability and Glutamic-Acid Rotation

For each hydrogen molecule to be formed at the catalytic center in the [FeFe] hydrogenase enzyme, two protons have to be transported from the solvent to the H-cluster via the proton transfer channel. Note, however, that after each delivery of a proton via Grotthuss mechanistic jumps along the hydrogen-bonded network in the channel (i.e., after the passing of a “hole” in the opposite direction) all hydrogen bonds have switched direction as the hydrogen-bond donors have become the acceptors and vice versa. Panel C in Figure shows the H-bonded network ready for forward transport. After all protons have made one Grotthuss jump forward, only backward transport is possible unless several structural rearrangements take place to reset the network to the initial “forward state”. In particular, the two glutamic acids have to undergo a 180° rotation and also water molecule H2O615 has to reorient.[40] Clearly, a high free energy barrier for these conformational changes to reset the network to the forward state will slow down the catalytic proton reduction process.

The energy barrier of the glutamic-acid rotation is known to be small, unless the rotamers are stabilized by strong hydrogen bonds to the neighboring residues. However, our 100 ns long force field MD simulation reveals that the H-bond network is not particularly stable, showing large fluctuations in the side-chain orientations of all amino acids involved in the proton transport channel with the H-bonds as often broken as intact. This is illustrated in Figure for the Glu282–Ser319–Glu279 segment. Both acid groups are seen to make spontaneous rotations on the simulation time scale, and also the serine moves in and out of conformations in which it can make hydrogen bonds with the glutamic acids. Interestingly, one might have expected that an intraprotein channel for efficient proton transfer would require a very stable H-bonded network, but that would hamper fast resetting of the network back to the forward configuration. Moreover, our QM/MM simulations suggest that the H-bonded network is (locally) stabilized during the actual proton transfer due to the negative charge of the passing hole. This finding also corroborates a previous study, which shows that the energy barrier of this conformational change is negligible compared to the proton transfer energy barrier.[40]

Figure 2

(Left) Distances of selected hydrogen bonds during the force field MD simulation. The O1–H1 and O2–H2 distances are ≤2 Å in their “forward” conformer state. The distances increase as they change to the “backward” conformer. (Right) Superposition of the OH and O atomic positions of the Glu282, Ser319, and Glu279 residues at small time intervals during the MD simulation shown by transparent red (O) and white (H) spheres. The blue and red ovals group positions belonging to the forward and backward states, respectively. The superimposed ball-and-stick structure illustrates the H-bond network in the forward state.

Water Penetration and the Effects of Cubane Reduction

In our previous study of the electron transport process via subsequent reduction of the Fe4S4 and Fe2S2 cubane clusters (shown in left panel of Figure ), we observed large fluctuations in the protein structure and penetration of solvent water, especially in the vicinity of the cubane clusters that had received an electron.[5] Such outer sphere reorganization is an integral part of Marcus’ theory of electron transfer that helps to stabilize the oxidized and reduced states. Here, we examine the protein scaffold fluctuations and water penetration in the vicinity of the proton transfer channel, and the influence of the location of an electron in the pathway outlined by the cubane clusters FS2, FS4D, FS4C, FS4B, and FS4A, here ordered by closeness to the proton transfer channel, FS2 being farthest away.

Water density maps computed from the classical force field MD simulations are shown in Figure . They reveal several small pockets of water, including the two encapsulated water molecules (H2O615 and H2O616) taking part of the main proton transport channel (encircled with a yellow oval). However, most larger pockets are located at the protein exterior (denoted by blue ovals), and little variation upon cubane reduction is seen. Also no continuous water channel from the protein surface to the H-cluster was observed, through which protons could directly diffuse to the catalytic site. Different from the cubane clusters vicinity, solvent reorganization is seen to be minimal in the neighborhood of the proton transfer channel. Water density maps at t = 50 and t = 100 ns show that the water percolation in the enzyme is established within tens of nanoseconds and that the water pockets remain stable with little variation on this time scale (shown in Figure S1 in the Supporting Information).

Figure 3

Water density maps in the vicinity of the proton transfer pathway at different oxidation states of the five cubane clusters. In the top left panel, all cubane clusters are oxidized; in the other panels, one cubane cluster is in the reduced state as indicated on top. The red isosurfaces represent the average density of water molecules during the simulation at an isovalue of 0.05. The yellow oval denotes the encapsulated water molecules in the proton transfer channel. Blue ovals indicate solvent water in contact with the outermost residues of the different proton pathways.

Proton Transfer from Cys299 to DTMA

When the proton transfer channel is in the initial “rest” state, i.e., no protons are being transferred, the two glutamic acids, the serine, and the cysteine are protonated, and the two encapsulated water molecules taking part in the hydrogen-bonded network are not protonated to hydronium ions. This initial state is shown in panel A of Figure , with all hydrogen bonds in the “forward” configuration, as described above. At the left side of the channel, the unprotonated DTMA ligand is shown, which is considered to function as a proton relay in the catalytic H-cluster that shuttles the protons to the di-iron site. Because of the presence of the lone pair of electrons on the nitrogen, DTMA can act as a “proton hole” to initiate the proton transfer process.[41,42]

The first stage of the proton transport that we model here with constrained QM/MM MD simulation is the transfer of the proton on the nearby Cys299 ligand to the DTMA ligand. The QM region includes all atoms shown in panel A of Figure , whereas the rest of the protein scaffold and its cofactors, the solvent, and the ions are included in the MM representation. The constrained CV is the difference between the S–H and N–H distances, q = dS–H – dN–H. Six constrained simulations are performed with the CV ranging from q = −1.2 Å   in the initial reactant state (labeled A) to q = 0.3 Å   in the product state (C).

The free energy profile shows that this first proton transfer is somewhat uphill by ∼2.5 kcal/mol, with a small barrier of ∼3 kcal/mol (see top-left panel in Figure ). Note that the proton donation by the cysteine residue is accompanied by a spontaneous strengthening of the hydrogen bond from the water molecule H2O615 to the sulfer atom of the cysteine residue (see the green graph in the bottom-left panel). However, creation of the protein hole at the cysteine does not trigger further spontaneous proton transfer reactions in the channel on our short simulation time scale.

Proton Transfer from Glu279 to Cys299

In the second proton transfer step, we set out to transfer a proton to the deprotonated Cys299 ligand from the nearby H2O615 molecule. Note, however, as shown in the panels A, B, and C of Figure , that this reaction is concerted with a spontaneous proton transfer from Glu279 to the water molecule. Therefore, we effectively observe proton transfer from Glu279 to Cys299 via the intermediary H2O615, and thus hole transfer in the opposite direction from Cys299 to Glu279. The constrained CV is q = dO1–H1 – dS–H1. The resulting free energy profile is very flat, with an even smaller reaction barrier than for the first transfer of ∼1 kcal/mol. The measured very low activation barrier and minute overall free energy increase indicates that indeed the second transfer from Glu279 to the H2O615 is spontaneous and reversible, although this transfer was not directly biased through our CV.

Figure 5

Top-left: free energy profile of the proton transfer from Glu279 to the cysteine via intermediary water molecule H2O615. The labels A, B, and C denote the initial, transition, and final states, respectively. Center-left: computed average force of constraint. Bottom-left: selected average distances. Right: molecular structures of states A, B, and C. The local hydrogen-bond network is highlighted with green lines in structure B; the average hydrogen-bond distances are shown in red.

An alternate proton transfer pathway has previously been proposed in which residue Glu279 does not take part in the transfer process.[40] In that case, when the proton is transferred to Cys299 from the H2O615 molecule, another proton is transferred to H2O615 from the H2O616 molecule. But our results indicate that such a pathway is less favorable. The spontaneous transfer of a proton from Glu279 and the subsequent stabilization of the intermediate structure by the hydrogen-bond network is more likely. Moreover, previous experimental work has shown that mutating the Glu279 residue will reduce the hydrogen production rate by a factor of 1/2500,[6] thus illustrating the importance of this residue in the proton transfer process. Our results show that the second water molecule aids in forming a favorable hydrogen-bond network with the unprotonated Glu279 residue, thereby stabilizing the negative charge on the residue and enhancing the proton transfer process. This H-bond network is illustrated by green lines in structure B in Figure .

Proton Transfer from Glu282 to Glu279

The third segment of the proton transport process is the transfer of a proton from Glu282 to the Glu279 via the intermediary serine residue. We perform a series of eight constrained QM/MM MD simulations, in which the QM region includes only the atoms of these three residues, as shown in panel A of Figure . Also here, we control the first transfer step from Ser319 to Glu279 and observe that simultaneously the second protonation from Glu282 to Ser319 occurs spontaneously. The constrained CV is q = dO1–H1 – dO2–H1. The free energy profile shows that also this concerted transfer is somewhat endergonic by ∼1 kcal/mol. The free energy barrier is, with ∼5 kcal/mol, the largest barrier, and thus rate limiting, although it is still a very small barrier.

Figure 6

Top-left: free energy profile of the proton transfer from Glu282 tot Glu279 via Ser319. The labels A, B, and C denote the initial, transition, and final states, respectively. Center-left: computed average force of constraint. Bottom-left: selected average hydrogen-bond distances. Right: molecular structures of states A, B, and C, with atom labels indicated in structure A.

Protonation of Glu282 from the Solvent

In the constrained QM/MM simulations of the previous proton transfer segment, the proton hole arrived at the glutamic acid 282, which is located at the protein surface in contact with the water solvent. To complete the proton transfer process, here we examine the protonation of the Glu282 residue from a proton donor in the solvent. The Glu282 residue and all water molecules within 4 Å of this residue are taking part of the QM region. For the proton donor, we simply introduce a hydronium ion, H3O+, by protonating one of the water molecules. After only 3 ps of (unconstrained) QM/MM MD simulation, the proton from the hydronium ion is seen to jump to the Glu282 residue. Figure shows two simulation snapshots before and after the proton transfer. This shows that the proton transfer from the solvent to the deprotonated Glu282 is spontaneous and can be observed within the picosecond time scale of our DFT-MD simulation.

Figure 7

Snapshots from the 3 ps QM/MM simulation showing the protonation of the Glu282 residue from a nearby hydronium ion. The proton is shown in orange. Initially, the proton is at the solvent (A), and after 3 ps, the proton is transferred to the Glu282 residue (B).

Conclusions

In this study, we have used force field and hybrid QM/MM molecular dynamics simulations to probe the proton transfer in a di-iron hydrogenase enzyme. The force field simulations show that the protein scaffold does not allow for significant water permeation in the neighborhood of the main proton transfer channel, other than the two incorporated crystallographic water molecules taking part of the channel itself, and minor fluctuations at the protein surface. This is in contrast with the substantial protein reorganization and water penetration seen previously in the vicinity of the cubane cofactors that line out the electron transfer pathway in the hydrogenase.[5]

On the other hand, we observe significant fluctuations in the protein transfer channel hydrogen-bond network itself. We believe that this flexibility is important to allow for rapid rotation of the glutamic-acid rotamers and resetting of the hydrogen-bond directions in the forward direction after each proton transfer. The actual proton transport through the channel follows a “proton hole” mechanism, in which a proton deficiency travels in the opposite direction all the way from the catalytic site to the protein surface, while all protons only make a single Grotthuss jump forward to the next amino acid or water molecule in the channel.

The overall process of proton transport in the di-iron hydrogenase in the oxidized state, that is, the process to transfer a proton hole from the DMTA ligand at the catalytic site to Glu282 at the solvent exposed protein exterior, is slightly endergonic by 4 kcal/mol. The rate-limiting free energy barrier for the proton transfer, computed with constrained QM/MM MD simulations, is only ∼5 kcal/mol. This rather low barrier thus allows for fast proton transfer kinetics. The flat free energy landscape agrees with the experimentally observed bidirectionality and reversibility of the process. The conformational changes of the glutamate residues are important for maintaining a favorable hydrogen-bond network. Classical molecular dynamics simulation shows that the energy barrier for the conformational changes is negligible when compared to the proton transfer steps. However, a consequence of the highly flexible hydrogen-bond network is that the observed mechanism is stepwise rather than a fast delocalized proton shuttle along a continuous wire spanning the entire channel. Formation of local hydrogen-bond interactions helps in the proton transfer process by stabilizing the excess charge on the residues that act as a proton acceptor. One of the two crystallographic water molecules present in the cavity between residues Cys299 and Glu279 is involved in the proton transfer reaction, whereas the other water molecule stabilizes the local hydrogen-bond network to facilitate the proton transfer. Our study also shows that the most favorable proton transfer route is via Glu279, thus explaining the observed slowdown of the hydrogen production when mutating the Glu279 residue. These findings may also be helpful for the ongoing efforts to design artificial di-iron hydrogenase catalysts in porous materials for hydrogen fuel production.[43]