Understanding 2D-IR Spectra of Hydrogenases: A Descriptive and Predictive Computational Study

: [NiFe] hydrogenases are metalloenzymes that catalyze the reversible cleavage of dihydrogen (H 2 ), a clean future fuel. Understanding the mechanism of these biocatalysts requires spectroscopic techniques that yield insights into the structure and dynamics of the [NiFe] active site. Due to the presence of CO and CN − ligands at this cofactor, infrared (IR) spectroscopy represents an ideal technique for studying these aspects, but molecular information from linear IR absorption experiments is limited. More detailed insights can be obtained from ultrafast nonlinear IR techniques like IR pump -IR probe and two-dimensional (2D-)IR spectroscopy. However, fully exploiting these advanced techniques requires an in-depth understanding of experimental observables and the encoded molecular information. To address this challenge, we present a descriptive and predictive computational approach for the simulation and analysis of static 2D-IR spectra of [NiFe] hydrogenases and similar organometallic systems. Accurate reproduction of experimental spectra from a ﬁrst-coordination-sphere model suggests a decisive role of the [NiFe] core in shaping the enzymatic potential energy surface. We also reveal spectrally encoded molecular information that is not accessible by experiments, thereby helping to understand the catalytic role of the diatomic ligands, structural differences between [NiFe] intermediates, and possible energy transfer mecha-nisms. Our studies demonstrate the feasibility and beneﬁts of computational spectroscopy in the 2D-IR investigation of hydrogenases, thereby further strengthening the potential of this nonlinear IR technique as a powerful research tool for the investigation of complex bioinorganic molecules.


Introduction
Molecular hydrogen (H 2 ) is an ideally clean fuel that releases large amounts of energy but no greenhouse gases upon combustion. Harnessing this potential for future energy conversion approaches requires green and sustainable strategies for the catalytic evolution and activation of H 2 . Hydrogenases are complex metalloenzymes that accomplish these reactions by using protein-embedded catalytic metal sites equipped with biologically uncommon diatomic ligands. In the case of [NiFe] hydrogenases, the active-site core contains two metal ions, nickel and iron, which are coordinated by a total of four cysteines ( Figure 1A) [1]. Two of them serve as terminal ligands to the Ni, while the other two bridge the two metals ions. In addition, the Fe ion is coordinated by one CO and two CN − ligands that tune the catalytic properties of the active site in ways that are not fully understood yet ( Figure 1A,B) [2][3][4][5][6][7]. Hydrogenases can inspire the rational design of compact and power (bio-)synthetic catalysts for the transformation of H 2 , given that their structural, d namical, and functional properties are thoroughly understood. As a consequence, sp troscopic techniques play a key role in hydrogenase research. Due to the presence of C and CN − ligands, infrared (IR) spectroscopy has been extensively utilized to study h drogenases since the first experimental observation of these unexpected diatomic ligan with this very technique [2,3,5,6,8,9]. In this respect, IR spectroscopy represents a part ularly powerful method for studying hydrogenases for two major reasons. (1) The d tomic ligands at the catalytic site give rise to bond-stretching vibrations-νCO, νCN and νCN s -that yield sharp and well detectable IR absorption bands in a spectral regim that is free from other vibrational transitions of the protein (ca. 1850 − 2150 cm −1 [NiFe] hydrogenases; Figure 2). The exact frequencies of these stretching vibrations d pend on various properties of the catalytic metal site. Thus, they represent site-speci and highly sensitive IR-spectroscopic probes for nuclear geometry and electronic stru ture [10]. (2) IR spectroscopy is a highly flexible technique that can be used to study h drogenases in various sample forms (e.g., in liquid and frozen solution, crystallized, ophilized, within living cells, or bound to electrodes) , under tight experimen control (temperature, illumination, electrochemical potential) [2,11,13,14,16,18,21-23,2 47], and with high time resolution [48][49][50][51][52][53][54][55]. Hydrogenases can inspire the rational design of compact and powerful (bio-)synthetic catalysts for the transformation of H 2 , given that their structural, dynamical, and functional properties are thoroughly understood. As a consequence, spectroscopic techniques play a key role in hydrogenase research. Due to the presence of CO and CN − ligands, infrared (IR) spectroscopy has been extensively utilized to study hydrogenases since the first experimental observation of these unexpected diatomic ligands with this very technique [2,3,5,6,8,9]. In this respect, IR spectroscopy represents a particularly powerful method for studying hydrogenases for two major reasons. (1) The diatomic ligands at the catalytic site give rise to bond-stretching vibrations-νCO, νCN as , and νCN s -that yield sharp and well detectable IR absorption bands in a spectral regime that is free from other vibrational transitions of the protein (ca. 1850 − 2150 cm −1 for [NiFe] hydrogenases; Figure 2). The exact frequencies of these stretching vibrations depend on various properties of the catalytic metal site. Thus, they represent site-specific and highly sensitive IR-spectroscopic probes for nuclear geometry and electronic structure [10]. (2) IR spectroscopy is a highly flexible technique that can be used to study hydrogenases in various sample forms (e.g., in liquid and frozen solution, crystallized, lyophilized, within living cells, or bound to electrodes) , under tight experimental control (temperature, illumination, electrochemical potential) [2,11,13,14,16,18,[21][22][23], and with high time resolution [48][49][50][51][52][53][54][55].
Despite these advantages, conventional IR absorption spectroscopy suffers from a lack of structural and dynamical insight, inter alia, because variations in CO and CN stretch frequencies can reflect different types of structural and electronic changes. Thus, they are typically difficult to interpret. Simulating IR spectra of different structural scenarios by quantum chemical methods can facilitate the understanding of the experimental data [24,33,[56][57][58], but it cannot solve the fundamental problem, i.e., the inherent limitation of linear IR spectroscopy to few experimental observables (here: the stretching frequencies of the CO and CN − ligands). To overcome this limitation, we recently introduced nonlinear IR techniques to hydrogenase research [55]. In a nutshell, these nonlinear techniques differ from a linear absorption experiment in the type and number of light-matter interactions [59]. A conventional IR absorption experiment involves a single interaction of the sample with the incoming continuous IR beam from a conventional light source. The amount of absorbed probe light depends linearly on the incident intensity, and the result is an absorbance spectrum. Nonlinear IR techniques, in contrast, involve multiple interactions of the sample with ultrashort laser pulses that lead to a sequence of vibrational transitions. The molecular response does not scale linearly with the total incident intensity, and the resultant spectrum depends on the pulse sequence. Two IR pulses yield an IR pump -IR probe (PP) spectrum that reflects the change in probe light absorbance after interaction of the sample with the pump pulse. Adding a third pulse represents one way to record a 2D-IR spectrum, which resolves the PP spectrum in terms of the pump frequency. Both experiments yield detailed and otherwise inaccessible information on structure and dynamics without sacrificing the inherent advantages of IR spectroscopy. For instance, by probing multiple transitions of a single vibrational mode (e.g., the CO stretch vibration) solution-phase insights into bond lengths, strengths, and energies can be obtained from PP and 2D-IR spectra. In addition, transitions involving multiple vibrational modes (e.g., CO and CN stretch vibrations) lead to cross peaks in 2D-IR spectra that inform on the coupling between the involved degrees of freedom and the underlying molecular groups. Besides a fundamental understanding of the shape of the potential energy surface (PES), which ultimately governs the reactivity of the active site, signals in PP and 2D-IR spectra also provide information on the rate and direction of intramolecular energy transfer. This type of information may be relevant for understanding regulatory aspects and possible non-classical effects in enzymatic catalysis. Despite these advantages, conventional IR absorption spectroscopy suffers from a lack of structural and dynamical insight, inter alia, because variations in CO and CN stretch frequencies can reflect different types of structural and electronic changes. Thus, they are typically difficult to interpret. Simulating IR spectra of different structural scenarios by quantum chemical methods can facilitate the understanding of the experimental data [24,33,[56][57][58], but it cannot solve the fundamental problem, i.e., the inherent limitation of linear IR spectroscopy to few experimental observables (here: the stretching frequencies of the CO and CN − ligands). To overcome this limitation, we recently introduced nonlinear IR techniques to hydrogenase research [55]. In a nutshell, these nonlinear techniques differ from a linear absorption experiment in the type and number of light-matter interactions [59]. A conventional IR absorption experiment involves a single interaction of the sample with the incoming continuous IR beam from a conventional light source. The amount of absorbed probe light depends linearly on the incident intensity, and the result is an absorbance spectrum. Nonlinear IR techniques, in contrast, involve multiple interactions of the sample with ultrashort laser pulses that lead to a sequence of vibrational transitions. The molecular response does not scale linearly with the total incident intensity, and the resultant spectrum depends on the pulse sequence. Two IR pulses yield an IRpump-IRprobe (PP) spectrum that reflects the change in probe light ab- hydrogenase and associated fundamental transitions in the IR absorption spectrum, exemplarily shown for the Ni a −S state (spectrum simulated from experimental transition frequencies of ReRH). Partial localization of symmetric and antisymmetric νCN vibrations to cyanide ligands I and II is indicated by orange and green shading, respectively.
While nonlinear IR techniques have proven to monitor the CO/CN stretch reporter vibrations of [NiFe] hydrogenase in detail, information encoded in the complex and unique patterns of signals is not fully understood yet. In the current account, we will demonstrate how computational strategies can be utilized for accurately simulating the vibrational transitions reflected by PP and 2D-IR spectra, thereby allowing to exploit these powerful techniques to their full potential. Of note, our approach differs from previous vibrational analyses performed on [NiFe] hydrogenases in two respects. (1) While linear IR absorption spectra of [NiFe] hydrogenases have been calculated within the harmonic limit by us and others [24,33,[56][57][58]60], this approach is not applicable to nonlinear IR spectra since the signals in these data are a direct consequence of the anharmonicity of the PES. Thus, we will use a perturbative treatment instead, which utilizes a quartic expansion in the normal mode basis [61][62][63][64][65]. (2) Previous vibrational analyses on hydrogenases have focused on a 'fingerprinting' approach that aims to extract structural information by matching experi-mental spectra with those calculated for certain structural scenarios. Here, we rather aim for a 'computational spectroscopy' approach that provides a fundamental understanding of the studied spectroscopic observables, transferrable concepts for interpreting future data, and information on molecular aspects that are experimentally inaccessible. As in a previous experimental study on a regulatory model hydrogenase from Ralstonia eutropha (ReRH) [55], we will focus on the H 2 -accepting intermediate of [NiFe] hydrogenases ( Figure 1A), called Ni a −S, which reflects the core features found in all redox-structural states of the [NiFe] site. To discriminate structurally sensitive observables and highlight spectroscopic trends, the hydride-containing Ni a −C state (the second, yet to be studied intermediate of ReRH; Figure 1B) and molecular complexes mimicking the enzymatic Fe 2+ (CO) CN − 2 moiety ( Figure 1C,D) [10,66,67] will also be analyzed in a predictive manner.

Computational Models
To calculate vibrational transition energies beyond the harmonic limit, we have utilized a small first-coordinate-sphere model of the [NiFe] active site ( Figure 1A,B; bottom; Materials and Methods). Compared to a previously employed model [24,56,60,68], cysteine residues where not modelled as ethanethiolate but reduced in size to methanethiolate. Limiting the size of the model has two important advantages for the current study. (1) A smaller number of atoms N in the computational model reduces the number of 3N − 6 internal degrees of freedom and, in particular, the number of 6N − 11 Hessian matrices required for describing the quartic potential. This way, the high cost for computing anharmonic frequencies can be reduced to an acceptable level. In addition, limiting the size of the model reduces the number of highly anharmonic large-amplitude low-frequency modes (often associated with multi-well potentials) that are not well treated by a second-order perturbative approach employing a quartic expansion of the PES [64]. (2) Focusing on the first coordination sphere of the active site, spectroscopic properties that are determined by the inherent architecture of the [NiFe] core can be distinguished from those that may be dictated by steric constraints or electrostatic effects imposed by the protein.

Harmonic and Fundamental Anharmonic Frequencies
We start our analysis with the Ni a −S state, which is characterized by the consensus architecture of the [NiFe] site, as described in the Introduction section, and a free bridging site for substrate binding between nickel and iron ( Figure 1A). With the chosen level of theory, we obtain harmonic frequencies, ω H , of 1911 and 2068/2078 cm −1 for the CO and CN stretching vibrations, respectively ( Figure 3B: grey line; Table S1). As observed previously [24,56,68], CN stretching frequencies agree well with the experimental values, e.g., of our model enzyme ReRH (2071/2080 cm −1 ), but the CO stretch frequency is clearly underestimated (1942 cm −1 ; Figure 3A). While this observation may be partly related to the lack of the protein environment in the computational model used, we like to stress that harmonic frequencies are, in general, highly sensitive to the level of theory. Thus, calculated absolute frequencies-whether fitting to experimental values or not-are generally not representative, and analyses of calculated vibrational spectra usually focus on empirically scaled values or frequency differences [21,24,33,[56][57][58].
values, e.g., of our model enzyme ReRH (2071/2080 cm ), but the CO stretch frequency is clearly underestimated (1942 cm −1 ; Figure 3A). While this observation may be partly related to the lack of the protein environment in the computational model used, we like to stress that harmonic frequencies are, in general, highly sensitive to the level of theory. Thus, calculated absolute frequencies-whether fitting to experimental values or not-are generally not representative, and analyses of calculated vibrational spectra usually focus on empirically scaled values or frequency differences [21,24,33,[56][57][58]. Signals in harmonic (grey) and anharmonic (black) absorption spectra are represented as Lorentzian lines. Two-dimensional-IR spectra are shown as contour plots for the entire νCX spectral region (center) and the region characteristic for νCN vibrations (bottom). These spectra were obtained by two-dimensional Fourier transform of the third-order response for 〈 〉 polarisation, as detailed in the Materials and Methods section. Since the simulation of the Signals in harmonic (grey) and anharmonic (black) absorption spectra are represented as Lorentzian lines. Two-dimensional-IR spectra are shown as contour plots for the entire νCX spectral region (center) and the region characteristic for νCN vibrations (bottom). These spectra were obtained by two-dimensional Fourier transform of the third-order response for ZZZZ polarisation, as detailed in the Materials and Methods section. Since the simulation of the 2D-IR spectra employs the harmonic approximation for scaling intensities of excited-state transitions, spectra throughout the entire figure are based on harmonic transition dipole strengths. In case of the 2D-IR contour plots, transition dipole strengths of νCN fundamental transitions were multiplied by a factor of 10 for better cognition. Relevant intra-and intermode anharmonicities, ∆ω a 1−2 and ∆ω a/b , are indicated for νCO and νCN vibrations.
As expected, calculated anharmonic frequencies reflecting the fundamental transitions between ground and first vibrationally excited states, ω 0→1 , are lower than their harmonic counterparts, yielding values of 1891 and 2041/2050 cm −1 for the CO and CN stretching vibrations, respectively ( Figure 3B: black line; Table S1). Accordingly, the difference between harmonic and anharmonic frequencies, ∆ω H−1 = ω H − ω 0→1 , which reflects the anharmonicity close to the bottom of the vibrational potential, assumes values of 20 and 27/28 cm −1 , respectively ( Figure 3B; Table S1). On the one hand, analyzing these frequency differences is advantageous since they are not affected by inaccuracies in the harmonic frequencies. On the other hand, these quantities cannot be directly determined experimentally (in a sense, the harmonic frequency is a hypothetical quantity rather than a real observable), and they are not specifically related to nonlinear IR spectroscopy. Moreover, both ω 0→1 and ∆ω H−1 depend on the accuracy of several cubic and quartic force constants, which may be limited if the molecular model features low-frequency large-amplitude vibrations (vide supra). Thus, we will not dwell on these quantities in the following.

Reproduction of Experimental Transitions: The Role of the [NiFe] Core
Signals in PP and 2D-IR spectra reflect vibrational transitions beyond single-mode first excited states, and two general scenarios can be distinguished in third-order nonlinear spectroscopy ( Figure S1) [59]. (1) After being excited from the ground to the first excited vibrational state by the pump beam, further excitation of the pumped vibrational mode a (e.g., νCO) by the probe beam will populate its second excited state. This process will yield signals at transition frequencies ω 0a→1a = ω a 0→1 and ω 1a→2a = ω a 1→2 , where indices indicate the number of vibrational quanta in mode a before and after the probed vibrational transition. (2) Alternatively, the probe beam may excite another vibrational mode b (e.g., νCN s )-not only from the shared vibrational ground state but also from the first excited state of the initially pumped mode a. This type of process will yield signals at transition frequencies ω 0a0b→0a1b = ω 00→0b = ω b 0→1 and ω 1a0b→1a1b = ω a0→ab . Both types of signals can be extracted from 2D-IR spectra, and the mixed-mode transitions are typically only accessible via this technique ( Figure S1). In the following, we will analyze all these transition frequencies, and spectra simulations will focus on the 2D-IR spectroscopy. To this end, we will evaluate frequency differences ∆ω a 1−2 = ω a 0→1 − ω a 1→2 (e.g., ∆ω νCO 1−2 ) and ∆ω a/b = ω 00→0b − ω a0→ab (e.g., ∆ω νCO/νCN s ), summarized in Table S1, because both ∆ω a 1−2 and ∆ω a/b can be directly read from 2D-IR (and PP) spectra ( Figure 3), and they have clear physical meanings. ∆ω a 1−2 is the intramode anharmonicity of mode a, which reflects the shape of the associated vibrational potential, and ∆ω a/b is the intermode anharmonicity, which reflects the coupling between modes a and b. Furthermore, like ∆ω H−1 , both ∆ω a 1−2 and ∆ω a/b are insensitive towards inaccuracies of the computed harmonic frequencies. In addition, and in contrast to ∆ω H−1 , calculation of ∆ω a 1−2 cancels out many contributions from potentially spurious cubic and quartic force constants that could affect the accuracy of absolute frequencies ω a 0→1 and ω a 1→2 . For the Ni a −S state of the [NiFe] active site, we obtain a value of 23 cm −1 for the intramode anharmonicity ∆ω νCO 1−2 of the CO stretch mode (mode number #64; Figure S2), which is very close to the experimental value of 25 cm −1 ( Figure 3A,B and Figure 4A,B) [55]. Likewise, we calculate ∆ω νCN as 1−2 = 19 cm −1 and ∆ω νCN s 1−2 = 9 cm −1 (modes #65 and #66, Figure S2), in good agreement with the experimental values of 18 cm −1 and 8 cm −1 , respectively ( Figure 3A,B and Figure 4A,B) [55,69]. This indicates that our calculations yield an accurate description of the individual mode potentials, and the real vibrational eigenstates are sufficiently well approximated by the normal modes, justifying the perturbative approach used here. This statement is further supported by an accurate reproduction of the pronounced coupling between the two CN stretch modes, ∆ω νCN as /νCN s = 21 cm −1 , again close to the experimental value of 24 cm −1 ( Figure 3A,B and Figure 4A,B) [55,69]. In total, these findings demonstrate that the computational model used and the employed perturbative approach are able to capture the anharmonicities that are extractable from PP and 2D-IR spectra of [NiFe] hydrogenase, in particular the unusually different intramode anharmonicities of the two CN stretch modes and the strong coupling between them. Thus, a good overall agreement of 2D-IR spectra simulated from experimental and calculated transition energies is obtained ( Figure 3A,B). If a large basis set is used for all atoms (see the Materials and Methods), the agreement is even better, with all mentioned anharmonicities being reproduced with an error ≤ 1 cm −1 . Except for the CO stretch anharmonicity, these quantities were computed prior to their extraction from experimental data, highlighting the predictive powers of this computational spectroscopy approach.
accurate description of the individual mode potentials, and the real vibrational eigenstates are sufficiently well approximated by the normal modes, justifying the perturbative approach used here. This statement is further supported by an accurate reproduction of the pronounced coupling between the two CN stretch modes, ∆ CN as / CN s = 21 cm −1 , again close to the experimental value of 24 cm −1 (Figures 3A,B and 4A,B) [55,69]. In total, these findings demonstrate that the computational model used and the employed perturbative approach are able to capture the anharmonicities that are extractable from PP and 2D-IR spectra of [NiFe] hydrogenase, in particular the unusually different intramode anharmonicities of the two CN stretch modes and the strong coupling between them. Thus, a good overall agreement of 2D-IR spectra simulated from experimental and calculated transition energies is obtained ( Figure 3A,B). If a large basis set is used for all atoms (see the Materials and Methods), the agreement is even better, with all mentioned anharmonicities being reproduced with an error ≤ 1 cm −1 . Except for the CO stretch anharmonicity, these quantities were computed prior to their extraction from experimental data, highlighting the predictive powers of this computational spectroscopy approach.   From the agreement between experimental and gas-phase calculated anharmonicities, we conclude that the general signal pattern in the 2D-IR spectrum of the [NiFe] active site is predetermined by the intrinsic architecture of the [NiFe] core, as is the case with the general shape of the linear IR absorption spectrum. On the one hand, this indicates that the marked difference between the intramode anharmonicities of the two CN stretch modes and the strong coupling between them cannot be explained by the entatic-state nature of Ni a −S or hydrogen-bonding interactions of the active-site CN − ligands (neither is included in the computational model used here) [12,55,60,68,70,71]. On the other hand, the exact values of experimentally observed anharmonicities may well be affected by these aspects as well as structural details of different [NiFe] intermediates (vide infra) and other subtle interactions with the protein environment, as is also the case with the linear IR absorption spectrum [3,10,70]. Given the solvent sensitivity of 2D-IR spectra recorded from other cyanido metal complexes [72], the unexpectedly exact agreement between enzyme-derived experimental anharmonicities and gas-phase calculated values may point towards a protein environment that is tailored to preserve the intrinsic properties of the [NiFe] center.
As noted previously [55,69], direct coupling between the CO and CN stretch modes is very weak and mainly evident for νCN s . Due to partial cancelling of positive and negative signals, the exact coupling strength is hard to extract from experimental data but probably on the order of about 2 cm −1 for the latter mode ( Figures 3A and 4A) [69]. In line with this statement, both ∆ω νCO/νCN as and ∆ω νCO/νCN s are calculated to be negligible within the accuracy of the employed approach ( Figures 3B and 4B). At this point, we cannot exclude, though, that the weak coupling observed in the experiment is affected by subtle interactions between the [NiFe] cofactor and the protein environment. Support for this hypothesis comes from the observation that this particular coupling is not reproduced by the firstcoordination-sphere model even if a large basis set is used for all atoms (vide supra).

Experimentally Inaccessible Insights into CO Bonding
Given the successful reproduction of experimentally observable quantities, we will next focus on aspects that are not directly accessible with experiments. Since their interaction is weak, we will separately analyze CO and CN stretching modes in order to understand in detail what kind of molecular information can be extracted from the associated vibrational transitions.
We will start our analysis by focusing on the CO stretching mode, which-according to harmonic vibrational analyses and (site selective) isotope exchange studies-is strongly localized to the CO bond coordinate and, thus, assumed to be weakly affected by other vibrational degrees of freedom [5,6,24,56,60,73,74]. To test this hypothesis, we performed anharmonic analyses in a reduced normal mode space [62] with all cubic and quartic force constants excluded except for those associated with the CO and CN stretch mode coordinates (∆ω νCO 1−2 ) or the CO stretch mode coordinate only (∆ω νCO 1−2 ). Surprisingly, we found that the resulting partial intramode anharmonicity ∆ω νCO 1−2 = ∆ω νCO 1−2 = 15 cm −1 differs considerably from the experimental value and the value that was obtained by the fully dimensional calculation (cf. Figure 4A-D respectively). This indicates that the experimentally determined value does not just reflect the anharmonicity of the PES in direction of the CO bond coordinate but a more complex quantity related to the real vibrational eigenstate transitions that is not fully captured in the normal mode picture. To understand this phenomenon in more detail, we next analyzed which of the other vibrational degrees of freedom have to be reincluded in the anharmonic analysis to reproduce the experimental intramode anharmonicity. Interestingly, we found that this value can be almost completely reproduced (∆ω νCO 1−2 = 23 cm −1 ; Figure 4E) by including lowerfrequency ( 350 ∼ 650 cm −1 ) metal-ligand modes associated with νFe-CO/CN stretching and δFe-CO/CN bending (#31−#39, Figure S3) in addition to νCO and νCN vibrations. While including all modes with significant contributions from Fe-CO/CN coordinates to the respective vibrational eigenvectors was necessary to fully capture the ∆ω νCO 1−2 value from the fully dimensional analysis, surprisingly good agreement was also found by in-  Figure S3). Similar observations were also made for the Ni a −C state of [NiFe] hydrogenase and, to a lesser extent, the model complex (Table S1). This indicates that the observed phenomenon represents a general feature of heteroleptic monocarbonyl complexes, as long as there is no considerable mixing with other ligand vibrational modes.
Based on the illustrated findings, a couple of important conclusions regarding the interpretation of vibrational transitions associated with the CO stretching mode of [NiFe] hydrogenase and similar molecular systems containing diatomic ligands can be drawn.
(1) Based on experimental data, we previously demonstrated that the sequence of vibrational transitions (up to the fourth exited level) follows the trend expected for a diatomic Morse potential. As a consequence, key bond properties like the equilibrium bond strength (force constant f ), the total bond strength (dissociation energy D 0 ), and the non-linear elasticity of the bond (anharmonicity constant a) can be extracted from these transition energies [55]. In the light of our computational findings, it should be stressed, though, that all these quantities reflect the properties of the bond within its dynamic molecular context, i.e., the ad hoc Morse potential represents an effective potential including some contributions from other molecular degrees of freedom (vide supra) rather than a true single-bond-coordinate potential. (2) While the quantities extracted from this effective Morse potential, as obtained experimentally, are those we are typically interested in when analyzing chemical bonds, complementary insights into the isolated properties of the bond (not convoluted with those of other molecular coordinates) are beneficial in cases where a more detailed understanding of bonding and intramolecular interactions is required. Naturally, these insights are impossible to obtain experimentally since the probed vibrational transitions reflect the total intramode anharmonicity ∆ω νCO 1−2 , which is a property of the real vibrational eigenstates. As demonstrated above, however, the partial intramode anharmonicity ∆ω νCO 1−2 , which is a property of a selected vibrational coordinate (here: νCO), can be calculated, and quantification of individual contributions to the experimentally accessible ∆ω νCO 1−2 is possible. This finding highlights the potential of computational 2D-IR spectroscopy as a genuine research tool that goes beyond a mere reproduction of experimental spectra. (3) Our calculations indicate a tight interplay between intra-ligand νCO and metal-ligand Fe-CO/CN vibrations. This interplay appears to be most pronounced for νCO and νFe-CO, supporting the general idea that CO bonding in metal carbonyl complexes is tightly controlled by Fe-CO bonding. Since intermode coupling between νCO and νFe-CO is negligible (vide infra), the effect of νFe-CO on ∆ω νCO 1−2 must be conveyed by the cubic force constant f aab , where a = νCO and b = νFe-CO. In other words, the contribution of νFe-CO to ∆ω νCO 1−2 is a measure for the dependence of the CO equilibrium bond strength, f aa , on the length of the Fe-CO bond, d b .

Computational Spectroscopy of the CN Stretch Vibrational Manifold
The two CN − ligands of the [NiFe] site give rise to two bond-stretch vibrational modes. Based on pioneering experimental work involving partial isotope labelling [6] and supported by harmonic vibrational analyses [24,56,60], these modes are generally assigned to a higher-frequency symmetric and a lower-frequency antisymmetric CN stretch mode, νCN s and νCN as , respectively. As a consequence (and in contrast to νCO), the associated transition frequencies cannot be ad hoc analyzed in terms of individual CN − bond properties, and their structural interpretation is generally more complicated.
While all studied dicyanido complexes exhibit a similar pattern of intramode and intermode CN stretch anharmonicities, the exact values are not identical. Specifically, we observe that the difference in calculated intramode anharmonicities, ∆∆ω νCN 1−2 = ∆ω νCN as 1−2 − ∆ω νCN s 1−2 , increases from the Ni a −C state (∆∆ω νCN 1−2 = 5 cm −1 ) via Ni a −S (∆∆ω νCN 1−2 = 10 cm −1 ) to model complex 1 (∆∆ω νCN 1−2 = 13 cm −1 ). Likewise, the intermode anharmonicity, i.e., the coupling between the two CN stretch modes, ∆ω νCN as /νCN s , increases in the order Ni a −C < Ni a −S < 1 (vide supra). We therefore conclude that both ∆ω νCN as /νCN s and ∆∆ω νCN 1−2 represent measures for the interaction between the two CN stretch modes. To establish a fundamental understanding of the strength of this interaction and the encoded structural information, we next take a closer look at the normal mode eigenvectors of νCN as and νCN s , which describe the atom movements in these vibrational modes. While the two CN stretch modes of [NiFe] hydrogenase can be well described as symmetric and antisymmetric combinations of the two CN − ligands (vide supra) [6,24,56,60], this representation is an approximation of the real mode eigenvectors (and even more so of the real eigenstates beyond the normal mode picture). In reality, each of the two modes is somewhat localized on each one of the two cyanide ligands, i.e., CN − I in cases of νCN as and CN − II in cases of νCN s (Figures 2 and S2). To quantify this partial localization, we calculated the scalar projection of the normal mode eigenvectors on idealized single-bondstretch unit vectors (see the Computational Details section). Here, a projection value of 1 corresponds to a perfectly localized single-bond-stretch vibration, while a value of ca. 0.7 indicates a vibration that is fully delocalized across the ligands (idealized symmetric or antisymmetric stretching of both ligand bonds). Using this analysis, we obtain a projection of 0.  Figure 4D and Table S1). As expected, partial intramode anharmonicities differ clearly from their experimentally accessible counterparts, further confirming the strong interaction between νCN as and νCN s . In particular, ∆ω νCN as 1−2 is negative because an antisymmetric stretching mode features a symmetric potential with steep walls on both sides, in contrast to the Morse potential applicable to single-bond stretching. By contrast, ∆ω νCN s 1−2 is positive but still smaller than expected for single-bond stretching (cf. Figure 4B-E) since the shape of the potential for symmetric two-bond stretching is similar but steeper than a single-bond-stretch Morse potential. Thus, increasing mode localization (delocalization) should increase (decrease) both values by approaching (diverging from) a Morse-like potential. Indeed, we obtain ∆ω νCN as  (Table S1). In the semi-localized limit, ∆ω νCN as  1−2 , thereby allowing to assess the bond potentials of the individual ligands. Of note, these partial intramode anharmonicities are not accessible experimentally, highlighting the importance of computational 2D-IR spectroscopy, e.g., in hydrogenase research.
In the preceding paragraphs, we have illustrated the structural information that can be extracted from the intermode anharmonicity ∆ω νCN as /νCN s and the difference in intramode anharmonicities ∆∆ω νCN 1−2 . So far, however, it is unclear why the two intramode anharmonicities, ∆ω νCN as 1−2 and ∆ω νCN s 1−2 , are different in the first place. To address this question, we had another look at individual contributions to these quantities, focusing on vibrational resonances, which are treated variationally in the generalized second order vibrational perturbation approach used here (GVPT2; see the Materials and Methods) [61][62][63][64][65]. Interestingly, we found a 2-2 Darling-Dennison (DD) resonance that involves the second excited states of the two νCN modes. If this 2-2 DD resonance is excluded from the vibrational analysis (together with other resonances of this type), the difference between ∆ω νCN as 1−2 and ∆ω νCN s 1−2 vanishes for all systems studied in this work (∆∆ω νCN 1−2 < 1 cm −1 ; Figure 4I-L). This demonstrates that the non-resonantly interacting modes have equal anharmonicities, and it is the resonance that increases (decreases) the energy of the second excited state of νCN s (νCN as ), thereby decreasing (increasing) the difference between ω 0→1 and ω 1→2 . This way, ∆ω νCN s 1−2 becomes smaller and ∆ω νCN as 1−2 becomes larger, leading to the characteristic difference in intramode anharmonicities. While other effects on these anharmonicities cannot be excluded, our results therefore indicate that the unequal values are primarily a consequence of quantum-mechanical mixing rather than a direct reflection of the shape of the PES. Thus, ∆∆ω νCN 1−2 will be large if there is a strong resonant interaction between the second excited states of the two νCN modes, which requires near degeneracy of these states. Consistently, we observe that ∆∆ω νCN 1−2 is a strictly monotonic function of the splitting of the non-resonantly interacting second excited states as well as the splitting of the fundamental frequencies. By the same logic, we can also understand why equal anharmonicities were observed for a previously studied dicarbonyl complex (vide supra) [75][76][77]: For this complex, the splitting of the fundamental νCO as and νCO s frequencies (ca. 70 cm −1 ) is more than seven times larger than observed for the νCN as and νCN s modes of [NiFe] hydrogenase (ca. 10 cm −1 ), which results in a large energy difference between the corresponding second excited states of ca. 140 cm −1 . This large splitting excludes resonant interactions, thereby preserving equal anharmonicities in the two νCO modes in that case. We also note that intermode anharmonicities ∆ω νCN as /νCN s and resonance-free intramode anharmonicities ∆ω νCN 1−2 are strictly anticorrelated for all studied systems ( Figure 4I-L). This observation reflects the transition between weak and strong coupling regimes [59], which is, again, a consequence of changes in molecular symmetry. Of note, the resonance-free anharmonicities cannot be determined experimentally, but they can be reconstructed by taking the mean of the two unequal anharmonicities measured for νCN as and νCN s (compare Figure 4B,F-H and I-L, respectively).

Vibrational Networks and Energy Transfer
So far, we have focused on the experimentally probed CO and CN stretch vibrations, their interactions, and the encoded molecular information. Given that both νCO and (to a lesser extent) νCN transition frequencies are affected by other normal coordinates, we will also explore their interactions with these vibrational modes by analyzing the corresponding intermode anharmonicities.
Focusing on the Ni a −S state, we find that both νCO and νCN are coupled to several lower-frequency modes between ca. 350 and 650 cm −1 (#31-#39; Figure S3) that are dominated by νFe-CO/CN and δFe-CO/CN coordinates (vide supra). Maximum coupling strengths of > 2 cm −1 are predicted ( Figure 5), but since intermode coupling seems to be slightly underestimated in our calculations (vide supra), we assume that the real coupling strengths may reach values close to 5 cm −1 . As expected, νCO exhibits the strongest coupling to modes #35-#38, which are dominated by Fe-CO motions, while the νCN vibrations interact most strongly with modes #31, #32, #35, and #36, which exhibit increased contributions from Fe-CN motions ( Figure S3). Of note, Fe-CO/CN modes also interact with each other, with maximum coupling strengths of > 3 cm −1 (Figure 5), confirming that the Fe 2+ (CO) CN − 2 moiety represents a tightly interacting unit that may facilitate energy redistribution at the [NiFe] active site. In this respect, normal mode #35 plays a key role since this highly mixed vibrational degree of freedom ( Figure S3) exhibits nonnegligible coupling (ca. 1 − 3 cm −1 ) with νCO, νCN, and all other δ/νFe-CO/CN modes ( Figure 5). We therefore propose that Fe-CO/CN vibrations, especially mode #35, are responsible for ultrafast energy transfer between the weakly coupled CO and CN − ligands (vide supra). Remarkably, mode #35 is also coupled to both Ni-S metal ligand vibrations and S-C side chain modes of the cysteine ligands (ca. 1 cm −1 ; Figures S4 and S5). This finding supports our previous proposal that energy transfer from the [NiFe] site towards the protein environment could be accomplished via a distinct route involving the cysteine ligands [55]. This could also explain the observation of level-independent vibrational lifetimes, which conflicts with the behavior expected for a (near) harmonic oscillator linearly coupled to a harmonic bath, i.e., undirected energy dissipation [78]. pled to both Ni-S metal ligand vibrations and S-C side chain modes of the cysteine ligands (ca. 1 cm −1 ; Figures S4 and S5). This finding supports our previous proposal that energy transfer from the [NiFe] site towards the protein environment could be accomplished via a distinct route involving the cysteine ligands [55]. This could also explain the observation of level-independent vibrational lifetimes, which conflicts with the behavior expected for a (near) harmonic oscillator linearly coupled to a harmonic bath, i.e., undirected energy dissipation [78].

General Aspects
Quantum mechanical calculations were performed in Gaussian16 (revision A.03) [79] on the density functional level of theory. The BP86 functional was employed in all cases [80,81], using the triple-zeta valence basis set including polarization functions, def2-tzvp [82], for the two metal ions and the 6-31 g(d) [83,84] split-valence double-zeta basis set for all other atoms. In selected cases, highlighted in the main text, the def2-tzvp basis set was used for all atoms. Further computational choices were made with the high-accuracy demands of anharmonic frequency analyses in mind (vide infra). A 'superfine' pruned integration grid was used for computing two-electron integrals and their derivatives throughout all calculations (175 / 250 radial shells for first-row / higher-row atoms and 974 angular points per shell in both cases). The convergence criteria for the self-consistent-field procedure (root-mean-square and maximum change of the

General Aspects
Quantum mechanical calculations were performed in Gaussian16 (revision A.03) [79] on the density functional level of theory. The BP86 functional was employed in all cases [80,81], using the triple-zeta valence basis set including polarization functions, def2tzvp [82], for the two metal ions and the 6-31 g(d) [83,84] split-valence double-zeta basis set for all other atoms. In selected cases, highlighted in the main text, the def2-tzvp basis set was used for all atoms. Further computational choices were made with the high-accuracy demands of anharmonic frequency analyses in mind (vide infra). A 'superfine' pruned integration grid was used for computing two-electron integrals and their derivatives throughout all calculations (175 / 250 radial shells for first-row / higher-row atoms and 974 angular points per shell in both cases). The convergence criteria for the self-consistent-field procedure (root-mean-square and maximum change of the density matrix) were set to 10 −9 and 10 −7 , respectively. According to the Gaussian16 documentation, this corresponds to an energy change of 10 −18 Ha. Calculations were performed on eight parallelized processors with a total of 10 GB of allocated memory.

Computational Models and Geometry Optimization
Starting geometries for all computational models were constructed from a previously established first-coordination-sphere model of the [NiFe] active site [24,56,60,68]. For reasons described in the Results and Discussion section, the size of the models was slightly reduced compared to previous studies by modelling cysteine sidechains as methanethiolate rather than ethanethiolate. All geometries were fully optimized using an analytic Hessian matrix in the final optimization steps and 'very tight' convergence criteria (maximum force = 2.0 × 10 −6 , root − mean − square force = 1.0 × 10 −6 , maximum displacement = 6.0 × 10 −6 , and root − mean − square displacement = 4.0 × 10 −6 ; all in atomic units).

Vibrational Analysis
Harmonic vibrational frequencies were calculated by diagonalization of the massweighted Hessian matrix. Anharmonic transition energies of fundamental, overtone, and combination transitions were calculated using an automatized implementation of the GVPT2 approach [61][62][63][64][65]. Cubic and semi-diagonal quartic force constants for these analyses were obtained by numerical derivation of analytic Hessian matrices with respect to the normal coordinates. Individual low-frequency large-amplitude modes that were found to be problematic in the treatment of resonances were generally 'inactivated', i.e., they were not included in the final GVPT2 analysis but associated cubic and quartic force constants were retained as long as there were no indications for detrimental influences on the results [62]. In calculations with reduced dimensionality (only selected vibrations included in the harmonic analysis), all other vibrational modes were 'frozen', i.e., associated cubic and quartic force constants were discarded altogether.
Harmonic and anharmonic transition frequencies, ω H and ω 0→1 , were directly extracted from the Gaussian output. The single mode sequence transition frequency ω 1→2 was calculated by subtracting ω 0→1 from the overtone frequency ω 0→2 . Mixed-mode sequence transition frequencies ω a0→ab were calculated accordingly by subtracting ω a 0→1 = ω 00→a0 from the combination transition frequency ω 00→ab . All other quantities were calculated as described in the Results and Discussion section.

Simulation of IR Spectra
To allow direct visual comparison of calculated and experimentally derived transition frequencies (unaffected by technical details of the experiment or other observables), simulated IR spectra are shown throughout the manuscript. Harmonic and anharmonic IR absorption spectra were simulated based on calculated or experimentally derived fundamental transition frequencies and calculated harmonic transition probabilities (IR intensities). Harmonic transition probabilities were chosen to facilitate direct comparability across all spectra and ensure consistency with signal intensities in the 2D-IR data (vide infra). Individual transitions were modelled as Lorentzian lines with a full width at half maximum of 5.3 cm −1 . Two-dimensional-IR spectra were simulated from third-order response functions as described in reference [59]. Briefly, rephasing and non-rephasing diagrams for the ZZZZ polarisation were constructed from anharmonic frequencies of fundamental, overtone, and combination transitions. The dipole terms were obtained from harmonic transition dipole moments for fundamental transitions (computed from the dipole strengths reported in the Gaussian log files), and the harmonic approximation was applied to calculate transition dipole moments associated with sequence transitions between singly and doubly excited states. Transition dipole moments of both fundamental and sequence transitions of all three ligand stretching modes were assumed to be orthogonal, and dipole strengths of the cyanide stretch modes were multiplied by a factor of 10 for better visibility in the contour plots. Purely absorptive 2D-IR spectra were obtained as the real part of summed rephasing and non-rephasing diagrams after two-dimensional Fourier transform of the corresponding time domain data. Spectra shown in the manuscript were simulated for a pump-probe waiting time of T W = 0 ps and an empirical dephasing time of T 2 = 2 ps. The latter value was chosen to obtain spectra with well recognizable signals rather than to reproduce experimental lineshapes.

Scalar Normal Mode Projection
To quantify the partial localization of symmetric and antisymmetric cyanide stretch modes to individual cyanide bond-stretch coordinates, we calculated the magnitude of scalar products of each one of these normal mode eigenvectors (unit vectors by definition) and normal mode vectors representing idealized single-bond-stretch vibrations for each one of the two cyanide ligands (all in Cartesian displacement coordinates). The latter were set up as unit vectors that describe the displacement of carbon and nitrogen atoms of the considered cyanide ligand from their equilibrium position-in opposite directions along the bond axis and with relative amplitudes that ensure that the center of mass is retained. Displacement amplitudes of all other atoms were set to zero for these idealized bond-stretch normal mode vectors.

Conclusions and Outlook
We have presented an accurate and computationally feasible approach for analyzing molecular quantities encoded in PP and 2D-IR spectra of [NiFe] hydrogenase and other cyanido metal carbonyls, in particular intramode and intermode anharmonicities. Our calculations demonstrate that the characteristic pattern of anharmonicities, as observed experimentally for this class of enzymes, can be reproduced by a first-coordination-sphere model of the hydrogenase active site, which highlights a decisive role of the [NiFe] core in shaping the active-site PES and associated vibrational spectroscopic properties. This finding is relevant for understanding functional determinants of hydrogenase catalysis, e.g., the role of the Fe 2+ (CO) CN − 2 moiety, and the information that can be inferred from advanced vibrational spectra. The fact that these subtle details of the PES are surprisingly well modelled without including the protein environment indicates that the intrinsic chemical properties of the Fe 2+ (CO) CN − 2 are largely retained inside the active-site cavity. This implies that the protein may be evolutionarily optimized to preserve these intrinsic properties, and hydrogen-bonding interactions of the CN − ligands seem to anchor the active site inside this cavity without considerably affecting its chemical properties, in line with previous suggestions [70].
Analyzing CO and CN stretch vibrations in detail, we found that a proper understanding of nonlinear IR spectra of [NiFe] hydrogenase requires a decomposition of experimental anharmonicities with respect to different contributing molecular coordinates. Utilizing reduced-dimensionality schemes [62] for the analysis of partial intramode anharmonicities, we have demonstrated that computational spectroscopy can yield detailed insights into CO and CN − bond properties that are encoded in but not directly accessible from 2D-IR experiments. On the one hand, we have shown that the experimental CO stretch anharmonicity reflects both CO and Fe-CO bond properties; on the other hand, our studies demonstrate how computational analysis of the strongly interacting CN stretch modes allows insights into active-site symmetry and individual CN − bond properties. In this context, taking vibrational resonances into account was found to be crucial for the reproduction and understanding of the experimental data.
Both CO and CN stretch modes were furthermore found to tightly interact with a network of metal-ligand Fe-CO/CN vibrations that are also coupled to other degrees of freedom, e.g., δ/νNi-S. This feature could explain ultrafast energy transfer between weakly coupled CO and CN − ligands and level-independent rates in heat dissipation towards the protein matrix [55]. These unusual features resemble observations made for carboxyheme species [85][86][87][88][89][90][91], and they may point towards highly directional energy transfer that could be relevant for aspects like signaling and/or non-classical catalysis in hydrogenase. In this respect, we would also like to highlight that the δ/νFe-CO/CN and δ/νNi-S vibrations, predicted to be coupled to IR-active νCO and νCN modes, can be probed with resonance Raman spectroscopy [11,12,15,68]. Thus, our calculations demonstrate that further experimental insights into these interactions and energy transfer phenomena may be obtained by nonlinear experiments that include both infrared and Raman transitions [92].
Future studies will extend our current approach in two respects. On the one hand, we will move on to larger computational models, higher accuracy, and more efficient computational strategies that will allow detailed insights, e.g., into the influence of the protein matrix on the analyzed observables. In concert with experimental studies on synthetic models of the [NiFe] active site, this will allow more quantitative insights and reveal how the spectroscopic properties of different intermediates are determined by the interplay of structural and environmental factors. On the other hand, we also aim to model the shape and time-evolution of 2D-IR signals to reveal their relation to functional dynamics and unusual energy transfer phenomena in more detail. Finally, the introduced concepts will be generalized and transferred to other molecular targets of interest, e.g., [FeFe] hydrogenase. This way, we plan to firmly establish computationally assisted 2D-IR spectroscopy as a powerful tool for the analysis of complex metalloenzymes and bioinorganic catalysis.

Data Availability Statement:
The data presented in this study are available in the article and the Supplementary ,aterial.