PCCP View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 24026
View Journal | View Issue
Modulation of active site electronic structure by the protein matrix to control [NiFe] hydrogenase reactivity† Dayle M. A. Smith,*a Simone Raugeia and Thomas C. Squierb Control of the reactivity of the nickel center of the [NiFe] hydrogenase and other metalloproteins commonly involves outer coordination sphere ligands that act to modify the geometry and physical properties of the active site metal centers. We carried out a combined set of classical molecular dynamics and quantum/classical mechanics calculations to provide quantitative estimates of how dynamic fluctuations of the active site within the protein matrix modulate the electronic structure at the catalytic center. Specifically we focused on the dynamics of the inner and outer coordination spheres
Received 6th August 2014, Accepted 30th September 2014
of the cysteinate-bound Ni–Fe cluster in the catalytically active Ni-C state. There are correlated
DOI: 10.1039/c4cp03518f
the electron affinity at the active site and the proton affinity of a terminal cysteinate. On the basis of
movements of the cysteinate ligands and the surrounding hydrogen-bonding network, which modulate these findings, we hypothesize a coupling between protein dynamics and electron and proton transfer
www.rsc.org/pccp
reactions critical to dihydrogen production.
Introduction The increased demand for energy from intermittent renewable sources such as solar and wind power requires the ability to store the excess energy produced and retrieve it on demand. In this context, electrocatalysts for rapid and efficient interconversion between electrical energy and chemical energy will be of critical importance. H2 production and oxidation represent central reactions that can be used for energy storage. Hydrogenase enzymes have been widely studied in order to gain insight into H2 activation by metal centers, and coupled electron and proton transfer reactions. This family of enzymes catalyzes the reversible oxidation of dihydrogen at a bimetallic active site ([FeFe] and [NiFe] hydrogenases) or a monometallic ([Fe]-only hydrogenases) active site. [FeFe]-hydrogenase enzymes exhibit, in general, superior rates of hydrogen evolution, but are irreversibly inactivated by the presence of even trace amounts of oxygen. In comparison, the active site structure of [NiFe]hydrogenase enzymes are reversibly inactivated by molecular oxygen, forming catalytically inactive ready (Ni-B) and unready a
Pacific Northwest National Laboratory, P.O. Box 999, MSIN J4-33, Richland, Washington 99352, USA. E-mail: [email protected] b Western University of Health Sciences, 200 Mullins Drive, Lebanon, Oregon 97335, USA † Electronic supplementary information (ESI) available: [NiFe] per-residue rootmean-squared fluctuations, solvent-accessible surface area and secondary structure representations are tabulated in Table S1. [NiFe] trajectories can be obtained by contacting the corresponding author. See DOI: 10.1039/c4cp03518f
24026 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
(Ni-A) states that respectively possess a bridging hydroxide and peroxide between the metal centers.1 The resistance of the active site structure of [NiFe]-hydrogenases to molecular oxygen,2 such that reductive reactivation occurs following inactivation, provides a central design principle important to the creation of robust molecular catalysts. As the active site structure is maintained between all [NiFe]-hydrogenases, the commonly studied enzyme from Desulfovibrio fructosovorans (1YRQ.pdb)3 provides a solid basis to understand the design principles of the hydrogenase protein that allows catalysis at substantially lower overpotentials than observed for currently available molecular catalysts.4 [NiFe] hydrogenase enzymes (Fig. 1a) are heterodimers, consisting of a 34 kDa small subunit and 61 kDa large subunit. The small subunit contains three iron–sulfur (Fe–S) clusters forming an electron transfer chain that shuttles electrons from cellular redox donors (e.g., membrane bound cytochromes) to the bimetallic active site in the large subunit. The active site consists of a nickel center, an iron center, three diatomic ligands (CO and CN ) bound to the iron, two cysteinates (deprotonated cysteines) bridging the two metal atoms, and two cysteinates bound only to the nickel atom (terminal cysteinates). Fig. 1b shows the active site structure including hydrogen bonded groups as inferred from the X-ray crystal structure from Desulfovibrio fructosovorans.3 The inner coordination sphere is defined as the active site and cysteinate atoms Sg and Cb. Cys75 and Cys546 bind to both Fe and Ni atoms (bridging cysteinates CysB1 and CysB2) and Cys72 and Cys543 bind only to
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
Fig. 1 (a) Structure of the ready, oxidized [NiFe] hydrogenase from Desulfovibrio fructosovorans PDB ID: 1YRQ. Protein residues are shown as ribbons (cyan: small subunit; orange: large subunit), metal clusters are shown as van der Waals spheres, and hydrogen bonded residues are shown as sticks. (b) [NiFe] active site, bound cysteinates and hydrogen bonded protein residues. The inner sphere (QM) region is drawn with sticks; residues hydrogen bonded to the NiFe core are drawn as lines; atoms hydrogen bonded to the QM region are spheres.
Ni (terminal cysteinates CysT1 and CysT2). Ser499 hydrogen bonds to one cyanide ligand (CN1), Arg476 hydrogen bonds to the other cyanide (CN2), and His79 hydrogen bonds to the sulfur atom of bridging cysteinate CysB2. There are also hydrogen bonds between bridging cysteinate backbone NH groups and their nearest terminal cysteinate sulfur atoms, which have the potential to drive rotation of cysteine sidechains, as described in Results. Given the high degrees of conservation between the active site structures of all [NiFe] hydrogenases,3,5–16 differences in catalytic turnover rates for H2 production suggest an important role for the surrounding protein matrix in modulating catalysis. In this respect, Pandelia and coworkers1 proved that the nickel is the redox active metal, and its electronic configuration is sensitive to the geometry and charge state of the cysteine sulfurs, which is modulated by the protein backbone by ‘‘affect[ing] the spin and charge density delocalization to modify the electronic properties of the metal site’’. These observations are consistent with empirical correlations that indicate the geometry around the nickel [i.e., the dihedral angle Fe–Ni–CysT2(Sg)–CysT2(Cb), Fig. 1] are responsible for differences in rates of catalysis, such that a smaller dihedral angle corresponds to a faster rate of hydrogen production.17 Indeed, recent density functional theory (DFT) calculations suggest that the Ni spin density and nickel–sulfur bond order depend on the terminal cysteinate orientation.18 Further, rotation of CysT2 modulates the relative stability of the low- and high-spin Ni electron configuration,19 consistent with suggestions that a terminal cysteinate (i.e., CysT2) acts as the proton acceptor20 and that the surrounding
This journal is © the Owner Societies 2014
protein matrix may control dihedral movements to gate proton transfer reactions within the catalytic core and modulate rates of H2 production. As part of the reaction mechanism, which is likely to involve redox-dependent changes in geometry and flexibility of the atoms around the active site metals, a bridging ligand can bind between the Fe and Ni metal centers. The most important catalytically active redox states of [NiFe] hydrogenase are Ni-SIa, Ni-C, and Ni-R (Scheme 1). The Ni-SIa state contains Fe(II)–Ni(II) and a loosely-bound water molecule on the Ni center; the Ni-C state contains Fe(II)–Ni(III) and a bridging hydride between the two metals; finally the Ni-R state contains Fe(II)–Ni(II) and features two isoenergetic, doubly-protonated structures. In one Ni-R structure, H2 is bound to the Ni ion; the other structure contains a bridging hydride and a protonated terminal cysteinate.21 The equilibrium between Ni-C and Ni-R is therefore a critical step for H2 cleavage or H2 production. Our goal is to elucidate the role of protein dynamics (and distortions of the active site geometry) in modulating the electronic properties of the nickel and terminal cysteinate ligand that receives a proton during the catalytic splitting cycle. To explore how protein fluctuations modulate the electron density at the active site in [NiFe] hydrogenases, we have carried out a combined set of classical molecular dynamics (MD) and hybrid quantum/classical mechanics (QM/MM) calculations on the catalytically active Ni-C state (Scheme 1). Our approach is complementary to prior approaches that used DFT calculations to model the active site. DFT calculations of gas-phase active site models have been instrumental in
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24027
View Article Online
PCCP
Paper
for the coordination sphere of the NiFe core derived from the classical force field for electronic structure calculations.
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Construction of Ni-C [NiFe] hydrogenase in aqueous solution
Scheme 1 Active catalytic cycle between Ni-SIa, Ni-C and Ni-R states of [NiFe] hydrogenase, including the isoenergetic, H2 splitting Ni-R conformations.21
identifying reaction intermediates alongside experimental spectroscopic measurements.1,15,22,23 More recent work has extended the model to include hydrogen-bonded protein residues and shed light on the effects of hydrogen bonding on the electron distribution of the catalytically-important active site atoms.18 Our approach is unique in that we used a customized force field17 and molecular dynamics simulations to collect an in situ ensemble of active site structures for QM/MM electronic structure calculations. We find that hydrogen bonding between bridging and terminal cysteine ligands modulate the dynamics of the cysteinate ligands and, in turn, the electron affinity of the Ni center and the proton affinity of one of the cysteinates. On the basis of these findings, we suggest a dynamic coupling between the protein matrix and the catalytic core that is critical in mediating electron density changes associated with catalysis (i.e., driving the Ni-C to Ni-R reaction necessary for hydrogen production).
Methods Overall approach We performed a sub-microsecond (0.1 ms) MD simulation of the Ni-C state of [NiFe] hydrogenase to sample long time-scale protein fluctuations in aqueous solution by using an accurate force field recently parameterized by us from DFT calculations.17 From this simulation, a statistically representative number of configurations were extracted for subsequent electronic structure analysis. For each configuration an analysis of the electron distribution of the catalytic core was performed using hybrid density functional theory (DFT)/molecular mechanics calculations. The DFT region includes the same atoms as in the modelcluster for the force field parameterization (Fig. 1b), whereas for the surrounding protein environment was treated at the force field level. This choice was adopted to avoid possible artifacts introduced by inconsistencies due to use of structures
24028 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
The initial model of the Ni-C [NiFe] hydrogenase was built starting from a X-ray structure Desulfovibrio fructosovorans at 1.83 Å resolution (protein data bank entry 1YRQ). This structure is in the oxidized (hydroxide-bridging) Ni-B state, and the bridging ligand was replaced by hydrogen to build the hydride-bound Ni-C active site. Histidine protonation states were determined using PropKa program24 and visual inspection. This analysis showed that His184 of the small subunit and His79 in the large subunit are protonated at Ne; His61 and His192 of the small subunit and His204, His210, and His549 of the large subunit are protonated at both Nd and Ne; the remaining histidine residues are protonated at Nd. The protein was inserted into a cubic box with minimum 12 Å distance between the protein and the side of the box, then the box was filled with water molecules. Sodium ions were added for charge neutrality. Molecular dynamics simulation Prerequisite to the classical MD simulations is an empirical force field that reliably describes the geometry and motion of the active site within the protein matrix. The force field parameters we derived and published17 reproduce the geometries and frequencies of gas-phase model clusters of Ni-B and Ni-C active sites (and iron– sulfur clusters) with ligated cysteinates modeled by methionates, and also reasonably reproduce the hydrogen bond distances in the initial protein structure. Using the AMBER0325 force field and TIP3P26 water model with the embedded metal cluster parameters for the Ni-C active site and reduced Fe–S clusters, the initial structure was relaxed employing conjugate gradient energy minimization to a tolerance of 50 kJ mol 1 Å 1 and was then used as a starting point for the MD simulation. Equilibration was done in 250 ps long consecutive runs, in which the coordinates and velocities of each equilibration step were used to initialize each subsequent step. In the first equilibration step, non-hydrogen atoms were restrained using harmonic position restraints with a force constant of 10 kJ mol 1 Å 2, and the temperature was set at 100 K. Next, the temperature was increased in 50 K increments to 300 K; then, restraint forces were reduced to zero in 2 kJ mol 1 Å 2 decrements. Finally, the system was left free to evolve for 100 ns. All ´–Hoover thermoof the simulations were performed using a Nose stat27 and a modified Berendsen barostat.28 A 1 fs time step was used without imposing constraints on the chemical bonds. Electrostatic interactions were calculated using the Particle Mesh Ewald method with a cutoff of 12 Å. Analysis of the root mean square displacement (RMSD) reveals that equilibration is reached in about 20 ns (Fig. 2); therefore average quantities were calculated from 20 ns and 100 ns. All classical mechanics simulations and the subsequent analysis of the trajectory were done using the doubleprecision version of GROMACS 4.0 and related analysis programs.29 Hybrid density functional theory/molecular mechanics calculations Hybrid DFT/MM calculations were carried out to describe the catalytic core on a statistically relevant number of configurations
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
configurations taken from the Ni-C [NiFe] hydrogenase trajectory were used for single-point DFT calculations. A total of 1600 configurations 50 ps apart were extracted from the trajectory between 20 and 100 ns. One set of calculations was run without the protein environment, and another set was run with the full protein environment included using a QM/MM scheme, in order to separate contributions from inner-sphere and outer-sphere fluctuations. The electrostatic coupling between the QM region and MM region was described in ref. 38. All of the QM calculations were done using the NWChem DFT module via the QMMM module using default parameters.39 Dynamical charge density analysis
Fig. 2 Black line, top: root-mean-squared deviation of [NiFe] hydrogenase Ca atoms from the 100 ns MD trajectory, relative to the initial crystallographic structure. Grey line, bottom: root-mean-squared deviation of [NiFe] active site atoms and hydrogen bonded Arg476, Ser499 and His79 residues with respect to the initial structure.
extracted from the force field-based MD simulations. All DFT calculations used the BLYP exchange and correlation functional30,31 and a double-z basis set with diffuse functions.32,33 BLYP/double-z has been used previously to successfully describe structures and experimental magnetic properties of the [NiFe] active site.15 Furthermore, the choice of functional and basis set was dictated by consistency considerations with DFT setup used for the [NiFe] force field parameterization17 based on model active site geometry optimization and harmonic frequency calculations (also the same approach used for the DFT-derived [FeFe] hydrogenase force field34). It is important to remark that our previous calculations showed that Mulliken populations and vibrational frequencies are very similar from the double- and triple-z basis sets,17 and different functionals have been shown to produce similar structural properties in [NiFe] active site models.35 As with other transition metal complexes, the high-versus-low spin energy difference of the nickel ion is sensitive to the choice of functional, relative to more accurate post-Hartree–Fock results;36,37 however, the inadequacy of BLYP in this regard is a separate issue from results presented here since our approach explores the potential energy surface of the more stable low-spin state. Since we are using a combined MD/QMMM approach to model the dependence of net electron densities on active site fluctuations within the protein, consistency with the classical force field is paramount. The QM region was selected to include the active site, where the cysteinates coordinated to the metal centers were treated as methionates, again for consistency with the cluster DFT calculations employed for the force field parameterization. [NiFe] hydrogenase
This journal is © the Owner Societies 2014
Mulliken charges and spin densities were used as descriptors of the electron density as a proxy for chemical reactivity (Scheme 1). Specifically, Ni electron deficiency (positive atomic Mulliken charge on the Ni ion) was used to describe the relationship between active site dynamics and Ni(III) - Ni(II), based on the identification of nickel as the redox-active metal,1 and excess CysT2(Sg) electron density was used to describe CysT2(Sg) CysT2(Sg–H), based on previous quantum chemical calculations that found a correlation between Mulliken population and pKa.40 Covariance analysis and functional mode analysis method (FMA)41 were used to calculate the modes of motion that best describe fluctuations of selected Mulliken populations. FMA calculates a single mode of motion representative of a timedependent observable property, in this case electron densities. For this analysis we used all covariance eigenvectors and snapshots between 20 and 80 ns for model building and 80–100 ns for cross-validation. The model building and cross-validation correlation coefficients (RM and RC) are the assessment quantities used to validate that the resulting motion is maximally correlated to the quantity of interest.
Results and discussion Flexibility of the protein backbone and residues near the catalytic core The RMSD of the protein backbone with respect to the starting crystal structure (Fig. 2) indicates that the protein is relatively rigid, as discussed in previous reports, with fluctuation of about 0.1 Å around an equilibrium average value 1.5 Å. Analysis of the root mean square fluctuation of each single residue reveals that the flexible protein elements are residues near the solvent accessible surface and residues in unstructured regions such as chain termini (per-residue RMSF, solvent-accessible surface area and secondary structure classifications are tabulated in the ESI†). The fluctuations apparent in the RMSD of the inner and outer coordination sphere atoms (Fig. 2) are primarily due to His79 and the terminal cysteinates, which are known to be important for regulating spin density.18 Comparison of Mulliken populations and contributions from the inner and outer coordination spheres Before discussing the relationships between fluctuations and electron distributions, it is useful to show how the geometries
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24029
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PCCP
Paper
and Mulliken populations compare using the isolated cluster model and the same atom selection extracted from the MD snapshots for purposes of assessment and also to see how fluctuations of the active site within the protein matrix influence the electronic distribution at the catalytic center. The Ni, CysB2(Sg) and CysT2(Sg) atoms are the most important atoms for consideration because of their known sensitivity to hydrogen bonding [CysB2(Sg)] and direct participation in conversion to the Ni-R state when dihydrogen splitting occurs [Ni and CysT2(Sg)].21 Therefore, our discussion will focus mainly on these centers and, for sake of completeness, on the hydride and the Fe center. In the isolated cluster model, the bulk of the spin density is on the Ni ion (0.64e). Additional spin density is on one of the bridging cysteinate sulfur atom CysB2(Sg) (0.20e) and the terminal cysteinate sulfur atom CysT2(Sg) (0.06e). A very small amount of spin density is present also on the Fe atom (0.02e), which is in the low-spin Fe(II) oxidation state. A previous study, which used the same DFT functional and a similar basis set, reported 0.51e and 0.29e for the Ni and CysB2(Sg) spin densities.15 Since their reported geometrical properties are similar to ours and since the same BLYP/double-z approach was applied, we attribute this difference to their more extended definition of the QM region, with cysteinates represented as [S–CH2–CH3] 1 instead of [S–CH3] 1. The Ni spin density in Table 2 is more similar to the value of 0.67 calculated with a more extensive model including the hydrogen-bonded protein residues, although the spin density on CysB2(Sg) is higher owing to the explicit inclusion of the histidine hydrogen bond in the QM region.18 Tables 1 and 2 list selected internal coordinates and Mulliken populations. The most noticeable difference between the structures of the isolated model and the MD snapshots is apparent in the compression of the bonds and angles between Ni and the bridging cysteinates due to the protein matrix, specifically, Ni–CysB1(Sg), Ni–CysB2(Sg), Ni–CysB1(Sg)–Fe, and Ni–CysB2(Sg)–Fe. This compression explains the increased Ni electron density (lower charge) in the isolated model relative to the MD snapshots: there is less electron sharing between Ni and sulfur atoms in the gas phase so the positively-charged Ni ion has more electron density. Analysis of the Mulliken spin densities reveals that only the bridging cysteinate CysB2(Sg) is appreciably influenced by hydrogen bonding (the effect of His on CysB2(Sg) is discussed in Table 1 Selected internal coordinates of the Ni–C active site calculated from the isolated cluster model, trajectory snapshots, and from the Desulfovibrio fructosovorans crystal structure
Coordinate
Isolated cluster model MD snapshots 1YRQ.pdb
Ni–Fe (Å) Ni–H (Å) Fe–H (Å) Ni–CysB1(Sg) (Å) Ni–CysB2(Sg) (Å) Ni–CysT1(Sg) (Å) Ni–CysT2(Sg) (Å) Fe–CysB1(Sg) (Å) Fe–CysB2(Sg) (Å) Ni–H –Fe (1) Ni–CysB1(Sg)–Fe (1) Ni–CysB2(Sg)–Fe (1)
2.64 1.60 1.73 2.34 2.44 2.26 2.29 2.38 2.34 105.1 78.9 88.4
2.64 0.05 1.59 0.07 1.75 0.07 2.25 0.06 2.38 0.07 2.28 0.08 2.23 0.07 2.35 0.06 2.37 0.07 104.4 3.2 70.0 1.8 67.7 1.9
24030 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
2.82 — — 2.29 2.54 2.15 2.21 2.29 2.28 — 76.1 71.4
Table 2 Mulliken populations calculated from the isolated cluster model and from trajectory snapshots. All quantities are elementary charge units (electrons)
Mulliken population Ni charge Fe charge CysB2(Sg) charge CysT2(Sg) charge H charge Ni spin Fe spin CysB2(Sg) spin CysT2(Sg) spin H spin
Isolated MD snapshots MD snapshots cluster model (protein excluded) (protein included) 1.043 0.758 0.444 0.358 0.552 0.639 0.017 0.204 0.062 0.010
1.142 0.534 0.438 0.205 0.599 0.634 0.043 0.233 0.074 0.013
0.142 0.043 0.082 0.082 0.099 0.039 0.021 0.032 0.029 0.020
1.150 0.514 0.642 0.250 0.572 0.646 0.018 0.207 0.121 0.019
0.143 0.044 0.118 0.086 0.101 0.040 0.014 0.032 0.037 0.018
detail in ref. 18). In contrast, the net electron populations (charges) are more strongly influenced by hydrogen bonding with outer coordination sphere, and we consider this quantity to be a proxy for the electron affinity [Ni(III) - Ni(II)] and proton affinity [CysT2(Sg) - CysT2(Sg-H)] aspects of Ni-R formation from Ni-C. The major effect is observed for CysB2(Sg), which is hydrogen bonded to His79; removal of the His–Cys interaction causes a +0.2e shift of the charge of CysB2(Sg). The terminal cysteinate CysT2(Sg) is also noticeably affected by hydrogen bonding to the backbone NH group of CysB2, which lowers the CysT2(Sg) charge by +0.05e. The Mulliken population on the bridging hydride is quite similar in all of the models studied here, being differences well smaller than the standard deviation of thermally averaged charges calculated from MD snapshots. Since the electron redistribution due to outer-sphere hydrogen bonding is primarily at Ni, CysB2(Sg) and CysT2(Sg), we applied multi-variable linear regression to assess the relative contributions of CysB2(Sg) and CysT2(Sg) to fluctuations in Ni electron densities. From this analysis, only CysT2(Sg) is associated with Ni charge fluctuations: q(Ni) = 0.05q(CysB2) 0.71q(CysT2). Charge exchange between Ni and CysT2(Sc) is driven by cysteinate rotation From the FMA analysis we have determined modes of motion modulating the Mulliken populations. Specifically we focused on the Ni and CysT2 charges because of their relevance to formation of the Ni-R state and the multi-variable regression analysis result. The model-building correlation coefficients RM are 0.96 and 0.92 and cross-correlation coefficients RC are 0.93 and 0.92 for the Ni and CysT2(Sg) charges, respectively, indicating a reasonable compromise between model-building and predictive power of the FMA modes and satisfactory correlation of this motion with the time-dependent property. Fig. 3 shows the Mulliken population FMA modes as vectors superimposed onto the average geometry of the inner coordination sphere. The most interesting feature of the Mulliken population-modulating modes is the nearly anti-parallel vectors corresponding to the Ni and CysT2(Sg) charges. The dihedral angles defining terminal cysteinate rotation [f1 = Fe–Ni–CysT1(Sg)–CysT1(Cb) and f2 = Fe–Ni–CysT2(Sg)–CysT2(Cb)] are the internal coordinates contributing most to the Ni and CysT2(Sg) charge fluctuations, with a
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
Sg atoms and the outer coordination sphere bridging cysteinate NH groups are the main elements modulating the dihedral f2 (and f1). This is evident form Fig. 4 (right panel) where a clear correlation between the hydrogen bond distance between CysB2(N) and CysT2(Sg) and the dihedral f2 is evident. Within one standard deviation of f2, the Ni charge decreases from 1.25e to 1.05e (a charge shift D = 0.2e) and the CysT2(Sg) charge increases from 0.31e to 0.22e (D = +0.09e), while the CysB2(N)–CysT2(Sg) hydrogen bond decreases from 3.40 to 3.17 Å (the distance between donor and acceptor atoms; the actual hydrogen bond distance is B1 Å shorter). The remaining Ni charge shifts to CysB2(Sg) (D = +0.04e), the hydride (D = +0.05e), and other atoms (0.02e).
Conclusions Fig. 3 Visualization of inner coordination sphere motion associated with increasing Ni charge (red vectors) and increasing CysT2(Sg) charge (black vectors) from functional mode analysis. See Fig. 1 caption for atom names and coloring scheme. The dihedral angle corresponding to CysT2 rotation is the greatest contributor to the charge distribution between Ni and CysT2(Sg).
small dihedral angle corresponding to a shorter distance to the nearest bridging cysteinate. To further investigate the relationship between terminal cysteinate rotation, hydrogen bonding and charge exchange between Ni and CysT2(Sg), we have plotted these charges and hydrogen bond lengths, bin-averaged over the dihedral angle f2 using Scott’s rule for the bin width42 (Fig. 4). The Mulliken charges show a clear sinusoidal dependence upon f2 with an increasing charge on CysT2(Sg) and a decreasing charge on Ni as f2 increases. The hydrogen bond between terminal cysteinate
The [NiFe] hydrogenase features several different catalytically active redox states (Scheme 1). The equilibrium between Ni-C and Ni-R states is believed to be the critical step for H2 cleavage or H2 bond formation. These two states contain Ni(III) and Ni(II) centers, respectively.21 Toward a better understanding of the structural, dynamical and electronic factors regulating the Ni-C - Ni-R equilibrium, we have carried out an analysis of how the protein environment and its thermal motion influence the electronic charge distribution of the Ni-C state. To this end, we performed hybrid DFT/MM calculations on an ensemble of active site structures with and without the protein environment around the NiFe catalytic core. Our results suggest a role for protein dynamics (and distortions of the active site geometry) in modulating the electron density on the nickel and terminal cysteinate ligand. We have found that the CysB2 backbone amide group induces the rotation of the terminal CysT2(Sg), identified by the dihedral angle f2 in Fig. 4. This rotational motion causes
Fig. 4 Left: Mulliken charge distribution between Ni and CysT2(Sg) follows the rotation of CysT2, defined by the dihedral angle f2 = Fe–Ni–CysT2(Sg)–CysT2(Cb). Dotted line is the average f2 and the grey shaded area is the standard deviation from the average. Circles: charges bin-averaged over f2; error bars: standard deviations of bin-averaged charges; lines: cosine fit (correlation coefficient R2 E 1.0 for Ni and 0.9763 for CysT2(Sg), respectively). Right: dihedral angle f2 is driven by the hydrogen bond between the backbone of the neighboring bridging and terminal cysteinates. Circles: hydrogen bond distance bin-averaged over f2; error bars: standard deviations of bin-averaged distances; line: cosine fit (R2 = 0.98). The dihedral f2 is 34.31 in the initial crystal structure.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24031
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PCCP
a flow of charge between the Ni center and CysT2(Sg) atom. As the CysT2 side chain moves closer to CysB2 there is a decrease of electron density on the Ni center and an increase of electron density on CysT2(Sg) atom. Consequently, the electron affinity of the Ni center increases, making it thermodynamically favorable to reduce Ni(III) to Ni(II), while simultaneously increasing the proton affinity of CysT2(Sg), thereby stabilizing the H2-cleaved Ni-R which itself is iso-energetic with the Ni-bound dihydrogen structure.21 Furthermore, this cysteinate rotation exposes the nickel ion to the gas channel, thereby facilitating H2 egress.43,44 This finding suggests that this motion prepares the inner coordination sphere of the Ni-C state for conversion to Ni-R state. The scenario described here is consistent with the observed relationship between the dihedral angle f2 and the overall H2 production rate, in which a smaller f2 value coincides with increased H2 production.17 Our results support previous studies,1,18 which emphasized that the polypeptide backbone (and 4 cysteine ligands) plays a major role in controlling the electronic properties and catalytic function of the [NiFe] hydrogenase, and that the sulfur atoms ‘‘affect the spin and charge density delocalization to modify the electronic properties of the metal site’’.1,18 Further, these previous studies emphasize that the nickel is the redox active metal, whose electronic configuration is sensitive to the geometry and charge state of the cysteine sulfurs. In summary, our study provides support for the role of protein dynamics in controlling the electron density, and presumably the reactivity, at the active site nickel center in the [NiFe] hydrogenase. Additionally, in future studies, the methods and findings reported here could be expanded to explicitly include the Ni-SIa and Ni-R catalytic intermediates.
Acknowledgements This work was funded by the US DOE Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences and Biosciences. Computer resources were provided by the Environmental Molecular Sciences Laboratory at the Pacific Northwest National Laboratory. The authors thank Aaron Appel and Marat Valiev of PNNL for helpful discussions.
References 1 M. E. Pandelia, H. Ogata and W. Lubitz, ChemPhysChem, 2010, 11, 1127–1140. 2 B. Friedrich, J. Fritsch and O. Lenz, Curr. Opin. Biotechnol., 2011, 22, 358–364. 3 A. Volbeda, L. Martin, C. Cavazza, M. Matho, B. W. Faber, W. Roseboom, S. P. J. Albracht, E. Garcin, M. Rousset and J. C. Fontecilla-Camps, J. Biol. Inorg. Chem., 2005, 10, 239–249. 4 D. L. DuBois, Inorg. Chem., 2014, 53, 3935–3960. 5 S. B. Dementin, F. Leroux, L. Cournac, A. L. de Lacey, ´ger, B. Burlat, N. Martinez, S. Champ, A. Volbeda, C. Le L. Martin, O. Sanganas, M. Haumann, V. M. Fernandez, B. Guigliarelli, J. C. Fontecilla-Camps and M. Rousset, J. Am. Chem. Soc., 2009, 131, 10156–10164.
24032 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
Paper
6 E. Garcin, X. Vernede, E. C. Hatchikian, A. Volbeda, M. Frey and J. C. Fontecilla-Camps, Structure, 1999, 7, 557–566. 7 Y. Higuchi, H. Ogata, K. Miki, N. Yasuoka and T. Yagi, Structure, 1999, 7, 549–556. 8 Y. Higuchi, T. Yagi and N. Yasuoka, Structure, 1997, 5, 1671–1680. 9 F. Leroux, S. Dementin, B. Burlatt, L. Cournac, A. Volbeda, S. Champ, L. Martin, B. Guigliarelli, P. Bertrand, J. FontecillaCamps, M. Rousset and C. Leger, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 11188–11193. 10 M. C. Marques, R. Coelho, A. L. De Lacey, I. A. C. Pereira and P. M. Matias, J. Mol. Biol., 2010, 396, 893–907. 11 P. M. Matias, C. M. Soares, L. M. Saraiva, R. Coelho, J. Morais, J. Le Gall and M. A. Carrondo, J. Biol. Inorg. Chem., 2001, 6, 63–81. 12 H. Ogata, S. Hirota, A. Nakahara, H. Komori, N. Shibata, T. Kato, K. Kano and Y. Higuchi, Structure, 2005, 13, 1635–1642. 13 H. Ogata, P. Kellers and W. Lubitz, J. Mol. Biol., 2010, 402, 428–444. 14 A. Volbeda, M.-H. Charon, C. Piras, E. C. Hatchikian, M. Frey and J. C. Fontecilla-Camps, Nature, 1995, 373, 580–587. 15 M. Stein and W. Lubitz, Phys. Chem. Chem. Phys., 2001, 3, 2668–2675. 16 R. Cammack, Nature, 1999, 397, 214–215. 17 D. M. A. Smith, Y. Xiong, T. P. Straatsma, K. M. Rosso and T. C. Squier, J. Chem. Theory Comput., 2012, 8, 2103–2114. 18 M. Kampa, W. Lubitz, M. van Gastel and F. Neese, J. Biol. Inorg. Chem., 2012, 17, 1269–1281. 19 M. P. Shaver, L. E. N. Allan, H. S. Rzepa and V. C. Gibson, Angew. Chem., Int. Ed., 2006, 45, 1241–1244. 20 R. L. Yson, J. L. Gilgor, B. A. Guberman and S. A. Varganov, Chem. Phys. Lett., 2013, 577, 138–141. ¨mer, M. Kampa, W. Lubitz, M. van Gastel and 21 T. Kra F. Neese, ChemBioChem, 2013, 14, 1898–1905. 22 P. E. M. Siegbahn, J. W. Tye and M. B. Hall, Chem. Rev., 2007, 107, 4414–4435. 23 M. Bruschi, G. Zampella, P. Fantucci and L. De Gioia, Coord. Chem. Rev., 2005, 249, 1620–1640. 24 H. Li, A. D. Robertson and J. H. Jensen, Proteins: Struct., Funct., Bioinf., 2005, 61, 704–721. 25 Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. M. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. M. Wang and P. Kollman, J. Comput. Chem., 2003, 24, 1999–2012. 26 M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys., 2000, 112, 8910–8922. 27 A. L. Cheng and K. M. Merz, J. Phys. Chem., 1996, 100, 1927–1937. 28 G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101. 29 E. Lindahl, B. Hess and D. van der Spoel, J. Mol. Model., 2001, 7, 306–317. 30 A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
31 C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789. 32 M. J. Frisch, J. A. Pople and J. S. Binkley, J. Chem. Phys., 1984, 80, 3265–3269. 33 A. Schafer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577. 34 C. H. Chang and K. Kim, J. Chem. Theory Comput., 2009, 5, 1137–1145. 35 M. Stein and W. Lubitz, Curr. Opin. Chem. Biol., 2002, 6, 243–249. 36 M. G. Delcey, K. Pierloot, Q. M. Phung, S. Vancoillie, R. Lindh and U. Ryde, Phys. Chem. Chem. Phys., 2014, 16, 7927–7938. 37 L. M. Lawson Daku, F. Aquilante, T. W. Robinson and A. Hauser, J. Chem. Theory Comput., 2012, 8, 4216–4231.
This journal is © the Owner Societies 2014
PCCP
38 C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280. 39 M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. Van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus and W. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489. 40 K. C. Gross, P. G. Seybold and C. M. Hadad, Int. J. Quantum Chem., 2002, 90, 445–458. 41 J. S. Hub and B. L. de Groot, PLoS Comput. Biol., 2009, 5, e1000480, DOI: 10.1371/journal.pcbi.1000480. 42 D. W. Scott, Biometrika, 1971, 66, 605–610. 43 Y. Montet, P. Amara, A. Volbeda, X. Vernede, E. C. Hatchikian, M. J. Field, M. Frey and J. C. FontecillaCamps, Nat. Struct. Biol., 1997, 4, 523–526. 44 V. H. Teixeira, A. M. Baptista and C. M. Soares, Biophys. J., 2006, 91, 2035–2045.
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24033
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 24026
View Journal | View Issue
Modulation of active site electronic structure by the protein matrix to control [NiFe] hydrogenase reactivity† Dayle M. A. Smith,*a Simone Raugeia and Thomas C. Squierb Control of the reactivity of the nickel center of the [NiFe] hydrogenase and other metalloproteins commonly involves outer coordination sphere ligands that act to modify the geometry and physical properties of the active site metal centers. We carried out a combined set of classical molecular dynamics and quantum/classical mechanics calculations to provide quantitative estimates of how dynamic fluctuations of the active site within the protein matrix modulate the electronic structure at the catalytic center. Specifically we focused on the dynamics of the inner and outer coordination spheres
Received 6th August 2014, Accepted 30th September 2014
of the cysteinate-bound Ni–Fe cluster in the catalytically active Ni-C state. There are correlated
DOI: 10.1039/c4cp03518f
the electron affinity at the active site and the proton affinity of a terminal cysteinate. On the basis of
movements of the cysteinate ligands and the surrounding hydrogen-bonding network, which modulate these findings, we hypothesize a coupling between protein dynamics and electron and proton transfer
www.rsc.org/pccp
reactions critical to dihydrogen production.
Introduction The increased demand for energy from intermittent renewable sources such as solar and wind power requires the ability to store the excess energy produced and retrieve it on demand. In this context, electrocatalysts for rapid and efficient interconversion between electrical energy and chemical energy will be of critical importance. H2 production and oxidation represent central reactions that can be used for energy storage. Hydrogenase enzymes have been widely studied in order to gain insight into H2 activation by metal centers, and coupled electron and proton transfer reactions. This family of enzymes catalyzes the reversible oxidation of dihydrogen at a bimetallic active site ([FeFe] and [NiFe] hydrogenases) or a monometallic ([Fe]-only hydrogenases) active site. [FeFe]-hydrogenase enzymes exhibit, in general, superior rates of hydrogen evolution, but are irreversibly inactivated by the presence of even trace amounts of oxygen. In comparison, the active site structure of [NiFe]hydrogenase enzymes are reversibly inactivated by molecular oxygen, forming catalytically inactive ready (Ni-B) and unready a
Pacific Northwest National Laboratory, P.O. Box 999, MSIN J4-33, Richland, Washington 99352, USA. E-mail: [email protected] b Western University of Health Sciences, 200 Mullins Drive, Lebanon, Oregon 97335, USA † Electronic supplementary information (ESI) available: [NiFe] per-residue rootmean-squared fluctuations, solvent-accessible surface area and secondary structure representations are tabulated in Table S1. [NiFe] trajectories can be obtained by contacting the corresponding author. See DOI: 10.1039/c4cp03518f
24026 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
(Ni-A) states that respectively possess a bridging hydroxide and peroxide between the metal centers.1 The resistance of the active site structure of [NiFe]-hydrogenases to molecular oxygen,2 such that reductive reactivation occurs following inactivation, provides a central design principle important to the creation of robust molecular catalysts. As the active site structure is maintained between all [NiFe]-hydrogenases, the commonly studied enzyme from Desulfovibrio fructosovorans (1YRQ.pdb)3 provides a solid basis to understand the design principles of the hydrogenase protein that allows catalysis at substantially lower overpotentials than observed for currently available molecular catalysts.4 [NiFe] hydrogenase enzymes (Fig. 1a) are heterodimers, consisting of a 34 kDa small subunit and 61 kDa large subunit. The small subunit contains three iron–sulfur (Fe–S) clusters forming an electron transfer chain that shuttles electrons from cellular redox donors (e.g., membrane bound cytochromes) to the bimetallic active site in the large subunit. The active site consists of a nickel center, an iron center, three diatomic ligands (CO and CN ) bound to the iron, two cysteinates (deprotonated cysteines) bridging the two metal atoms, and two cysteinates bound only to the nickel atom (terminal cysteinates). Fig. 1b shows the active site structure including hydrogen bonded groups as inferred from the X-ray crystal structure from Desulfovibrio fructosovorans.3 The inner coordination sphere is defined as the active site and cysteinate atoms Sg and Cb. Cys75 and Cys546 bind to both Fe and Ni atoms (bridging cysteinates CysB1 and CysB2) and Cys72 and Cys543 bind only to
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
Fig. 1 (a) Structure of the ready, oxidized [NiFe] hydrogenase from Desulfovibrio fructosovorans PDB ID: 1YRQ. Protein residues are shown as ribbons (cyan: small subunit; orange: large subunit), metal clusters are shown as van der Waals spheres, and hydrogen bonded residues are shown as sticks. (b) [NiFe] active site, bound cysteinates and hydrogen bonded protein residues. The inner sphere (QM) region is drawn with sticks; residues hydrogen bonded to the NiFe core are drawn as lines; atoms hydrogen bonded to the QM region are spheres.
Ni (terminal cysteinates CysT1 and CysT2). Ser499 hydrogen bonds to one cyanide ligand (CN1), Arg476 hydrogen bonds to the other cyanide (CN2), and His79 hydrogen bonds to the sulfur atom of bridging cysteinate CysB2. There are also hydrogen bonds between bridging cysteinate backbone NH groups and their nearest terminal cysteinate sulfur atoms, which have the potential to drive rotation of cysteine sidechains, as described in Results. Given the high degrees of conservation between the active site structures of all [NiFe] hydrogenases,3,5–16 differences in catalytic turnover rates for H2 production suggest an important role for the surrounding protein matrix in modulating catalysis. In this respect, Pandelia and coworkers1 proved that the nickel is the redox active metal, and its electronic configuration is sensitive to the geometry and charge state of the cysteine sulfurs, which is modulated by the protein backbone by ‘‘affect[ing] the spin and charge density delocalization to modify the electronic properties of the metal site’’. These observations are consistent with empirical correlations that indicate the geometry around the nickel [i.e., the dihedral angle Fe–Ni–CysT2(Sg)–CysT2(Cb), Fig. 1] are responsible for differences in rates of catalysis, such that a smaller dihedral angle corresponds to a faster rate of hydrogen production.17 Indeed, recent density functional theory (DFT) calculations suggest that the Ni spin density and nickel–sulfur bond order depend on the terminal cysteinate orientation.18 Further, rotation of CysT2 modulates the relative stability of the low- and high-spin Ni electron configuration,19 consistent with suggestions that a terminal cysteinate (i.e., CysT2) acts as the proton acceptor20 and that the surrounding
This journal is © the Owner Societies 2014
protein matrix may control dihedral movements to gate proton transfer reactions within the catalytic core and modulate rates of H2 production. As part of the reaction mechanism, which is likely to involve redox-dependent changes in geometry and flexibility of the atoms around the active site metals, a bridging ligand can bind between the Fe and Ni metal centers. The most important catalytically active redox states of [NiFe] hydrogenase are Ni-SIa, Ni-C, and Ni-R (Scheme 1). The Ni-SIa state contains Fe(II)–Ni(II) and a loosely-bound water molecule on the Ni center; the Ni-C state contains Fe(II)–Ni(III) and a bridging hydride between the two metals; finally the Ni-R state contains Fe(II)–Ni(II) and features two isoenergetic, doubly-protonated structures. In one Ni-R structure, H2 is bound to the Ni ion; the other structure contains a bridging hydride and a protonated terminal cysteinate.21 The equilibrium between Ni-C and Ni-R is therefore a critical step for H2 cleavage or H2 production. Our goal is to elucidate the role of protein dynamics (and distortions of the active site geometry) in modulating the electronic properties of the nickel and terminal cysteinate ligand that receives a proton during the catalytic splitting cycle. To explore how protein fluctuations modulate the electron density at the active site in [NiFe] hydrogenases, we have carried out a combined set of classical molecular dynamics (MD) and hybrid quantum/classical mechanics (QM/MM) calculations on the catalytically active Ni-C state (Scheme 1). Our approach is complementary to prior approaches that used DFT calculations to model the active site. DFT calculations of gas-phase active site models have been instrumental in
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24027
View Article Online
PCCP
Paper
for the coordination sphere of the NiFe core derived from the classical force field for electronic structure calculations.
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Construction of Ni-C [NiFe] hydrogenase in aqueous solution
Scheme 1 Active catalytic cycle between Ni-SIa, Ni-C and Ni-R states of [NiFe] hydrogenase, including the isoenergetic, H2 splitting Ni-R conformations.21
identifying reaction intermediates alongside experimental spectroscopic measurements.1,15,22,23 More recent work has extended the model to include hydrogen-bonded protein residues and shed light on the effects of hydrogen bonding on the electron distribution of the catalytically-important active site atoms.18 Our approach is unique in that we used a customized force field17 and molecular dynamics simulations to collect an in situ ensemble of active site structures for QM/MM electronic structure calculations. We find that hydrogen bonding between bridging and terminal cysteine ligands modulate the dynamics of the cysteinate ligands and, in turn, the electron affinity of the Ni center and the proton affinity of one of the cysteinates. On the basis of these findings, we suggest a dynamic coupling between the protein matrix and the catalytic core that is critical in mediating electron density changes associated with catalysis (i.e., driving the Ni-C to Ni-R reaction necessary for hydrogen production).
Methods Overall approach We performed a sub-microsecond (0.1 ms) MD simulation of the Ni-C state of [NiFe] hydrogenase to sample long time-scale protein fluctuations in aqueous solution by using an accurate force field recently parameterized by us from DFT calculations.17 From this simulation, a statistically representative number of configurations were extracted for subsequent electronic structure analysis. For each configuration an analysis of the electron distribution of the catalytic core was performed using hybrid density functional theory (DFT)/molecular mechanics calculations. The DFT region includes the same atoms as in the modelcluster for the force field parameterization (Fig. 1b), whereas for the surrounding protein environment was treated at the force field level. This choice was adopted to avoid possible artifacts introduced by inconsistencies due to use of structures
24028 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
The initial model of the Ni-C [NiFe] hydrogenase was built starting from a X-ray structure Desulfovibrio fructosovorans at 1.83 Å resolution (protein data bank entry 1YRQ). This structure is in the oxidized (hydroxide-bridging) Ni-B state, and the bridging ligand was replaced by hydrogen to build the hydride-bound Ni-C active site. Histidine protonation states were determined using PropKa program24 and visual inspection. This analysis showed that His184 of the small subunit and His79 in the large subunit are protonated at Ne; His61 and His192 of the small subunit and His204, His210, and His549 of the large subunit are protonated at both Nd and Ne; the remaining histidine residues are protonated at Nd. The protein was inserted into a cubic box with minimum 12 Å distance between the protein and the side of the box, then the box was filled with water molecules. Sodium ions were added for charge neutrality. Molecular dynamics simulation Prerequisite to the classical MD simulations is an empirical force field that reliably describes the geometry and motion of the active site within the protein matrix. The force field parameters we derived and published17 reproduce the geometries and frequencies of gas-phase model clusters of Ni-B and Ni-C active sites (and iron– sulfur clusters) with ligated cysteinates modeled by methionates, and also reasonably reproduce the hydrogen bond distances in the initial protein structure. Using the AMBER0325 force field and TIP3P26 water model with the embedded metal cluster parameters for the Ni-C active site and reduced Fe–S clusters, the initial structure was relaxed employing conjugate gradient energy minimization to a tolerance of 50 kJ mol 1 Å 1 and was then used as a starting point for the MD simulation. Equilibration was done in 250 ps long consecutive runs, in which the coordinates and velocities of each equilibration step were used to initialize each subsequent step. In the first equilibration step, non-hydrogen atoms were restrained using harmonic position restraints with a force constant of 10 kJ mol 1 Å 2, and the temperature was set at 100 K. Next, the temperature was increased in 50 K increments to 300 K; then, restraint forces were reduced to zero in 2 kJ mol 1 Å 2 decrements. Finally, the system was left free to evolve for 100 ns. All ´–Hoover thermoof the simulations were performed using a Nose stat27 and a modified Berendsen barostat.28 A 1 fs time step was used without imposing constraints on the chemical bonds. Electrostatic interactions were calculated using the Particle Mesh Ewald method with a cutoff of 12 Å. Analysis of the root mean square displacement (RMSD) reveals that equilibration is reached in about 20 ns (Fig. 2); therefore average quantities were calculated from 20 ns and 100 ns. All classical mechanics simulations and the subsequent analysis of the trajectory were done using the doubleprecision version of GROMACS 4.0 and related analysis programs.29 Hybrid density functional theory/molecular mechanics calculations Hybrid DFT/MM calculations were carried out to describe the catalytic core on a statistically relevant number of configurations
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
configurations taken from the Ni-C [NiFe] hydrogenase trajectory were used for single-point DFT calculations. A total of 1600 configurations 50 ps apart were extracted from the trajectory between 20 and 100 ns. One set of calculations was run without the protein environment, and another set was run with the full protein environment included using a QM/MM scheme, in order to separate contributions from inner-sphere and outer-sphere fluctuations. The electrostatic coupling between the QM region and MM region was described in ref. 38. All of the QM calculations were done using the NWChem DFT module via the QMMM module using default parameters.39 Dynamical charge density analysis
Fig. 2 Black line, top: root-mean-squared deviation of [NiFe] hydrogenase Ca atoms from the 100 ns MD trajectory, relative to the initial crystallographic structure. Grey line, bottom: root-mean-squared deviation of [NiFe] active site atoms and hydrogen bonded Arg476, Ser499 and His79 residues with respect to the initial structure.
extracted from the force field-based MD simulations. All DFT calculations used the BLYP exchange and correlation functional30,31 and a double-z basis set with diffuse functions.32,33 BLYP/double-z has been used previously to successfully describe structures and experimental magnetic properties of the [NiFe] active site.15 Furthermore, the choice of functional and basis set was dictated by consistency considerations with DFT setup used for the [NiFe] force field parameterization17 based on model active site geometry optimization and harmonic frequency calculations (also the same approach used for the DFT-derived [FeFe] hydrogenase force field34). It is important to remark that our previous calculations showed that Mulliken populations and vibrational frequencies are very similar from the double- and triple-z basis sets,17 and different functionals have been shown to produce similar structural properties in [NiFe] active site models.35 As with other transition metal complexes, the high-versus-low spin energy difference of the nickel ion is sensitive to the choice of functional, relative to more accurate post-Hartree–Fock results;36,37 however, the inadequacy of BLYP in this regard is a separate issue from results presented here since our approach explores the potential energy surface of the more stable low-spin state. Since we are using a combined MD/QMMM approach to model the dependence of net electron densities on active site fluctuations within the protein, consistency with the classical force field is paramount. The QM region was selected to include the active site, where the cysteinates coordinated to the metal centers were treated as methionates, again for consistency with the cluster DFT calculations employed for the force field parameterization. [NiFe] hydrogenase
This journal is © the Owner Societies 2014
Mulliken charges and spin densities were used as descriptors of the electron density as a proxy for chemical reactivity (Scheme 1). Specifically, Ni electron deficiency (positive atomic Mulliken charge on the Ni ion) was used to describe the relationship between active site dynamics and Ni(III) - Ni(II), based on the identification of nickel as the redox-active metal,1 and excess CysT2(Sg) electron density was used to describe CysT2(Sg) CysT2(Sg–H), based on previous quantum chemical calculations that found a correlation between Mulliken population and pKa.40 Covariance analysis and functional mode analysis method (FMA)41 were used to calculate the modes of motion that best describe fluctuations of selected Mulliken populations. FMA calculates a single mode of motion representative of a timedependent observable property, in this case electron densities. For this analysis we used all covariance eigenvectors and snapshots between 20 and 80 ns for model building and 80–100 ns for cross-validation. The model building and cross-validation correlation coefficients (RM and RC) are the assessment quantities used to validate that the resulting motion is maximally correlated to the quantity of interest.
Results and discussion Flexibility of the protein backbone and residues near the catalytic core The RMSD of the protein backbone with respect to the starting crystal structure (Fig. 2) indicates that the protein is relatively rigid, as discussed in previous reports, with fluctuation of about 0.1 Å around an equilibrium average value 1.5 Å. Analysis of the root mean square fluctuation of each single residue reveals that the flexible protein elements are residues near the solvent accessible surface and residues in unstructured regions such as chain termini (per-residue RMSF, solvent-accessible surface area and secondary structure classifications are tabulated in the ESI†). The fluctuations apparent in the RMSD of the inner and outer coordination sphere atoms (Fig. 2) are primarily due to His79 and the terminal cysteinates, which are known to be important for regulating spin density.18 Comparison of Mulliken populations and contributions from the inner and outer coordination spheres Before discussing the relationships between fluctuations and electron distributions, it is useful to show how the geometries
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24029
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PCCP
Paper
and Mulliken populations compare using the isolated cluster model and the same atom selection extracted from the MD snapshots for purposes of assessment and also to see how fluctuations of the active site within the protein matrix influence the electronic distribution at the catalytic center. The Ni, CysB2(Sg) and CysT2(Sg) atoms are the most important atoms for consideration because of their known sensitivity to hydrogen bonding [CysB2(Sg)] and direct participation in conversion to the Ni-R state when dihydrogen splitting occurs [Ni and CysT2(Sg)].21 Therefore, our discussion will focus mainly on these centers and, for sake of completeness, on the hydride and the Fe center. In the isolated cluster model, the bulk of the spin density is on the Ni ion (0.64e). Additional spin density is on one of the bridging cysteinate sulfur atom CysB2(Sg) (0.20e) and the terminal cysteinate sulfur atom CysT2(Sg) (0.06e). A very small amount of spin density is present also on the Fe atom (0.02e), which is in the low-spin Fe(II) oxidation state. A previous study, which used the same DFT functional and a similar basis set, reported 0.51e and 0.29e for the Ni and CysB2(Sg) spin densities.15 Since their reported geometrical properties are similar to ours and since the same BLYP/double-z approach was applied, we attribute this difference to their more extended definition of the QM region, with cysteinates represented as [S–CH2–CH3] 1 instead of [S–CH3] 1. The Ni spin density in Table 2 is more similar to the value of 0.67 calculated with a more extensive model including the hydrogen-bonded protein residues, although the spin density on CysB2(Sg) is higher owing to the explicit inclusion of the histidine hydrogen bond in the QM region.18 Tables 1 and 2 list selected internal coordinates and Mulliken populations. The most noticeable difference between the structures of the isolated model and the MD snapshots is apparent in the compression of the bonds and angles between Ni and the bridging cysteinates due to the protein matrix, specifically, Ni–CysB1(Sg), Ni–CysB2(Sg), Ni–CysB1(Sg)–Fe, and Ni–CysB2(Sg)–Fe. This compression explains the increased Ni electron density (lower charge) in the isolated model relative to the MD snapshots: there is less electron sharing between Ni and sulfur atoms in the gas phase so the positively-charged Ni ion has more electron density. Analysis of the Mulliken spin densities reveals that only the bridging cysteinate CysB2(Sg) is appreciably influenced by hydrogen bonding (the effect of His on CysB2(Sg) is discussed in Table 1 Selected internal coordinates of the Ni–C active site calculated from the isolated cluster model, trajectory snapshots, and from the Desulfovibrio fructosovorans crystal structure
Coordinate
Isolated cluster model MD snapshots 1YRQ.pdb
Ni–Fe (Å) Ni–H (Å) Fe–H (Å) Ni–CysB1(Sg) (Å) Ni–CysB2(Sg) (Å) Ni–CysT1(Sg) (Å) Ni–CysT2(Sg) (Å) Fe–CysB1(Sg) (Å) Fe–CysB2(Sg) (Å) Ni–H –Fe (1) Ni–CysB1(Sg)–Fe (1) Ni–CysB2(Sg)–Fe (1)
2.64 1.60 1.73 2.34 2.44 2.26 2.29 2.38 2.34 105.1 78.9 88.4
2.64 0.05 1.59 0.07 1.75 0.07 2.25 0.06 2.38 0.07 2.28 0.08 2.23 0.07 2.35 0.06 2.37 0.07 104.4 3.2 70.0 1.8 67.7 1.9
24030 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
2.82 — — 2.29 2.54 2.15 2.21 2.29 2.28 — 76.1 71.4
Table 2 Mulliken populations calculated from the isolated cluster model and from trajectory snapshots. All quantities are elementary charge units (electrons)
Mulliken population Ni charge Fe charge CysB2(Sg) charge CysT2(Sg) charge H charge Ni spin Fe spin CysB2(Sg) spin CysT2(Sg) spin H spin
Isolated MD snapshots MD snapshots cluster model (protein excluded) (protein included) 1.043 0.758 0.444 0.358 0.552 0.639 0.017 0.204 0.062 0.010
1.142 0.534 0.438 0.205 0.599 0.634 0.043 0.233 0.074 0.013
0.142 0.043 0.082 0.082 0.099 0.039 0.021 0.032 0.029 0.020
1.150 0.514 0.642 0.250 0.572 0.646 0.018 0.207 0.121 0.019
0.143 0.044 0.118 0.086 0.101 0.040 0.014 0.032 0.037 0.018
detail in ref. 18). In contrast, the net electron populations (charges) are more strongly influenced by hydrogen bonding with outer coordination sphere, and we consider this quantity to be a proxy for the electron affinity [Ni(III) - Ni(II)] and proton affinity [CysT2(Sg) - CysT2(Sg-H)] aspects of Ni-R formation from Ni-C. The major effect is observed for CysB2(Sg), which is hydrogen bonded to His79; removal of the His–Cys interaction causes a +0.2e shift of the charge of CysB2(Sg). The terminal cysteinate CysT2(Sg) is also noticeably affected by hydrogen bonding to the backbone NH group of CysB2, which lowers the CysT2(Sg) charge by +0.05e. The Mulliken population on the bridging hydride is quite similar in all of the models studied here, being differences well smaller than the standard deviation of thermally averaged charges calculated from MD snapshots. Since the electron redistribution due to outer-sphere hydrogen bonding is primarily at Ni, CysB2(Sg) and CysT2(Sg), we applied multi-variable linear regression to assess the relative contributions of CysB2(Sg) and CysT2(Sg) to fluctuations in Ni electron densities. From this analysis, only CysT2(Sg) is associated with Ni charge fluctuations: q(Ni) = 0.05q(CysB2) 0.71q(CysT2). Charge exchange between Ni and CysT2(Sc) is driven by cysteinate rotation From the FMA analysis we have determined modes of motion modulating the Mulliken populations. Specifically we focused on the Ni and CysT2 charges because of their relevance to formation of the Ni-R state and the multi-variable regression analysis result. The model-building correlation coefficients RM are 0.96 and 0.92 and cross-correlation coefficients RC are 0.93 and 0.92 for the Ni and CysT2(Sg) charges, respectively, indicating a reasonable compromise between model-building and predictive power of the FMA modes and satisfactory correlation of this motion with the time-dependent property. Fig. 3 shows the Mulliken population FMA modes as vectors superimposed onto the average geometry of the inner coordination sphere. The most interesting feature of the Mulliken population-modulating modes is the nearly anti-parallel vectors corresponding to the Ni and CysT2(Sg) charges. The dihedral angles defining terminal cysteinate rotation [f1 = Fe–Ni–CysT1(Sg)–CysT1(Cb) and f2 = Fe–Ni–CysT2(Sg)–CysT2(Cb)] are the internal coordinates contributing most to the Ni and CysT2(Sg) charge fluctuations, with a
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
PCCP
Sg atoms and the outer coordination sphere bridging cysteinate NH groups are the main elements modulating the dihedral f2 (and f1). This is evident form Fig. 4 (right panel) where a clear correlation between the hydrogen bond distance between CysB2(N) and CysT2(Sg) and the dihedral f2 is evident. Within one standard deviation of f2, the Ni charge decreases from 1.25e to 1.05e (a charge shift D = 0.2e) and the CysT2(Sg) charge increases from 0.31e to 0.22e (D = +0.09e), while the CysB2(N)–CysT2(Sg) hydrogen bond decreases from 3.40 to 3.17 Å (the distance between donor and acceptor atoms; the actual hydrogen bond distance is B1 Å shorter). The remaining Ni charge shifts to CysB2(Sg) (D = +0.04e), the hydride (D = +0.05e), and other atoms (0.02e).
Conclusions Fig. 3 Visualization of inner coordination sphere motion associated with increasing Ni charge (red vectors) and increasing CysT2(Sg) charge (black vectors) from functional mode analysis. See Fig. 1 caption for atom names and coloring scheme. The dihedral angle corresponding to CysT2 rotation is the greatest contributor to the charge distribution between Ni and CysT2(Sg).
small dihedral angle corresponding to a shorter distance to the nearest bridging cysteinate. To further investigate the relationship between terminal cysteinate rotation, hydrogen bonding and charge exchange between Ni and CysT2(Sg), we have plotted these charges and hydrogen bond lengths, bin-averaged over the dihedral angle f2 using Scott’s rule for the bin width42 (Fig. 4). The Mulliken charges show a clear sinusoidal dependence upon f2 with an increasing charge on CysT2(Sg) and a decreasing charge on Ni as f2 increases. The hydrogen bond between terminal cysteinate
The [NiFe] hydrogenase features several different catalytically active redox states (Scheme 1). The equilibrium between Ni-C and Ni-R states is believed to be the critical step for H2 cleavage or H2 bond formation. These two states contain Ni(III) and Ni(II) centers, respectively.21 Toward a better understanding of the structural, dynamical and electronic factors regulating the Ni-C - Ni-R equilibrium, we have carried out an analysis of how the protein environment and its thermal motion influence the electronic charge distribution of the Ni-C state. To this end, we performed hybrid DFT/MM calculations on an ensemble of active site structures with and without the protein environment around the NiFe catalytic core. Our results suggest a role for protein dynamics (and distortions of the active site geometry) in modulating the electron density on the nickel and terminal cysteinate ligand. We have found that the CysB2 backbone amide group induces the rotation of the terminal CysT2(Sg), identified by the dihedral angle f2 in Fig. 4. This rotational motion causes
Fig. 4 Left: Mulliken charge distribution between Ni and CysT2(Sg) follows the rotation of CysT2, defined by the dihedral angle f2 = Fe–Ni–CysT2(Sg)–CysT2(Cb). Dotted line is the average f2 and the grey shaded area is the standard deviation from the average. Circles: charges bin-averaged over f2; error bars: standard deviations of bin-averaged charges; lines: cosine fit (correlation coefficient R2 E 1.0 for Ni and 0.9763 for CysT2(Sg), respectively). Right: dihedral angle f2 is driven by the hydrogen bond between the backbone of the neighboring bridging and terminal cysteinates. Circles: hydrogen bond distance bin-averaged over f2; error bars: standard deviations of bin-averaged distances; line: cosine fit (R2 = 0.98). The dihedral f2 is 34.31 in the initial crystal structure.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24031
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
PCCP
a flow of charge between the Ni center and CysT2(Sg) atom. As the CysT2 side chain moves closer to CysB2 there is a decrease of electron density on the Ni center and an increase of electron density on CysT2(Sg) atom. Consequently, the electron affinity of the Ni center increases, making it thermodynamically favorable to reduce Ni(III) to Ni(II), while simultaneously increasing the proton affinity of CysT2(Sg), thereby stabilizing the H2-cleaved Ni-R which itself is iso-energetic with the Ni-bound dihydrogen structure.21 Furthermore, this cysteinate rotation exposes the nickel ion to the gas channel, thereby facilitating H2 egress.43,44 This finding suggests that this motion prepares the inner coordination sphere of the Ni-C state for conversion to Ni-R state. The scenario described here is consistent with the observed relationship between the dihedral angle f2 and the overall H2 production rate, in which a smaller f2 value coincides with increased H2 production.17 Our results support previous studies,1,18 which emphasized that the polypeptide backbone (and 4 cysteine ligands) plays a major role in controlling the electronic properties and catalytic function of the [NiFe] hydrogenase, and that the sulfur atoms ‘‘affect the spin and charge density delocalization to modify the electronic properties of the metal site’’.1,18 Further, these previous studies emphasize that the nickel is the redox active metal, whose electronic configuration is sensitive to the geometry and charge state of the cysteine sulfurs. In summary, our study provides support for the role of protein dynamics in controlling the electron density, and presumably the reactivity, at the active site nickel center in the [NiFe] hydrogenase. Additionally, in future studies, the methods and findings reported here could be expanded to explicitly include the Ni-SIa and Ni-R catalytic intermediates.
Acknowledgements This work was funded by the US DOE Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences and Biosciences. Computer resources were provided by the Environmental Molecular Sciences Laboratory at the Pacific Northwest National Laboratory. The authors thank Aaron Appel and Marat Valiev of PNNL for helpful discussions.
References 1 M. E. Pandelia, H. Ogata and W. Lubitz, ChemPhysChem, 2010, 11, 1127–1140. 2 B. Friedrich, J. Fritsch and O. Lenz, Curr. Opin. Biotechnol., 2011, 22, 358–364. 3 A. Volbeda, L. Martin, C. Cavazza, M. Matho, B. W. Faber, W. Roseboom, S. P. J. Albracht, E. Garcin, M. Rousset and J. C. Fontecilla-Camps, J. Biol. Inorg. Chem., 2005, 10, 239–249. 4 D. L. DuBois, Inorg. Chem., 2014, 53, 3935–3960. 5 S. B. Dementin, F. Leroux, L. Cournac, A. L. de Lacey, ´ger, B. Burlat, N. Martinez, S. Champ, A. Volbeda, C. Le L. Martin, O. Sanganas, M. Haumann, V. M. Fernandez, B. Guigliarelli, J. C. Fontecilla-Camps and M. Rousset, J. Am. Chem. Soc., 2009, 131, 10156–10164.
24032 | Phys. Chem. Chem. Phys., 2014, 16, 24026--24033
Paper
6 E. Garcin, X. Vernede, E. C. Hatchikian, A. Volbeda, M. Frey and J. C. Fontecilla-Camps, Structure, 1999, 7, 557–566. 7 Y. Higuchi, H. Ogata, K. Miki, N. Yasuoka and T. Yagi, Structure, 1999, 7, 549–556. 8 Y. Higuchi, T. Yagi and N. Yasuoka, Structure, 1997, 5, 1671–1680. 9 F. Leroux, S. Dementin, B. Burlatt, L. Cournac, A. Volbeda, S. Champ, L. Martin, B. Guigliarelli, P. Bertrand, J. FontecillaCamps, M. Rousset and C. Leger, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 11188–11193. 10 M. C. Marques, R. Coelho, A. L. De Lacey, I. A. C. Pereira and P. M. Matias, J. Mol. Biol., 2010, 396, 893–907. 11 P. M. Matias, C. M. Soares, L. M. Saraiva, R. Coelho, J. Morais, J. Le Gall and M. A. Carrondo, J. Biol. Inorg. Chem., 2001, 6, 63–81. 12 H. Ogata, S. Hirota, A. Nakahara, H. Komori, N. Shibata, T. Kato, K. Kano and Y. Higuchi, Structure, 2005, 13, 1635–1642. 13 H. Ogata, P. Kellers and W. Lubitz, J. Mol. Biol., 2010, 402, 428–444. 14 A. Volbeda, M.-H. Charon, C. Piras, E. C. Hatchikian, M. Frey and J. C. Fontecilla-Camps, Nature, 1995, 373, 580–587. 15 M. Stein and W. Lubitz, Phys. Chem. Chem. Phys., 2001, 3, 2668–2675. 16 R. Cammack, Nature, 1999, 397, 214–215. 17 D. M. A. Smith, Y. Xiong, T. P. Straatsma, K. M. Rosso and T. C. Squier, J. Chem. Theory Comput., 2012, 8, 2103–2114. 18 M. Kampa, W. Lubitz, M. van Gastel and F. Neese, J. Biol. Inorg. Chem., 2012, 17, 1269–1281. 19 M. P. Shaver, L. E. N. Allan, H. S. Rzepa and V. C. Gibson, Angew. Chem., Int. Ed., 2006, 45, 1241–1244. 20 R. L. Yson, J. L. Gilgor, B. A. Guberman and S. A. Varganov, Chem. Phys. Lett., 2013, 577, 138–141. ¨mer, M. Kampa, W. Lubitz, M. van Gastel and 21 T. Kra F. Neese, ChemBioChem, 2013, 14, 1898–1905. 22 P. E. M. Siegbahn, J. W. Tye and M. B. Hall, Chem. Rev., 2007, 107, 4414–4435. 23 M. Bruschi, G. Zampella, P. Fantucci and L. De Gioia, Coord. Chem. Rev., 2005, 249, 1620–1640. 24 H. Li, A. D. Robertson and J. H. Jensen, Proteins: Struct., Funct., Bioinf., 2005, 61, 704–721. 25 Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. M. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. M. Wang and P. Kollman, J. Comput. Chem., 2003, 24, 1999–2012. 26 M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys., 2000, 112, 8910–8922. 27 A. L. Cheng and K. M. Merz, J. Phys. Chem., 1996, 100, 1927–1937. 28 G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101. 29 E. Lindahl, B. Hess and D. van der Spoel, J. Mol. Model., 2001, 7, 306–317. 30 A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
This journal is © the Owner Societies 2014
View Article Online
Published on 30 September 2014. Downloaded by The University of Manchester Library on 31/10/2014 14:26:07.
Paper
31 C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789. 32 M. J. Frisch, J. A. Pople and J. S. Binkley, J. Chem. Phys., 1984, 80, 3265–3269. 33 A. Schafer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577. 34 C. H. Chang and K. Kim, J. Chem. Theory Comput., 2009, 5, 1137–1145. 35 M. Stein and W. Lubitz, Curr. Opin. Chem. Biol., 2002, 6, 243–249. 36 M. G. Delcey, K. Pierloot, Q. M. Phung, S. Vancoillie, R. Lindh and U. Ryde, Phys. Chem. Chem. Phys., 2014, 16, 7927–7938. 37 L. M. Lawson Daku, F. Aquilante, T. W. Robinson and A. Hauser, J. Chem. Theory Comput., 2012, 8, 4216–4231.
This journal is © the Owner Societies 2014
PCCP
38 C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280. 39 M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. Van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus and W. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489. 40 K. C. Gross, P. G. Seybold and C. M. Hadad, Int. J. Quantum Chem., 2002, 90, 445–458. 41 J. S. Hub and B. L. de Groot, PLoS Comput. Biol., 2009, 5, e1000480, DOI: 10.1371/journal.pcbi.1000480. 42 D. W. Scott, Biometrika, 1971, 66, 605–610. 43 Y. Montet, P. Amara, A. Volbeda, X. Vernede, E. C. Hatchikian, M. J. Field, M. Frey and J. C. FontecillaCamps, Nat. Struct. Biol., 1997, 4, 523–526. 44 V. H. Teixeira, A. M. Baptista and C. M. Soares, Biophys. J., 2006, 91, 2035–2045.
Phys. Chem. Chem. Phys., 2014, 16, 24026--24033 | 24033
![Photochemical dihydrogen production using an analogue of the active site of [NiFe] hydrogenase.](/assets/img/hold.png)
