## 2 Methods  
Three conformers of [NiFeSe] hydrogenase from _D. vulgaris_ Hildenborough were observed in the crystallographic data from PDB ID: 5]SH based on an oxidized species [1]. A model was therefore built for each of the three conformers by selecting the coordinates of the active site of each of the three conformers, which will henceforth be referred to as conformer I, II, and III (Figures 2 and 3). These conformers differ not only in the number of sulfur atoms and their bonding patterns at the Ni site, but also in the geometry of the selenocysteine. Hydrogen atoms were modeled using CHARMM with the charmm36 force field [19; 20]. His82 was protonated at position N\(\varepsilon\) due to possible interactions with the Fe[CN]\({}_{2}\)CO unit. Unless otherwise specified, protonation of amino acid side chains was according to the standard assignment for pH 7. Missing atoms were modeled with CHARMM using the charmm36 force field [19; 20]. The entire protein, including the catalytic center and three [Fe\({}_{4}\)S\({}_{4}\)] clusters, and crystal water molecules were then solvated in a water box with sodium chloride ions (0.1 M). Partial charges for the three [Fe\({}_{4}\)S\({}_{4}\)] clusters were adapted from previous work [13; 14; 15]. Energy minimizations were performed to optimize the hydrogen bond networks at the protein-solvent interface. In a final step, all protein atoms, crystal water, solvent water, and salt ions within a radius of 40 A centered on the Ni atom were selected to comprise the complete model (Figure 4A). Figure 3: The positions of atoms optimized with quantum mechanical (QM) methodology include the Fe-Ni active site and the surrounding amino acids shown in color (carbon atoms in cyan, nitrogen atoms in blue, oxygen atoms in red, and sulfur and selenium atoms in yellow). Hydrogen atoms of protein are omitted for clarity. The QM/MM geometry optimization of each of the six models was carried out using a model that divides the enzyme into three regions: QM, MM active, and MM inactive (Figure 4B). Only the positions of atoms within the QM and MM active regions were optimized, whereas atoms in the MM inactive region were held fixed. Atoms included in the QM region are all atoms involved in the Fe-Ni catalytic center, as well as atoms belonging to amino acids Ala81, Ala420, Ala444, Arg422, His82, and Leu425. Pro421, Ser443, and Ser445, two water molecules (w365 and w589) from the crystal structure, were also included in the QM region (Figure 3). Atoms belonging to protein and crystal water (oxygen) within a 15 A radius of the Ni atom were considered as part of the MM active region. Atoms in the QM region were treated with the Density Functional Theory (DFT) and the BP86 functional with the resolution of identity (RI) technique, as implemented in Turbomole6.2 [21, 22, 23]. A mixed basis set of def2-TZVP (Ni, Fe, S, and Se atoms) and 6-31G\({}^{*}\) (C, N, O, and H atoms) was used; this combination of basis sets has been tested and shown to be appropriate for geometry optimizations and vibrational frequency analyses on such enzymes [14]. Based on EPR studies indicating that [NiFeSe] hydrogenase is nickel-EPR silent in the oxidized state [1], the computations presented here treat Ni in the 2+ oxidation state. Furthermore, a common feature of hydrogenases is the presence in the active site of a low-spin Fe coordinated by CO- and CN-ligands [24, 25]. Therefore, the electronic states considered in the models are based on [Ni\({}^{2+}\)Fe\({}^{2+}\)] and [Ni\({}^{2+}\)Fe\({}^{3+}\)] configurations. For Fe\({}^{2+}\), the total charge of the QM system is \(-1\) for all three conformers and the total spin is S = 0 (multiplicity 1). For Fe\({}^{3+}\), the total charge of the system is 0 for all three conformers, giving a total spin of S = 1/2 and a spin multiplicity of 2. The energetic coupling between the first and second layers (QM and mobile MM) was treated using electrostatic embedding with a charge-shift scheme as implemented in ChemShell [26, 27, 28].