Aerobic Damage to [FeFe]-Hydrogenases: Activation Barriers for the Chemical Attachment of O2

[FeFe]-hydrogenases are the best natural hydrogen-producing enzymes but their biotechnological exploitation is hampered by their extreme oxygen sensitivity. The free energy profile for the chemical attachment of O2 to the enzyme active site was investigated by using a range-separated density functional re-parametrized to reproduce high-level ab initio data. An activation free-energy barrier of 13 kcal mol−1 was obtained for chemical bond formation between the di-iron active site and O2, a value in good agreement with experimental inactivation rates. The oxygen binding can be viewed as an inner-sphere electron-transfer process that is strongly influenced by Coulombic interactions with the proximal cubane cluster and the protein environment. The implications of these results for future mutation studies with the aim of increasing the oxygen tolerance of this enzyme are discussed.

1 Details of the QM calculations: geometry optimisation All geometry optimisations were carried out with BP86 functional [1] and def2-TZVP basis [2] as implemented in Turbomole [3] program. Additionally, the empirical dispersion corrections in form proposed by Grimme were included [4] with the B-J dumpig scheme. [5] Resulting functional is denoted as BP86+D3 throughout the manuscript. Geometry optimisations were performed with the default settings (SCF convergence 10 -6 , grid size m3). All single point energies were calculated with tightened SCF convergence (10 -7 ) and enlarged grid size ('m4' in the Turbomole nomenclature and 'fine' in NWChem). We took the advantage of the resolution of identity approximation [6] in all calculations. The activation was followed with constrained optimisations where Fe d …O 1 distance was varied between 3.5 Å to 1.8 Å and all C α and nitrogen atoms in NH 3 + groups of LYS 322 and 359 (Cp numbering, 201 and 238 in case of Dd) remained frozen.

Small model and CA1-B3LYP functional
The small model structure consist of an isolated [2Fe] H cluster where the [Fe 4 S 4 ] cubane was replaced with H + ion and bridging cysteine was substituted with a smaller CH 3 -Sfragment (see Figure S1). We followed the energy change with different methods in a series of singlepoint calculations on the top of optimised structures where the Fe…O distance was varied between 3.5 Å and 1.8 Å. In Figure S1 we show most the most important geometry changes that accompanies the oxygen binding. The oxygen-bound state was found at the Fe…O distance of 1.95 Å with O 1 -O 2 bond length of 1.29 Å. Such lengthening of a bond between two oxygen atoms together with calculated vibration frequency of 1186 cm -1 indicates that the oxygen moiety is best described as a superoxo entity O 2 -. [7] Figure S1. Geometry changes upon oxygen activation on the small cluster model.

2
The reference data were obtained with NEVPT2 method which base on the CASSCF reference wave function. [8] Because of higher basis-set demands of the wave function-based methods the def2-TZVPP basis set was employed which features additional polarization functions on Fe and S atoms. Although a full-valence CAS was not possible, ligands like CO/CN-create so strong ligand field that most of the doubly occupied and empty 3d orbitals of both iron atoms could remain inactive. We demonstrate this effect in Figure S2 where a qualitative scheme of orbital energies is presented. In fact, d z2 and d xz orbitals of Fe d , two π* orbitals of O 2 molecule and one π* orbital of the bridging CO along with five electrons [CAS(5,5)] were sufficient to describe the binding process (MOs enclosed in the teal rectangle in Figure S2; in Figure S3 the natural orbital isosurfaces are presented). In order to stabilise the active space over whole potential energy surface scan we employed stateaveraging approach in which the CASSCF wave function was optimised for three doublet and two quartet states. Because selected active space can be considered as minimal, we performed additional MR-DDCI2 with small def2-SVP basis and we found that at each Fe…O distance the weight of the CASSCF reference was not smaller than 0.85 with other weights lower than 0.01. Figure S2. Qualitative molecular orbitals scheme of the small model system with active space orbitals marked in teal rectangle. Figure S3. Isosurfaces of the natural orbitals obtained in the state-averaged CASSCF (5,5) calculations for Fe…O distance of 2.4 Å.
The reference calculations were compared with various functionals and in Figure S4a the one dimensional potential energy surface cut calculated with various DFT approaches is presented. Functionals tested include semi-local BP86, global hybrids B3LYP, [1a,9] BHandLYP [1a,b,d,8] with 20% and 50% of the Hartree-Fock exchange, respectively, and two range-separated (CAM-B3LYP, [10] LC-ωPBE [11] ). The latter functionals base on the Ewald split of the Coulomb operator: [12] 1 r ଵଶ = 1 − [α + β • er fሺμr ଵଶ ሻ] r ଵଶ + α + β • erf (μr ଵଶ ) r ଵଶ where r ଵଶ is the interelectronic distance, erf is an error function and μ is range-separation parameter that is typically between 0.1 and 1. The parameter α defines a fixed amount of Hartree-Fock exchange (HFX) at all values of r ଵଶ while β regulates the variable amount. The sum of both these parameters give the maximum portion of HFX at long distance. The main differences between CAM-B3LYP and LC-ωPBE, apart from the form of the correlation part, are in the values of α, β and μ parameters which are α = 0.19, β = 0.46, μ = 0.33 for the former and α = 0.0, β = 1.0, μ = 0.4 for the latter.
In case of global hybrid functionals the barrier increases with the amount of HFX but the binding becomes less favourable. BHandLYP functional shows that the process is endothermic by 0.2 kcal/mol. Some improvement over B3LYP is its range-separated counterpart for which the barrier is 2.8 kcal/mol. Unfortunately, the error binding energy (3.9 kcal/mol) is still large. LC-ωPBE performes similar to BHandLYP. Thus we decided to tune the CAM-B3LYP functional to reproduce reference NEVPT2 activation barrier and binding energy. To decrease the number of parameters we set α + β = 1 so in the limit of infinite interelectronic separation the functional features 100% of HFX. We calculated the energies of bound state, transition state at Fe…O distance of 2.4 Å and complex at long Fe…O distance of 3.5 Å with α values between 0.0 and 0.4 and μ between 0.1 and 0.9. For each combination we calculated the relative error in ∆E and ∆E ‡ with respect to reference NEVPT2 data (-7.6 kcal/mol and 5.6 kcal/mol, respectively). In the next step, we took the sum of absolute values of these errors as a measure of total error. For each value of μ we computed the average of such total errors for all five α values between tested. Obtained data can be found in Figure S5.
We see that for μ = 0.5 a minimum has been found. The lowest error for this value of rangeseparation parameter was calculated for α = 0.1. Thus the final optimised set of parameters are α = 0.1, β = 0.9, μ = 0.5. We also note that our relatively large range-separation parameter was reported to provide better activation barriers [13] what is now confirmed in our calculations. The new functional is denoted as CA1-B3LYP. We also carried out calculations on non-covalent complexation energies test set (NCCE31/04 test set of Zhao and Truhlar [14] ) in order to evaluate the performance of CA1-B3LYP functional for interaction energies. The statistical evaluation can be found in Table S1. In most cases the new functional outperforms BP86 and B3LYP, especially for charge-transfer-dominated complexes. It performs surprisingly well for weak interactions as well what will be important for calculations on the large clusters. We note however that we set α + β to 1 and possible improvement is offered by an independent optimisation of all three parameters.
5 Figure S4. Relative energy change upon oxygen binding to the distal iron of the small cluster calculated with various methods. For each method energy at R = 3.5 Å is taken as a reference. Figure S5. The dependency of the sum of unsigned relative errors with respect to rangeseparation parameter µ. The reference curve (black solid line) in Figure 2b was obtained with a doublet ground-state state-specific CAS configuration interaction (CASCI) wave function with the state-averaged CASSCF orbitals presented graphically on the right (orbital labels taken from Turbomole program).

Large cluster model and ONIOM calculations
The cluster models were constructed from X-Ray structures of Cp and Dd proteins (PDB codes: 3C8Y and 1HFE, respectively). Consistently with recent experimental and theoretical studies, [15,16] the di-µ-dithiolato bridge was modelled as SCH 2 (NH)CH 2 S. Residues included in the model are listed in Table S2. In both models all residues that can form a N-H…S(Cys) hydrogen bonds were included (cut-off 4 Å). Cuts through covalent bonds were terminated with hydrogen atoms. The H ox state has a total spin of 1/2 originating from strong antiferromagnetic coupling of 18 electrons within the cubane (two high-spin Fe III and two high spin Fe II ) that in turn couple weakly with one unpaired electron at the [2Fe] H site ( Figure   1b). The geometry optimisations were performed with usual set-up. In the ONIOM calculations described in the main body of the manuscript the hydrogen atoms were used as a link atoms between high-and low-layers. The distance scaling factors for C-C and C-N(amide) bonds were 0.709 and 0.729, respectively. The starting orbitals for calculations were converted with an in-home written Perl script from converged Turbomole alpha and beta ASCII files to input files that can be used with asc2mov program (part of official NWChem 6.3 distribution) to 8 obtain binary orbital files of NWChem. In this way we kept the same broken-symmetry coupling in all calculations.

Intermediate model set-up for Gibbs' free energy corrections calculations and explanation of the role of the [Fe 4 S 4 ] cubane
Because we kept some of the atoms fixed in our model and approximate nature of the transition states the direct usage of the harmonic approximation to account for Gibbs' free energy corrections was not possible. We decided to construct a model (blue and red atoms in Figure 1a) for which the Gibbs free energy corrections can be evaluated easily. The system comprised from entire H-cluster along with anchoring CYS residues that were replaced with CH 3 -Sfragments as usual. The oxygen molecule was optimised separately. We noted that the oxygen molecule already lost some translation entropy due to penetration of the protein and according to MD simulations [17a] this effect is about 4 kcal/mol. In our calculations this can be estimated by noting that O 2 molecule lost some translation entropy due to entering into the binding pocket (average radius of ~3 Å 17 ). In standard condition one oxygen molecule occupies ܸ = 3.73 • 10 ିଶ m ଷ while the protein pocket volume is ܸ = 1.13 • 10 ିଶ଼ m ଷ . We thus scaled the translation partition function to account for this reduced volume what resulted in the -TS contribution of 10.4 kcal/mol at 298.15 K. Further corrections for zero-point energy change and enthalpy change were small: 1.5 kcal/mol and -1.0 kcal/mol, respectively. Thus, we found that 9.9 kcal/mol need to be added to the electronic binding energy to account for the free energy change. The same shift was applied to the electronic activation energy what constitute an upper limit for ∆G ‡ .
We observed dramatic change (~10 kcal/mol) between binding energies of the small ( Figure   2) and large models ( Figure 3). Thus we looked at the binding energies in the intermediate model (one used for ∆G corrections) and found that ∆E goes down to -20.6 kcal/mol in comparison to the small model (to -5.2 kcal/mol). Moreover, the binding energy becomes less favourable by 3 -5 kcal/mol for the large models in comparison to the medium system. This is mainly due to inclusion of the counter-charges and partially due to some less pronounced screening effects.

Calculations of oxygen inactivation rates
Within the steady-state approximation, the rate for ligand binding is given by where k +1 is the bi-molecular rate constant for ligand diffusion from the solvent to the active site cavity, k -1 is the rate constant for the reverse process and k 2 is the rate constant for chemical attachment of the ligand initially located in the active site cavity.
Converting the computed free energies into k 2 using transition state theory and adopting values for k +1 and k -1 from previous MD simulations for [NiFe]-hydrogenase, [22] we obtain values for k in of 3.6 s -1 mM -1 and 1.2 s -1 mM -1 for Cp and Dd enzymes, respectively. We believe that this should provide good first approximation to the corresponding diffusion rates in [FeFe]-hydrogenases because for small k +2 k in is proportional to k +1 /k -1 which in turn only depends on the size of the active site cavity (for the same concentration).