Hydrogen-Activation Mechanism of [Fe] Hydrogenase Revealed by Multi-Scale Modeling

When investigating the mode of hydrogen activation by [Fe] hydrogenases, not only the chemical reactivity at the active site is of importance but also the large-scale conformational change between the so-called open and closed conformations, which leads to a special spatial arrangement of substrate and iron cofactor. To study H2 activation, a complete model of the solvated and cofactor-bound enzyme in complex with the substrate methenyl-H4MPT+ was constructed. Both the closed and open conformations were simulated with classical molecular dynamics on the 100 ns time scale. Quantum-mechanics/molecular-mechanics calculations on snapshots then revealed the features of the active site that enable the facile H2 cleavage. The hydroxyl group of the pyridinol ligand can easily be deprotonated. With the deprotonated hydroxyl group and the structural arrangement in the closed conformation, H2 coordinated to the Fe center is subject to an ionic and orbital push-pull effect and can be rapidly cleaved with a concerted hydride transfer to methenyl-H4MPT+. An intermediary hydride species is not formed.

[Fe] hydrogenase requires the substrate methenyl-H 4 MPT + (see Fig. 1), which acts as the hydride acceptor. The question is therefore why the iron cofactor iron-guanylylpyridinol (FeGP) (Fig. 1) is required for catalysis [11,15] even though it does not seem to be redox-active.  16. Both molecules are depicted as parametrized for molecular-dynamics simulations (the cysteinate ligand is modeled by a methylthiolate).
An accurate mechanistic description of the H 2 activation process must thus be able to account for the intriguing role of the metal cofactor. Yang and Hall were the first to investigate the mechanism computationally, using a truncated active-site model in an electrostatic continuum [17]. The first main step of the catalytic cycle is the heterolytic H 2 cleavage, with the proton transferred to either the oxypyridine ligand (deprotonated pyridinol) or the cysteinate ligand, which need to be present in the deprotonated form. The second main step is hydride transfer to methenyl-H 4 MPT + [17], which is the rate-limiting step [17]. However, the theoretical description with density-functional-theory (DFT) methods is sensitive to the incorporation of empirical dispersion corrections and the energetics of all elementary reaction steps can be manipulated by firstshell ligand modifications [18]. Yang and Hall formulated the product of H 2 cleavage as a bound dihydrogen species with an elongated, polarized H-H bond, Fe(II)· · · H δ− -H δ+ · · · − O. This species is the resting state in their catalytic cycle [17]. This intermediate could, however, also be described as a hydride complex [18,19].
The difficulties fully to reconcile the experimental observations with that mechanism, which was derived based on a small model, point to the need for an extended theoretical treatment of the reaction. Any mechanistic proposals should be compatible with the following experimental data: (i) If the key intermediate was a stable hydride species, the electronic structure of the Fe atom would change. However, experimental and theoretical Mössbauer spectra indicate that no stable hydride or H 2 -bound species exists under turnover conditions [12,19,20]. (ii) Model compounds that accurately mimic the first coordination sphere of the iron do not even bind H 2 [21][22][23], hence, the protein environment is likely to be crucial. (iii) In model compounds, protonating the thiolate ligand leads to dissociation of the ligand [24]. This indicates that also in the enzyme, the thiolate might not be a viable proton acceptor. (iv) Mutating histidine 14 to alanine reduces the catalytic activity to 1 % of the wild-type level [16,25]. The mechanism should thus explain why His14 is important for the catalytic activity.
Already in 2009, Hiromoto et al. suggested that a large-scale protein motion might play a role in catalysis [16]. The dimeric protein has three subunits (see Fig. 2): The central subunit, which is formed from the intertwined Cterminal domains of both monomers, and two identical peripheral subunits [26]. The peripheral subunits each harbor an FeGP cofactor [25]; methenyl-H 4 MPT + binds to the central subunit [16]. The protein can adopt two conformational states, referred to as open and closed, respectively. In the open conformation, there is a cleft between the central and the peripheral subunits, as shown in Fig. 2 [25,27]. Hiromoto  To model the reaction accurately, a method is required that incorporates the geometrical constraints imposed by the protein and, if possible, also the electronic polarization exerted by the environment. Combined quantum mechanics/molecular mechanics (QM/MM) fulfils both these requirements [28,29]. A further complication is that a crystal structure of the closed conformation is only available for the apoenzyme [26], while the holoenzyme could only be crystallized in the open conformation [25,27]. A crystal structure of the enzyme in complex with the substrate is available only for the C176A mutant in the open conformation. Thus, suitable starting structures for the wild-type holoenzyme-substrate complex in the open and closed conformations must be  To investigate the crucial H 2 cleavage step, we generated such starting structures of the enzyme-substrate complex in the open and closed conformations. In the subsequent molecular-dynamics (MD) simulations, the FeGP cofactor and methylene-H 4 MPT (compare Fig. 1) were parametrized with the General Amber Force Field (GAFF) [30]. We chose to simulate the enzyme-product complex because methylene-H 4 MPT is more straightforward to parameterize with GAFF than methenyl-H 4 MPT + . According to the principle of microscopic reversibility [31], this corresponds to the product structure directly after hydride transfer, and the sampled configurations are relevant for both reaction directions. The Fe atom, both CO ligands, the cysteinate S atom, the acyl CO of the pyridinol ligand, and the oxygen atom of the bound water were positionally restrained to avoid the need for Fe-L bonded parameters. As FeGP is strongly bound to the peripheral subunit, these restraints effectively lock the hinge motion that would interconvert open and closed conformations. However, this conformational change is likely to take place on a time scale much longer than the sampling times used here. The protein was described with the Amber ff03 force field [32,33]. MD simulations were run with Gromacs 4.5.5 [34][35][36][37] (open conformation: 100 ns, closed conformation: 95 ns). Starting from snapshots of the MD trajectory, the H 2 splitting reaction was investigated by QM/MM calculations. These were performed with ChemShell [38][39][40] interfaced to Turbomole [41,42] as QM back-end. QM calculations used the TPSS-D3 [43,44] DFT method with the def2-TZVP [45] basis set for iron and the def2-SVP [46] basis set on all other atoms. The effect of a larger basis set was assessed for one reaction step and the differences in energies and structures were found to be negligible. Details on model construction and computational methods can be found in the electronic supporting information (ESI).
Herein, we use a full model of the dimeric enzyme in molecular-dynamics simulations and QM/MM calculations to address two central questions relating to the H 2 -activation mechanism in [Fe] hydrogenase. These are the protonation state of the FeGP cofactor and the possible H 2 -activation pathways in the closed conformation.

Molecular dynamics simulations 2.1 Open conformation
The MD simulations of the dimer with both monomers in the open conformation yield insights into the dynamics around the cofactor in this non-reactive conformation. The methylene-H 4 MPT molecules (one bound to each monomer) stay attached to the central subunit of the protein throughout the simulation, mainly due to hydrogen-bonding interactions. Most of these interactions were already identified in the crystal structure [16] and remained largely stable throughout the MD trajectory. Important hydrogen bonds are formed between the 2a-amino group of the pterin unit (see Fig. 3) to the backbone carbonyls of Thr317 and Cys250. The carbonyl group of Cys250 also forms a hydrogen bond with the pterin N3-H. The hydroxyl of Ser320 occasionally forms a hydrogen bond to the pterin 2a-amino group and the carbonyl of Ser320 with the pterin N8-H. The hydroxyl group of Ser254 occasionally also engages in a hydrogen bond to the pterin 2a-amino group. The tail of methylene-H 4 MPT is highly flexible and mainly involved in hydrogen bonds to surrounding water molecules. One relatively stable hydrogen bond is formed between Lys151 and either of the glutarate carboxylates (but also to the phosphate). The tail can adopt an extended conformation, and occasionally the glutarate carboxylates form hydrogen bonds with the distant residues Asn153, Lys154, Lys182 or Lys131. However, the predominant conformation of the tail is U-shaped, with the bend at the ribose. Snapshots from the MD simulations with methylene-H 4 MPT in either conformation are shown in Fig. 3.
The different conformational behavior of the head and tail parts of bound methylene-H 4 MPT are mirrored in the RMSD (root-mean-square deviation) with respect to the starting structure. The evolution of the RMSDs of the head and tail parts of both independent methylene-H 4 MPT molecules over the whole trajectory are plotted in Fig. 4. The RMSD of the head fluctuates between 1 to 2Å, so this part of the molecule remains essentially fixed. In contrast, the RMSD of the conformationally flexible, very mobile tail is around 5Å (U-shaped tail), with values up to 10Å corresponding to the extended conformation.
In the open conformation, the Fe center is exposed to the solvent. The hydroxyl group of the pyridinol ligand mainly hydrogen-bonds to the water molecule coordinated to Fe (whose oxygen atom was positionally restrained). It  can also form hydrogen bonds to bulk water molecules. Interestingly, there is a relatively abundant conformation where the hydroxyl group forms a hydrogen bond to His14. His14 is known to be crucial for high catalytic rates, since a H14A mutation reduces the turnover rate to 1 % of the wild-type level [25]. The presence of this hydrogen-bonded conformation suggests that His14 may act as a base to deprotonate the pyridinol ligand, as previously suggested [16,25]. The distance between the hydroxyl proton and N of His14 for both monomers is plotted in Fig. 5. In monomer 1, this hydrogen bond is formed frequently at the beginning of the trajectory but is no longer present beyond 36 ns, whereas in monomer 2, it was mainly observed later during the simulation (see Fig. 5).

Closed conformation
The simulation of the closed conformation samples the conformations of the protein in the state that is believed to be reactive [16]. The reduced substrate methylene-H 4 MPT is bound in the active-site cleft with the hydride-accepting C14a atom in spatial proximity to the Fe atom of FeGP. The geometrical arrangement of methylene-H 4 MPT and the iron center is stable throughout the simulation. The mean distance between Fe and C14a increased from approximately 3.8Å (0 to 30 ns) to around 4.3Å (30 to 90ns) in monomer 1, while it remained constant at approximately 4.3Å throughout the simulation of monomer 2. Notably, methylene-H 4 MPT blocks a water channel identified in the crystal structure [16]. However, a few water molecules are still able to enter the active site by passing along the cofactor already during the first few nanoseconds of the MD simulation. This fast access of a few water molecules indicates that water is able to enter the active-site region through the cavity in the closed conformation.
The hydrogen bond between the pyridinol hydroxyl and His14, already observed in the open conformation, is formed in the closed conformation as well. The distance between the hydroxyl proton and N of His14 for both monomers is plotted in Fig. 6. In monomer 2, this hydrogen bond (OH-N distance of about 2Å) is frequently formed and broken over the course of the simulation. In monomer 1, the OH-N distance plot shows three stages (see Fig. 6). In the first phase, up to 30 ns, a direct hydrogen bond is often formed. From 30 to 51 ns, the OH-N distance remains at around 4.5Å. In this phase, a water molecule that entered the active site bridges the hydroxyl group and N (OH-HOH-N ). Thus, water-mediated proton transfer should still be possible. In the third phase, from 51 to 88 ns, the hydrogen bond is lost and re-forms again only after 88 ns with one bridging water molecule. Hence, in both monomers, deprotonation of the pyridinol with His14 as the base should be viable. The proton transfer may be mediated by a water molecule bridging between the hydroxyl and the proton-accepting N . To investigate possible H 2 activation mechanisms, we first need to clarify the protonation state of the active site. The experimentally verified importance of His14 for a high turnover rate [16,25] and the hydrogen bond between His14 and the pyridinol OH observed in the MD simulations point to a crucial role of His14 as the base in the proton transfer pathway. Deprotonation of the hydroxyl group results in a potent proton acceptor (oxypyridine) for heterolytic H 2 cleavage. To investigate the energetics of pyridinol deprotonation, we chose two representative MD snapshots that feature the OH-His14 hydrogen bond: One snapshot from the open conformation (at 10.78 ns) and one from the closed conformation (at 13.2 ns). Both snapshots were prepared for QM/MM optimization, i.e., the full protein plus a water shell around one of the active sites was extracted (see ESI for details). The QM region contained the FeGP cofactor up to the phosphate linker (with an Fe-bound water in the open conformation) and the His14 side chain; see  In the closed conformation, the proton transfer is endothermic by +2.3 kcal/mol. From a potential-energy surface (PES) scan along the proton-transfer coordinate (defined as the difference between O-H and H-N bond lengths), we estimate an upper bound for the proton-transfer barrier of +4.6 kcal/mol (+2.3 kcal/mol for the back reaction). Thus, proton transfer between the pyridinol OH and His14 is facile, with the OH/His form being favoured. Considering that His14 is connected to the bulk solvent through a proton-transfer chain, we conclude that the less favoured oxypyridine (O − /HisH + ) form is still present in significant amounts under equilibrium conditions.
In the open conformation, the proton transfer is thermoneutral (∆E = 0.0 kcal/mol). The oxypyridine form is thus equally likely. Although we did not calculate the reaction barrier for this case, it is reasonable to assume that it will be similar to the barrier in the closed conformation as the proton transfer reactions are the same in both cases, except for a slight change in the environment. The stabilization of the oxypyridine form in the open conformation compared to the closed conformation arises because the water molecule coordinated to iron can form a hydrogen bond to the oxypyridine oxygen, stabilizing the anion. Note that the active site in the open conformation is exposed to the bulk solvent and thus filled with water (see Fig. 8). For the H 2 activation to proceed, the bound water must be displaced by H 2 . Based on all available data, we cannot assess with any certainty if this happens while the enzyme is in the open or the closed conformation. As we have found that the active site is still accessible to water in the closed conformation (see Sect. 2.2), it is certainly possible for H 2 to enter the active site only after the closed conformation has formed. What is clear, however, is that the prevailing protonation state of the pyridinol/His14 pair will critically depend on the external pH.

H 2 activation
To investigate hydrogen cleavage and hydride transfer to methenyl-H 4 MPT + , we chose two representative snapshots from the closed conformation: one with a short Fe-C14a distance of 3.7Å (at 11 ns) and one with a longer distance of 4.3Å (at 56.5 ns). Because the hydride is transferred to C14a of methenyl-H 4 MPT + , one might expect the reaction to be facilitated by a short Fe-C14a distance, and our discussion thus focuses first on the former snapshot. The QM region included again FeGP up to the phosphate linker, together with the chemically relevant part of the substrate and H 2 (see Fig. 7). In the selected snapshot, the pyridinol-OH-His14 hydrogen bond is not present. Note that His14, in the neutral form, is in the MM region and does not directly participate in the reaction. Considering that the pyridinol-His14 hydrogen bond is frequently formed and broken (see Fig. 6) and that the proton transfer is kinetically facile (see Sect. 3.1), this choice of setup sustains two scenarios: (i) The pyridinol ligand has been deprotonated via His14, and the proton removed from the active site through the proton-transfer chain, leaving behind oxypyridine and neutral His14. (ii) The pyridinol ligand remains neutral, without hydrogen-bonding to (also neutral) His14.

H 2 activation via oxypyridine.
For the scenario with oxypyridine, we studied several possible hydrogen coordination modes to the open coordination site at the Fe center: end-on and two rotamers of side-on coordination. All initial structures converged to a side-oncoordinated hydrogen molecule, which is thus the reactant for the hydrogen cleavage reaction. The structure is shown in Fig. 9. The coordinated H 2 is activated, its bond being elongated to 0.83Å (from 0.74Å in free H 2 ). There is only one reasonable pathway to cleave H 2 in this configuration: In a concerted heterolytic cleavage step, the proton is transferred to the oxypyridine oxygen and the hydride to C14a of methenyl-H 4 MPT + . This reaction is exothermic by −18.7 kcal/mol. A PES scan (along the difference of O-H and H-H bond lengths) provided an upper bound for the barrier of about +1.0 kcal/mol. Despite various attempts, we were unable to locate a stable minimum on the PES that would correspond to an iron hydride species. Hence, we find that the iron is involved in H 2 binding and activation, but does not bind a hydride species. The H 2 cleavage mechanism we have identified here thus complies with the first of the requirements formulated in Sect. 1.
In the reactant complex, the coordinated H 2 is subjected to an electronic push-pull effect from the negatively charged oxypyridine oxygen and the positively charged carbocation of methenyl-H 4 MPT + . This is reflected in the relevant frontier orbitals (Fig. 10): The LUMO has a strong contribution from the p z orbital on C14a, which is oriented perpendicular to the ring plane. The HOMO−2 is delocalized over the oxypyridine ring and the thiolate S atom, with a strong contribution from the oxypyridine oxygen. (Note that HOMO and HOMO−1 are strongly localized on the phosphate linker and thus do not contribute to the reactivity at the iron center).
In the optimized product structure resulting directly from the PES scan, the O-H bond of the pyridinol hydroxyl points towards the empty coordination site of the iron center (Fig. 9). The newly formed C14a-H bond is pointing towards the Fe atom. It is slightly elongated (1.19Å compared to 1.10Å for the other C14a-H bond), indicating a weak interaction with the iron center. Another conformer, where the O-H bond has turned away from the iron, was found to be 2.3 kcal/mol more stable.
Similar results were obtained for the second snapshot, where the substrate was positioned slightly further away from the iron center (4.3 vs. 3.7Å). The cleavage reaction in that case is even more exothermic (by −26.3 kcal/mol, compared to −18.7 kcal/mol for the first snapshot), and again no iron hydride species could be optimized. We thus find that fluctuations of the substrate position of the order as they were observed in the MD simulations have only a minor effect on the reactivity of [Fe] hydrogenase in the H 2 cleavage step.

H 2 activation via thiolate.
In the second scenario, we consider the cleavage reaction with a neutral pyridinol ligand. There are two relevant reactant conformers, which differ in the orientation of the pyridinol O-H bond (Fig. 11). In the following, we quote energies relative to the favoured conformer with the O-H pointing away from the iron. The second conformer, where the O-H is pointing towards the co-ordinated H 2 , is 5.7 kcal/mol less stable. Direct H 2 splitting in the favoured conformer, with the pyridinol OH acting as the proton acceptor, is not possible: Product-like starting structures, with one hydrogen atom already transferred to C14a of methenyl-H 4 MPT + , are not stable minima but converged back to reactant structures during optimization. For the second conformer, this pathway is precluded in the first place by the orientation of the pyridinol OH bond.
However, we find for both reactant conformers that hydrogen cleavage can occur with the thiolate ligand, rather pyridinol OH, as the proton acceptor. When the OH group is oriented away from the Fe center, the resulting iron hydride structure is not stable but the hydride is directly transferred to methenyl-H 4 MPT + to form the product. This reaction is exothermic by −4.4 kcal/mol, significantly less so than hydride cleavage to oxypyridine-O − . For the reactant conformer with the OH bond oriented towards the Fe center, a stable iron hydride intermediate could indeed be located (Fig. 11). It is only +0.3 kcal/mol higher in energy than the favoured reactant conformer. The hydride Fe-H bond length is 1.61Å, in excellent agreement with Fe-H bonds in comparable hydride complexes optimized in vacuo [19] (1.60Å). OH rotation, which is likely to have a low activation barrier, triggers the transfer of the hydride from iron to methenyl-H 4 MPT + , yielding the same product as direct H 2 splitting from the preferred conformer (see Fig. 11). The thiolate is thus able to act as the base, which may provide an explanation for the 1 % remaining activity of the H14A mutation [16,25].
Remarkably, the Fe-SH bond in the product (2.34Å) is even slightly shorter than the Fe-S − bond in the reactant (2.36Å). In the product, the thiol proton forms a short hydrogen bond (1.37Å) to the pterin carbonyl group of methylene-H 4 MPT (see Fig. 11), which is a very good hydrogen-bond acceptor because of its conjugation with the guanidine moiety in the pterin ring. This in turn weakens the S-H bond (elongated to 1.51Å, compared to 1.39Å in the hydride intermediate), which may be described as "partial deprotonation" of the thiol. The formal thiol ligand in the product is similar in character to a thiolate in terms of its interaction with the metal center. The thiol-pterin hydrogen bond thus makes the thiolate a better proton acceptor in the H 2 splitting step, stabilizes the thiol product, and also makes the thiol a better ligand, preventing it from dissociating like in model complexes. The hydrogen bond is enabled by the exact positioning of the FeGP cofactor and the methenyl-H 4 MPT + substrate in the active site. As water molecules are still able to access the active site in the closed conformation (see Sect. 2.2), we can envisage that the excess proton on the thiol ligand is removed from the active site via water.

Discussion and Conclusions
[FeFe] and [NiFe] hydrogenases cleave or form H 2 by redox chemistry [7][8][9]; a basic group close to the active iron atom in [FeFe] hydrogenases is important to donate or accept protons. The mechanism of hydrogen activation in [Fe] hydrogenase is different. The enzyme has two large-scale conformations, which  [47]. However, the oxypyridine plays an additional, crucial role in activating H 2 : It is close to the iron atom and represents an ideal Lewis base. On the other side of the iron is the carbocationic C14a of the substrate methenyl-H 4 MPT + , which is an ideal Lewis acid. Furthermore, both groups are ionic. When a hydrogen molecule coordinates to the iron, it is polarized by these charges and subjected to an electronic push-pull effect exerted by the Lewis pair. The spatial arrangement in the closed conformation is exactly such that the coordinated H 2 lies in-between C14a + and O − . This leads to facile, exothermic heterolytic H 2 cleavage, without involving electron transfers to/from the metal center. The Lewis acid C14a + is only present in proximity to FeGP when methenyl-H 4 MPT + is bound in the closed conformation. H 2 cleavage in the open conformation is thus unlikely.
The activation mechanism we have described is reminiscent of hydrogen activation by frustrated Lewis pairs [48]. The hydrogen-bound adduct does not need to be very stable since the H 2 cleavage barrier is extremely low (about 1 kcal/mol). Hence, any H 2 binding event can directly lead to H 2 cleavage, without requiring a long-lived H 2 -bound intermediate.
When the pyridinol ligand is not deprotonated, it is still possible to split H 2 via proton transfer to the thiolate ligand. However, we have found this pathway to be much less favorable. This is consistent with observations in biomimetic model complexes that thiol is a poor ligand [24]. The pyridinol/oxypyridine equilibrium must be strongly affected by the pH, so we would expect the reactivity to depend critically on pH as well, which is indeed the case [1,49].
The atomistic mechanism of H 2 activation in [Fe] hydrogenase we have proposed herein satisfies the criteria set out in Sect. 1: No stable hydride intermediate; no occurrence of, or requirement for, a long-lived H 2 adduct; no involvement of the thiolate ligand as a proton acceptor; a crucial role for His14. In our preferred mechanism, the pyridinol hydroxyl group and His14, together with the stable placement of the substrate carbocation in the active site, are the essential players, which is in accord with the observation that the enzyme looses 99 % of its activity upon H14A mutation [25]. The residual activity of the H14A mutant can be explained by the alternative, less favorable activation pathway via the thiolate.
In the open conformation, which might be prevailing in solution, the waterbound FeGP cofactor is the most probable form, which agrees with the results of a theoretical Mössbauer study [19]. To be able to study [Fe] hydrogenase by means of molecular-dynamics (MD) simulations and quantum-mechanics/molecular-mechanics (QM/MM) calculations a force-field description of the whole system is necessary. We chose the Amber force field (version ff03 1,2 ) for the molecular-mechanics (MM) part. This force field can be combined with the general Amber force field (GAFF), 3 which allows for an easy parametrization of organic molecules. We followed the GAFF parametrization procedure to derive parameters for the iron-guanylylpyridinol (FeGP) (see Fig. 1 in the main paper) cofactor and the substrate methylene-tetrahydromethanopterin (methylene-H 4 MPT, see Fig. 1 in the main paper). The parametrization can be divided into two steps. As first step, partial charges derived from the electrostatic potential (ESP) for all atoms were calculated. Then, GAFF topology files for the organic parts were created. The parameters used are given in Figs. 1, 2 and Tables 1, 2.

Calculation of partial charges
For the calculation of partial charges, we chose a procedure similar to that applied in Ref.  15 and an aug-cc-pVTZ basis set 16-18 on all atoms was employed. Unconstrained structure optimization of methenyl-H 4 MPT (starting structure from PDB file 3H65 19 ) leads to the formation of internal hydrogen bonds of the tailing carboxylate groups to hydroxo groups of either the furanose ring or the glycerin-derived part. These internal hydrogen bonds are not found in the protein structure because the carboxylates can form intermolecular hydrogen bonds. Hence, intramolecular hydrogen bonds do not resemble the bonding situation in the crystal structure and need to be avoided. This was achieved by constraining several dihedral angles during optimization.

Generation of GAFF topologies
The generation of GAFF topologies for methylene-H 4 MPT was straightforward. Topology files were created with the program ACPYPE, 20 which interfaces Antechamber 21 of AmberTools 13 22 to create Amber and Gromacs topology files. The whole methylene-H 4 MPT molecule was treated as one MM-residue.
Generation of suitable parameters for the FeGP cofactor was more involved. 23 Since metal complexes are not parametrizable with GAFF the cofactor was split up into the guanylylpyridinol ligand, 2 CO molecules, the iron ion (each treated as one MM-residue) and the cysteine rest. For the iron ion a new ion type was introduced. The two CO molecules were parametrized with GAFF. The guanylylpyridinol ligand could be parametrized with GAFF, however, the parameters of the acyl part coordinating the iron atom had to be adjusted (equilibrium value of the acyl O-C-C angle) (see Fig. 1 in the main paper for the structure). In the crystal structure, the open coordination site is coordinated by water. For this water molecule, the standard parameters of the TIP3P water model 24 were used as the ESP-calculated partial charges differed insignificantly from the standard charges. For the cysteinate rest, a modified cysteinate residue was defined, which had the same bonded and van der Waals parameters as a standard cysteinate, but ESP-derived charges. With a view on terminating the QM region at the cysteinate C β -C α bond in the subsequent QM/MM calculations, while maintaining integer charges for the QM and MM regions, a charge of +0.069588e was distributed equally among the C β H 2 atoms of the new cysteinate residue to obtain a total integer charge of −1e for FeGP. Since the Fe atom, both CO molecules, the cysteine sulfur atom, and the CO atoms of the acyl-ligand were positionally restrained during the classical MD simulations, no metal-ligand bonded parameters had to be derived.     could not be refined. They form an unordered tail, far away from the reactive centers, and thus were simply omitted. To prepare the structure for the simulations, several modifications had to be made. FeGP is ligated by dithiothreitol in the crystal structure, which leads to a slight displacement of the whole cofactor compared to its position in the wild-type structure (without methylene-H 4 MPT, PDB code: 3F47 25 ). To obtain a wild-type-like structure Ala176 was mutated  back to cysteinate and the dithiothreitol ligand was removed from the structure. To place the cofactor in the correct (wild-type) position, the structure of the wild type (3F47) was aligned with the structure of the C176A mutant (3H65) using PyMOL. 26 The coordinates of the FeGP moiety (including the cysteine S atom) in the aligned wild-type structure were then inserted into the modified mutant structure, replacing those of the misplaced FeGP. With this procedure, the position of FeGP in the modified mutant structure resembled the position in the wild-type structure. Only the C β -S bond of Cys176 (formerly Ala176) was stretched from 1.79Å in the wild-type structure to 2.21Å in the modified mutant structure. The correct bond length is restored in the energyminimization step at the start of the MD simulations. The protonation states of titratable residues were determined with the program reduce, 27 which can also correct flipped Asn/Gln/His residues (none were in this case), and verified with propKa. 28 All titratable residues were in their standard protonation states; His protonation states are summarized in Table 4. Finally, the whole protein dimer was created from the monomer chain with PyMOL. 26

Closed conformation
The closed conformation was generated from the modified crystal structure of the open conformation (see previous section), which we call starting structure here. The only crystal structure available in the closed conformation is for the wild-type apoenzyme (PDB code: 2B0J 29 ). To build a model of the complete enzyme in the closed conformation, the central subunit (residues 253-345) of 2B0J was first aligned to the central subunit of the starting structure. In the second step, the peripheral subunit (residues 1-241) of the starting structure was aligned to the peripheral subunit of the now aligned 2B0J. The FeGP cofactor was aligned together with the peripheral subunit. The structure thus obtained for the closed conformation of the holoenzyme-substrate complex had no significant atom overlaps, except for the tail part of methylene-H 4 MPT (after the ribitol part, see Fig. 1 of the main paper). This tail was rotated with the help of the UCSF Chimera program 30 to remove atom overlaps. Given the high flexibility and mobility of the tail, any memory of the initial conformation will be lost during the MD sampling. In the structure thus generated, there is a gap in the backbone between residues 241 and 242, between the hinge region (residues 242-252) and peripheral subunit (see Fig. 3), which will be closed during energy minimization. Finally, the water molecule coordinated to Fe was removed because it would prevent the hydride transfer reaction that we intend to study. With the water molecule absent, the structure is exactly the product of the hydride transfer. The full dimer was again created with PyMOL. 26 A superposition of the starting structure, the final structure, and the apoenzyme in the open conformation is presented in Fig. 3. In the generated structure of the closed conformation, the distance between Fe and the hydride accepting carbon atom (C14a) is 3.23Å, which compares well to the distance of 3Å found by Hiromoto et al., who modelled the structure of the closed conformation in a similar fashion. 19

MD simulation protocol
All MD preparation steps and simulations were performed with the Gromacs molecular dynamics package version 4.5.5. [31][32][33][34] The protein was centered in a triclinic box with a minimal distance of 1 nm between solute and box border. The box was solvated and ions (Na + and Cl − ) were added to neutralize the system and to obtain an ion concentration of 0.15 mmol/L. The Fe atom, both CO ligands, the cysteinate S atom, and the CO group of the Fe-coordinating acyl ligand in the FeGP cofactor were kept frozen or positionally restrained at their positions in the prepared crystal structure. Positional restraints, rather than constraints were necessary for the pressure equilibration because pressure scaling with constrained (frozen) atoms leads to technical difficulties in the Gromacs implementation. The integration time step was 2 fs. The linear constraint solver (LINCS) algorithm to 4th order with 1 iteration was invoked to enforce constraints (all bonds after energy minimization). Water molecules were kept rigid with the SETTLE algorithm in all steps after energy minimization. For neighbor searching, a grid-based group cut-off scheme was used with a cut-off distance of 1 nm for short-range interactions and the neighbor list was updated every 10 steps. Coulomb interactions were calculated with a smooth particle-mesh Ewald algorithm with interpolation order of 4, a Coulomb cut-off of 1 nm, Fourier spacing of 0.12 nm, tolerance of 10 −5 and optimized Fourier transforms. Van der Waals interactions were calculated with a cut-off scheme (radius = 1 nm). Energy and pressure were corrected for long-range dispersion effects. Initial velocities were generated according to a Maxwell-Boltzmann distribution at 293 K. For temperature scaling the system was coupled to a v-rescale thermostat. 35 During equilibration, several subsystems were coupled to their own thermostats. For production, the entire system was coupled to one thermostat with a relaxation time of 2 ps and reference temperature of 293 K. Box equilibration was achieved with isotropic box rescaling by coupling the system to a Berendsen barostat with 1 ps relaxation time, compressibility of 4.5 · 10 −5 bar −1 , center-of-mass scaling of reference coordinates, and a target pressure of 1 bar.
The simulation stages are summarized in Table 5. The system in the open conformation was energy minimized in vacuum and in solvent. During 200 ps heating in an N V T ensemble, the protein, the cofactor, the substrate and the solvent were coupled to a temperature bath while all heavy atoms of the protein were positionally restrained. Thereafter, the box was equilibrated in an N P T simulation where the restraints on protein atoms were reduced in three steps (total 900 ps). In the first step, protein, cofactor, substrate and solvent were coupled to three thermostats. In the second and third steps, only two thermostats were utilized (solvent, rest of the system). After adjustment of the box vectors (see Table 6), the system was equilibrated in an N V T ensemble for 400 ps with two thermostats (solvent, rest of the system) and finally 2 ns with one thermostat. The production trajectory was 100 ns. Coordinates were saved to disk every 20 ps.
The simulation procedure for the closed conformation was similar. Differences are a reduced force threshold for the vacuum minimization, position restraints for the heavy protein atoms during solvent minimization and only two thermostats (solvent, rest of the system) in all box-equilibration simulations. The equilibrated box parameters are collected in Table 6. The final equilibration simulation (N V T ) had to be elongated for an additional 2 ns with a reduced temperature coupling constant and reduced restraints for Fe and its first-shell ligands (see Table 6) to avoid instabilities in the production run. Still, the production run terminated at 94.6 ns. We trace this back to the initial strain put into the system during crystal structure manipulation, which could not fully relax because the Fe atoms of the two FeGP cofactors were positionally restrained. Furthermore, the parameters for FeGP and methylene-H 4 MPT might be not optimal for long simulations. However, the trajectory of 94.6 ns is already longer than trajectories normally produced to prepare QM/MM calculations 36 and is sufficiently long to allow us to analyse visited conformations, structural behavior, and protein dynamics around the active site. Trajectories were analysed with the VMD program. 37

QM/MM setup
For the QM/MM calculations, several snasphots representing important conformations were selected, as described in the main paper. As QM/MM calculations under periodic boundary conditions are not supported by ChemShell, water molecules and ions outside a shell of 18Å around the quantum-mechanics (QM) region were therefore discarded to create a finite system. The entire protein dimer was retained. QM/MM optimizations were carried out with the ChemShell program 38-40 (version 3.5.0) with Turbomole (version 6.5) 12,41 as the QM back-end. The optimizations were performed in hybrid delocalized internal coordinates (HDLCs) using the HDLCOpt module. 42 The scaling factor for Cartesian coordinates when constructing HDLCs was set to 0.8; the interval to update the pair list and regenerate HDLCs was set to 100 steps; the convergence threshold was set to 0.001 a.u.
Several regions were defined for the QM/MM optimizations. The QM region contained all quantum-mechanically described atoms (52 to 84 atoms). The MM region contained all other atoms of the system (approximately 12350 to 12900 atoms). Only the atoms in the active region (around 670 to 1700 atoms) were allowed to move in structure optimizations. We used a QM/MM microiterative optimization scheme 43, in which the inner region (around 61 to 93 atoms) contained the QM atoms, the MM boundary atoms, and the MM atoms bonded to the them.
The QM part of the calculations was treated with the TPSS exchangecorrelation functional 44-46 plus Grimme's DFT-D3 dispersion correction. 47 Structures were optimized with the def2-TZVP basis set 48 on iron and the def2-SVP  Run was necessary to avoid system collapse at the production stage basis set 49 on all other atoms. The resolution of the identity approximation with corresponding auxiliary basis sets 50 was invoked to speed up the calculations. 51 The same force-field setup as in the MD simulations was used for the MM region. QM-MM electrostatic interactions were calculated fully (no cut-off), and the QM and MM systems were coupled with the charge-shift scheme. 39,[52][53][54] Where the QM-MM boundary cut through a covalent bond, the QM region was saturated with a H link atom. MM partial charges of residues cut by a QM-MM boundary were corrected such that the QM and MM parts had integer charges. The correction charge (which is usually of the order of 0.01e) was distributed equally over the entire residue involved, leading to very minor differences between the MM and QM/MM partial charges for these residues (see Tables 1, 2,  3). Two different QM regions were defined (see Fig. 7 of the main paper). The first QM region contains the Fe center, both CO ligands, the side chain of the iron-coordinating cysteine (Cys176), the gunaylylpyridinol ligand up to the phosphate linker and the hydride-acceptor up to the phenyl part. This first region was used to study H 2 activation reactions at the Fe center and hydride transfer to methenyl-H 4 MPT + . The second region did not contain the substrate, but instead the side chain of His14. This region was used for the investigation of proton transfer from the pyridinol hydroxyl group to His14.

Effect of the basis set
The H 2 cleavage reaction in the closed conformation with the snapshot at 11 ns was recalculated with a def2-TZVP basis set 48 on all atoms. Also with the larger basis, no stable hydride species could be located, and structure optimizations converged to the reduced substrate. The reaction energy was −16.0 kcal/mol compared to −18.7 kcal/mol with the smaller basis set. Hence, a larger basis set has no significant quantitative or qualitative effect on structures; the small effect on energies leaves conclusions and interpretation unaffected.

Potential energy surface scans
To estimate the reaction barriers for the H 2 cleavage reaction with proton transfer to oxypyridine (snapshot at 11 ns) and of the proton transfer from thy pyridi-nol OH group to His14 in the open conformation (snapshot at 10.78 ns), we calculated energy profiles along scan coordinates that closely resemble the reaction coordinates; the profiles are plotted in Fig. 4. The scans provide an upper bound for the reaction barrier. It was not possible to locate transition states for the two reactions. It is likely that the flat profile of the two reactions caused difficulties in numerical hessian calculations.