Article pubs.acs.org/JPCB
Proton Transport in Clostridium pasteurianum [FeFe] Hydrogenase I: A Computational Study Hai Long,* Paul W. King, and Christopher H. Chang National Renewable Energy Laboratory, MS ESIF301, 15013 Denver West Parkway, Golden, Colorado 80401, United States S Supporting Information *
ABSTRACT: To better understand the proton transport through the H2 production catalysts, the [FeFe] hydrogenases, we have undertaken a modeling and simulation study of the proton transfer processes mediated by amino acid sidechain residues in hydrogenase I from Clostridium pasteurianum. Free-energy calculation studies show that the side chains of two conserved glutamate residues, Glu-279 and Glu-282, each possess two stable conformations with energies that are sensitive to protonation state. Coordinated conformational changes of these residues can form a proton shuttle between the surface Glu-282 and Cys-299, which is the penultimate proton donor to the catalytic H-cluster. Calculated acid dissociation constants are consistent with a proton relay connecting the H-cluster to the bulk solution. The complete protontransport process from the surface-disposed Glu-282 to Cys-299 is studied using coupled semiempirical quantum-mechanical/ classical-mechanical dynamics. Two-dimensional free-energy maps show the mechanisms of proton transport, which involve Glu279, Ser-319, and a short internal water relay to connect functionally Glu-282 with the H-cluster. The findings of conformational bistability, PT event coupling with pKa mismatch, and water participation have implications in the design of artificial water reduction or general electrocatalytic H2-production catalysts.
■
INTRODUCTION Dihydrogen gas (H2) is an energy carrier that is widely used in industry and has been identified as a potential fuel for a carbonconstrained world.1,2 In the pursuit of catalysts that can reduce protons to H2, one consideration is the rate at which protons can be supplied to catalytically active centers. Although proton transfer (PT) functionality as a catalytic ligand design feature has already begun to receive attention, particularly with respect to accelerating electrocatalytic H2 production,3,4 evolutionarily optimized biocatalysts may offer guidance in the design of more highly tuned proton reduction catalysts.5,6 (We use “proton transfer” (PT) to imply the elementary chemical step, “proton transport” to imply the directional and collective net process arising from multiple transfers, and “PT pathway” to imply the structural feature(s) that host transport via multiple transfers.) In enzymes, complex PT pathways involve amino acid side chains as well as protein-bound water, and the dynamic motion of the protein is a crucial component in the overall proton transport process.7 It is commonly assumed that protons in enzyme are transported by the classical stepwise Grotthuss mechanism. In this mechanism, protons “hop” stepwise from the proton donor at the starting point to the proton acceptor at the end point with hydronium ions or acidic form of amino acid side chains formed in the intermediate states. However, if the proton acceptor at the end point has a more basic pKa value, a “proton hole” mechanism is also possible, in which the proton hops in a reverse manner with hydroxide or basic side chain formed in the intermediate states.8 Studies of two model enzyme systems illustrate certain PT features that may be generally important to biological, and biomimetic, catalysis. © 2014 American Chemical Society
Human carbonic anhydrase II (CAII) has a primary PT path involving a network of water molecules and a bistable histidine side chain that acts as a shuttle,9−12 although alternative routes may exist nearby.13 Computational simulations on the PT pathway in CAII indicated that a proton might transport using the “proton hole” mechanism.8,12,14 Cytochrome c oxidase involves a more complex catalytic cycle than CAII. A detailed kinetic study of this system has illustrated that protons are transported through gating mechanisms,15 both through bistable switching via the Glu-242 side chain by taking protonation state-dependent conformational changes,16,17 and through rearrangement of an internal water chain.18 In both of these systems, the concepts of bistable switching and internal water chains are found. For hydrogenases, which catalyze the interconversion between separated protons and electrons and hydrogen gas,19 a complete understanding of the PT process is still emerging. Hydrogenases are candidates as catalysts in water-splitting photosynthetic20 and H2 fuel-cell devices.21,22 Three classes of hydrogenases are defined by their active site metal complement. [NiFe] and [FeFe] hydrogenases contain binuclear metalloclusters, while the Fe-only hydrogenases contain a mononuclear iron center and catalyze hydride transfer to an organic cofactor.19 In vivo, [NiFe] and [FeFe] hydrogenases show a preferred direction of turnover based on the relative electrochemical potentials of redox centers, with [NiFe] hydrogenases Received: August 28, 2013 Revised: January 7, 2014 Published: January 9, 2014 890
dx.doi.org/10.1021/jp408621r | J. Phys. Chem. B 2014, 118, 890−900
The Journal of Physical Chemistry B
Article
generally catalyzing H2 uptake and [FeFe] hydrogenases catalyzing proton reduction.23 While several PT pathways have been proposed24 and investigated using computational methods25−29 for both of these enzymes, we have focused on describing in greater detail the PT process in model [FeFe] hydrogenase Clostridium pasteurianum [FeFe] hydrogenase I (CpI). To date, five X-ray crystallographic structures have been determined for different [FeFe] hydrogenases, representing various forms of CpI (PDB IDs 1FEH,30 3C8Y,31 and 1C4A32), periplasmic [FeFe] hydrogenase (DdH) from Desulfovibrio desulf uricans (PDB 1HFE33), and the apo-form of HydA1 lacking the [2Fe]H center from the eukaryotic green alga Chlamydomonas reinhardtii (3LX434). Although each enzyme comprises different electron-transfer cofactors and differs in tertiary structure, the proteinaceous region that coordinates the H-cluster shows high sequence and structural conservation among the [FeFe] hydrogenases.24 Figure 1a shows the
Figure 2. Proton-transfer pathways in CpI [FeFe] hydrogenase. The two pathways in CpI described in the text are labeled with green (Pathway A) and blue (Pathway B) arrows. Water oxygen atoms appear as red spheres, and peptide atoms are in stick representation. Hydrogen atoms have been omitted for clarity.
Pathway A includes several protein-coordinated water molecules and residues Glu-361 and Lys-358 and ends at the CN ligand of the distal Fe atom. (The residue numbers correspond to the CpI polypeptide sequence.) The importance of this pathway to the overall hydrogenase reactivity is still an open question and how the proton is transferred from Lys-358 to the H-cluster remains unknown. In this manuscript, we do not consider Pathway A further. The second pathway (Pathway B) transports protons to the H-cluster through residues Glu-282, Ser-319, and possibly Glu-279. These residues form a cavity referred to as the “central cavity”, and two water molecules reside inside it. This cavity is separated from the H-cluster by Cys-299 and backbone amide groups of Pro-324, Phe-417, and Val-423. Protons have been proposed to be transferred from Ser-319 by either water, Glu-279, or both to Cys-299 and then onto the DTMA ligand’s central N atom.24 All of the residues listed above are conserved in all known [FeFe] hydrogenases. Experimental site-directed mutagenesis studies have shown that mutations of Glu-279, Glu-282, Cys-299, or Ser-319 in CpI and CaI (Clostridium acetobutylicum [FeFe] hydrogenase I, 71% amino acid-identical and 84% -similar to CpI) slow down H2 evolution and oxidation.39−41 For example, all of the mutants of CpI involving these four residues made by Cornish et al. result in a >95% loss of H2 evolution.39 E279D, E279L, C299A, and C299S mutants of CpI have only 1/2500, 1/100, 1/5000, and 1/2500, respectively, of the H2 evolution rate of wild-type CpI.39 Hong et al. used QM/MM calculations to simulate PT trajectories along Pathway B in DdH, and the results agree that these four residues are important in proton transport.42 However, in Hong’s paper, dynamic contributions to the equilibrium protonation states for residues of Pathway B, which would play an important role in the PT process, have not yet been examined in detail. In addition, the free-energy barriers for PT along this pathway remain to be studied. In this study, we have investigated in greater detail the PT mechanism along Pathway B by both classic MD and QM/MM methods. We utilized MD simulations to investigate the conformations of the Glu-279 and Glu-282 residues. Freeenergy differences among different conformations were calculated using the umbrella sampling method. Acid dissociation constants (pKa) for Glu-279 and Glu-282 were calculated by alchemical free-energy perturbation (FEP). QM/
Figure 1. Structure of the CpI hydrogenase and its catalytic center. (a) Structure of CpI derived from PDB entry 3C8Y. The polypeptide backbone is shown as blue ribbon; metalloclusters are shown as van der Waals spheres. (b) Structure of the [4Fe4S]H-[2Fe]H H-cluster, with dithiomethylamine as the bridging ligand; hydrogen atoms have been omitted. An oxygen atom near the distal Fe, tentatively assigned as a water molecule,22 may either be bound to the distal Fe or hydrogen-bonded to the DTMA nitrogen atom (dashed lines). The color scheme for Figures in this manuscript: C: cyan; O: red; H: white; N: blue; S: yellow; Fe: green.
structure of CpI, and Figure 1b shows the H-cluster from CpI. The two Fe atoms of the binuclear center are bridged by a (κ2,κ2)-bound dithiolate ligand and a (μ2-C) carbonyl ligand. In addition, each iron atom has terminally bound CO and CN ligands. The central bridging atom in the dithiolate bridge (DTMA) is believed to be nitrogen.35−38 The proximal Fe atom of the H-cluster is attached to a cysteinate S atom, which also binds the [4Fe4S]H subcluster. The distal Fe atom of the H-cluster in the X-ray structure of oxidized CpI is coordinated definitively by five ligands (2 COs, 1 CN, and a dithiolate), with a proximal crystallographic water molecule either weakly bound to the distal Fe as a sixth ligand30 through a long Fe−O bond, hydrogen-bonded to a central atom in the dithiolate bridge, or both. In either case, a rearrangement with modest energetic requirements would leave a single vacant site on the distal Hcluster Fe atom that could bind to a hydrogen atom during the catalytic cycle. For [FeFe] hydrogenase, two PT pathways have been proposed following analysis of the crystal structures of CpI and DdH (Figure 2).24,30 Starting from the exterior of the protein, 891
dx.doi.org/10.1021/jp408621r | J. Phys. Chem. B 2014, 118, 890−900
The Journal of Physical Chemistry B
Article
way to determine pKa values and equilibrium protonation states. For the glutamate residues of interest, that is, Glu-279 and Glu-282, atomic charges on a protonated Glu (Glh) residue were gradually changed into those of deprotonated Glu (Glu−) during MD simulation, reducing the charge on the acidic proton gradually to 0.0. In this way, ΔG between Glu− and Glh states can be computed using MD trajectories. However, the Glh side chain is neutral, whereas Glu− hosts a single formal negative charge. This presents a challenge for periodic boundary conditions, where a neutral simulation box is required to avoid infinite charge and energy. To overcome this problem, a uniform background charge was applied to the simulation box to neutralize the overall molecular charge. The ΔG calculated from the FEP/MD simulations was thus not the absolute ΔG for Glh deprotonation. To cancel out this artificial effect and get the corrected free energies, we used a thermodynamic cycle, as shown in Scheme 1, similar to the
MM dynamics simulations were performed to obtain the 2D free-energy maps for PT along the pathway.
■
METHODS MD Simulations. The reference structure used in this study was PDB 3C8Y for CpI hydrogenase31 but with nitrogen remodeled as the central atom of the bridging dithiolate ligand of the H-cluster (Figure 1B). Amber force fields and atomic charges were used for protein residues.43 The force fields and atomic charges for the H-cluster and the Cys residues that bond to the H-cluster and [FeS] clusters were based on the values from Chang and Kim.44 Because we investigated the PT from Glu-282 only through Cys-299, the distal water was not included in the simulations. Test runs using both the [FeIFeII]H Hox and [FeIFeI]H Hred valence states showed no apparent distinctive proteinaceous conformational differences, so the Hcluster was kept in the oxidized state throughout, except for the QM/MM studies described later. Tight mechanistic coupling between electron transfer and PT was not considered here.28,29 The protonation states of protein residues were obtained by the PDB2PQR package45 at pH 5.0, close to the reported 3C8Y crystallization conditions (pH 5.1). Solvation energies for different protonation states were estimated using the Poisson− Boltzmann equation and pKa values calculated accordingly. The choices of protonation states for Glu-279, Glu-282, and Cys299 were different in different types of simulations and are summarized in Table S1 in the Supporting Information. For the MD equilibration simulations of this section, because Cys-299 exhibits a pKa of 12.8, it was left protonated. For the later QM/ MM studies, Cys-299 was necessarily considered in its deprotonated state to provide the necessary PT sink. The appropriately protonated enzyme was then solvated in a TIP3P water box, maintaining at least 10 Å between any protein atoms and the box boundaries. Molecular dynamics was simulated using the Amber 10.0 package,46−48 with a time step of 1 fs and periodic boundary conditions. Long-range electrostatics was treated with particle mesh Ewald summation using mesh density of ∼1 point/Å in each cell direction. The system was equilibrated at 298 K and 1 atm pressure, as previously described.49 Calculations were performed on the NREL’s Red Mesa supercomputer system hosted at Sandia National Laboratory. Umbrella Sampling of Conformers. To calculate Gibbs free-energy differences (ΔG) between amino acid conformers, the umbrella sampling method was used. To do this, the corresponding torsion angle was constrained to the values at the initial state by a force of 50 kcal mol−1 degree−2, and MD simulations were run for 100 ps, recording the torsion angles during the simulation. The torsion angle was then progressively constrained to larger or smaller values in 5° increments for successive 100 ps simulations until the target dihedral angle was reached. The stored torsion angle information from different windows was then connected together using the weighted histogram analysis method (WHAM),50−52 and ΔG between the initial and final states was computed. Simulation and analysis parameters were chosen so that the forward (increasing dihedral angle) and backward ΔG estimates were converged. Standard deviations estimated using WHAM and forward/ backward calculations were found to be
