Diffusion and Spectroscopy of H₂ in Myoglobin

Abstract

The diffusional dynamics and vibrational spectroscopy of molecular hydrogen (H₂) in myoglobin (Mb) is characterized. Hydrogen has been implicated in a number of physiologically relevant processes, including cellular aging or inflammation. Here, the internal diffusion through the protein matrix was characterized, and the vibrational spectroscopy was investigated using conventional empirical energy functions and improved models able to describe higher-order electrostatic moments of the ligand. Depending on the energy function used, H₂ can occupy the same internal defects as already found for Xe or CO (Xe1 to Xe4 and B-state). Furthermore, four additional sites were found, some of which had been discovered in earlier simulation studies. Simulations using a model based on a Morse oscillator and distributed charges to correctly describe the molecular quadrupole moment of H₂ indicate that the vibrational spectroscopy of the ligand depends on the docking site it occupies. This is consistent with the findings for CO in Mb from experiments and simulations. For H₂, the maxima of the absorption spectra cover ∼20 cm⁻¹ which are indicative of a pronounced Stark effect of the surrounding protein matrix on the vibrational spectroscopy of the ligand. Electronic structure calculations show that H₂ forms a stable complex with the heme iron (stabilized by ∼−12 kcal/mol), but splitting of H₂ is unlikely due to a high activation energy (∼50 kcal/mol).

1. Introduction

The interaction of myoglobin (Mb) with small molecules is of profound interest from a physiological perspective. Myoglobin is structurally related to one of the two $(\alpha ,\beta )$ subunits of tetrameric hemoglobin (Hb), both of which bind, store, and transport diatomic ligands (O₂, NO, CO). Due to their high biological relevance, both proteins interacting with all three diatomics have been intensely studied experimentally [1,2,3,4,5,6,7] and computationally [8,9,10,11,12].

The tertiary structure of Mb is characterized by a number of internal cavities. These were first mapped out by pressurizing the protein with xenon gas [13]. The physiological role these packing defects may play indicates that blocking sites “B” and/or “Xe4” have important consequences for the overall rates of oxygen binding to Mb [14]. The dynamics of ligands between these sites were also investigated using computer simulations [15]. Interestingly, proteins other than Mb also display such cavities, including neuroglobin or truncated hemoglobin [16,17,18,19]. Hence, such packing defects may be functionally relevant, and characterizing them and the dynamics between them for various ligands is of general and fundamental interest.

Molecular hydrogen, H₂, has been reported to play physiologically relevant roles in cell protection by reducing hydroxyl radicals [20], has shown therapeutic effects in carcinoma after hyperbaric hydrogen therapy in mice [21], and has been used as a medical therapeutic gas to treat brain disorders [22]. Furthermore, through its antioxidative effect, H₂ maintains genomic stability, mitigates cellular aging, and influences histone modification, telomere maintenance, and proteostasis. In addition, the diatomic may prevent inflammation and regulate the nutrient-sensing mTOR system, autophagy, apoptosis, and mitochondria, which are all factors related to aging [23]. Hence, H₂ is expected to play various roles in the human body.

Myoglobin and other heme-based proteins are known to interact with small molecules, in particular diatomics. Gas inhalation as disease therapy has been investigated and heme-containing proteins—in particular cytochrome c oxidase—have been found to be primary targets for small molecules such as O₂, NO, CO, or H₂S [24,25]. Furthermore, H₂ did not reduce the oxidized heme in cytochrome c which implicated that the primary target for H₂ appeared to be different from cytochrome c oxidase [20]. Interestingly, combined therapy with H₂ and CO demonstrated enhanced therapeutic effects [26]. This is akin to hyperbaric hydrogen treatment, which used 2.5% O₂ combined with 97.5% H₂ and showed regression of skin tumors in mice [21].

Given that (i) H₂ has been found to be physiologically relevant in general (see above) and (ii) myoglobin is known to favorably interact with a wide range of diatomics solvated in blood—which also applies to molecular hydrogen—characterizing the interaction between Mb and gaseous H₂ appears to be desirable even in the absence of an established relationship between the two. On the other hand, it is pointed out that the activation of heme oxygenase-1 by H₂ has been described in the literature [27]. Hence, first studying the interaction between H₂ and the “hydrogen atom of biology”—myoglobin—is a first step for obtaining molecular-level insights. Due to its small size and the fact that it is electrically neutral, it can be expected to diffuse easily into and within the protein. On the other hand, H₂ has an appreciable electric quadrupole moment [28]. Hence, from a spectroscopic perspective, it may behave in a similar fashion to CO, which is also electrically neutral, with a rather small permanent dipole moment but a sizable molecular quadrupole [29].

Vibrational spectroscopy is a powerful means to provide structural information for proteins. One example is the positioning of CO in the B-state and the Xe4 pocket of Mb, which has been characterized through experiments and simulations [5,30,31,32]. Ligand migration and spectroscopic features were also analyzed for CO in neuroglobin using experiments and simulations [33]. Computationally, vibrational spectroscopy from MD simulations has become a standard tool, specifically in the realm of empirical and physics-based energy functions [34] but also from ab initio MD simulations [35].

The present work characterizes the ligand diffusion and vibrational spectroscopy of H₂ in Mb from atomistic simulations. Furthermore, the binding and chemistry of H₂ attached to the heme iron is considered. First, the computational methods are presented. This is followed by characterizing the interaction between H₂ and the heme group, the structural dynamics of H₂ in the protein, and the vibrational spectroscopy in particular internal pockets. Finally, conclusions are drawn.

2. Material and Methods

2.1. Molecular Dynamics Simulations and Analysis

Molecular dynamics (MD) simulations of myoglobin and one or five H₂ molecules were performed using the CHARMM [36,37] molecular simulation package starting from the X-ray structure 1 MBC [38] prepared as described in a previous work [39]. H₂ was initially inserted within the known pockets for xenon atoms in myoglobin [12,13,15,32]. The simulations were carried out in a cubic box of size ${\left(61.3\right)}^3$ ų with explicit TIP3P water [40] and including 29 K⁺ and Cl⁻ ions. Bond lengths involving H-atoms were constrained using SHAKE [41] except for the H₂ molecule itself. A cutoff of 12 Å with switching at 10 Å was used for non-bonded interactions [42]. The system was initially heated in the $NVT$ simulation to $303.15$ K for 40 ps, followed by 40 ps equilibration simulation in the $NpT$ ensemble at $p=1$ bar using the leap-frog integrator and a Hoover thermostat [43]. Subsequently, production runs were performed for 5 ns again using the $NpT$ ensemble. This has been found sufficient to converge the vibrational spectroscopy of small freely diffusing molecules with high-frequency vibrations, in particular diatomics [32,44]. This differs from covalently bound spectroscopic probes, attached to the protein framework, for which not only the vibrational mode itself needs to be sampled sufficiently but also the protein conformational dynamics to which the small molecule is coupled [45,46,47].

To identify internal localization sites for the H₂ ligand, candidate sequences of simulation trajectories were selected. The H₂ ligand was considered occupying a pocket if the H₂ center of mass remained within 5 Å of the pocket center for a time ${\tau }_{\mathrm{dwell}}$. The dwell time chosen was ${\tau }_{\mathrm{dwell}}=20$ ps but is largely arbitrary. Iterating over all candidate sequences, the distances between the H₂ center of mass and all C_α atoms were computed for the trajectories. Within the search radius of 10 Å, a set of C_α atoms in amino acid residues was selected, so that their geometric center overlaps best with the average center of mass of H₂ moving inside the pocket. The procedure was repeated with additional simulations and candidate sequences, if no clear pocket could be identified. Pockets Xe1 to Xe4, the B-state, and pockets 6 to 9 identified by following this procedure are shown in Figure 1 and their definition by the residue number and name are documented in Table S1.

The line shape $I\left(\omega \right)$ of the power spectra for H₂ in myoglobin are obtained via the Fourier transform of the distance–distance autocorrelation function from the H₂ separation $r\left(t\right)$

$$I\left(\omega \right)n\left(\omega \right)\propto Q\left(\omega \right)·\mathrm{Im}{\int }_0^{\infty }dt{e}^{i\omega t}〈r\left(t\right)·r\left(0\right)〉,$$

(1)

where r is the H₂ bond length. A quantum correction factor $Q\left(\omega \right)=tanh(\beta \hslash \omega /2)$ was applied to the results of the Fourier transform [48]. It should be noted that Equation (1) yields the infrared absorption spectrum, if $r\left(t\right)$ is replaced by the dipole moment $\overrightarrow{\mu }\left(t\right)$ (which is a rescaled atom separation for a diatomic), or the Raman spectrum, if the polarizability tensor $\alpha$ is used instead of $r\left(t\right)$. A direct comparison between the power and infrared spectrum has been given, for example, for protonated water dimer in the gas phase, which was used to assign modes and characterize coupling between internal degrees of freedom [49]. The procedure used here (“power spectra”) yields H₂ vibrational signals at the respective frequencies but not the absolute intensities as would be observed in experimentally measured vibrational spectra. This is akin to using velocity autocorrelation functions to determine the vibrational density of states, which yields the positions of fundamentals, overtone, and combination bands but not the intensities of transition frequencies [50,51].

2.2. The Energy Function

The MD simulations were carried out using the all-atom force field CHARMM36 (CGenFF) [52], the corresponding TIP3P water model [40], and the Lennard-Jones (LJ) parameter to describe the non-bonded van-der-Waals interaction for the K⁺ and Cl⁻ ions [53]. For H₂, the bonded interaction was either a harmonic potential with a force constant of $k=350$ kcal/mol/Ų and $r_e=0.7414$ Å or a Morse potential fitted to energies determined at the CCSD(T)/aug-cc-pVQZ level of theory using the Gaussian16 program code [54]. For this, the H₂ bond was scanned between 0.5 and 1.5 Å in steps of 0.01 Å, and the energies were represented as $V\left(r\right)=D_e·{(1-{\mathrm{e}}^{-\beta ·(r-r_e)})}^2$. The fit yielded $D_e=111.76$ kcal/mol, $r_e=0.7477$ Å, and $\beta =1.9487$ Å⁻¹. The Lennard-Jones parameters for H₂ were those from the literature [55], which were fitted for accurate H₂ gas and the interface biomolecules. For the H₂ electrostatics, a minimal distributed charge model was developed as described below.

To validate the H₂ stretch potential, the anharmonic frequency was determined from solving the nuclear Schrödinger equation based on a discrete variable representation. The fundamental transition was found at $\nu =4202.4$ cm⁻¹ compared with the experimentally reported value of 4161.1 cm⁻¹ for the rotationless transition [56]. For the harmonic potential, the frequency is at $\omega =4047.0$ cm⁻¹.

2.3. Electronic Structure Calculations

For the potential energy surface (PES), for H₂ interacting with the heme unit, electronic structure calculations using the ORCA program were carried out [57]. The PES was scanned along the H₂ bond r and the z-direction between the Fe atom and the H₂ center of mass. A mixed quantum mechanics/molecular mechanics (QM/MM) approach was adopted to include electrostatic interactions between the surrounding protein and the His-Heme-H₂ subsystem with the iron atom in its Fe(II) oxidation state. For this, all protein atoms were assigned their CGenFF charges, and the His-Heme-H₂ subsystem was treated at the rPBE/def2-TZVP level of theory, including D4 dispersion corrections [58,59,60]. First, the structure of H₂ was optimized with otherwise constrained atom positions. Then, scans along directions z and r were carried out, see Figure S1. The z-direction was set as the normalized vector from the Fe to the center of mass of H₂ whereas, the r-direction was defined between both H atoms, and it was corrected to be orthogonal to the z-vector, also considering the optimized position and rotation of H₂ on Fe in myoglobin. The center of mass of H₂ was shifted along the z-direction to match the distances to the Fe atom from 0.28 Å to 3.18 Å in 0.1 Å steps. From the shifted center of mass, the H₂ atoms were set apart from 0 to 3.2 Å along the r-direction in 0.2 Å steps, conserving the center of mass.

2.4. MDCM Model for H₂

With a standard force field-based energy function, H₂ only interacts through van-der-Waals interactions with its environment, due to its neutrality and vanishing molecular dipole moment. The first non-vanishing permanent moment of H₂ is the molecular quadrupole moment. To capture this, the electrostatic potential (ESP) of the neutrally charged H₂ molecule is reproduced by the minimal distributed charge model (MDCM) [61] using 3 off-centered charges per hydrogen atom. The off-centered charges are restricted along the bond axis and symmetric to the horizontal mirror plane perpendicular to the H₂ bond. The reference ESP is computed at the CCSD(T)/aug-cc-pVQZ level of theory using Gaussian16 at the H₂ equilibrium conformation at the CCSD/Aug-cc-pVQZ level of theory ($r_e=0.7424$ Å) and represented in a cube file format with Gaussian16’s default coarse grid resolution [54].

The off-centered charge displacements and amplitudes are optimized to best fit the ESP grid points in the range of 1.44 ($1.2·r_{\mathrm{vdW}}$) to 2.64 Å ($2.2·r_{\mathrm{vdW}}$) around the closest hydrogen atom, which are related to the van-der-Waals radius of the hydrogen atom $r_{\mathrm{vdW}}=1.2$ Å [62]. The MDCM reproduces the ESP grid points within the range with a root mean square error (RMSE) of $0.73$ kcal/mol per grid point. In comparison, the RMSE within the same range for a point charge model with charges set to zero is $2.31$ kcal/mol. Contour plots of the reference ESP ($V_{\mathrm{ref}}$) on grid points at distances larger than $1.44$ Å from the closest hydrogen atom and the deviation from the model ESP ($V_{\mathrm{MDCM}}$) with $\Delta \mathrm{ESP}=V_{\mathrm{MDCM}}-V_{\mathrm{ref}}$ are shown in Figure 2A and B, respectively.

The computed components of the molecular quadrupole tensor $\mathit{Q}$ of H₂ are $\{Q_{xx},Q_{yy},Q_{zz}\}=\{-0.22,-0.22,0.45\}$ DÅ, close to the measured experimental results of $\{-0.26,-0.26,0.52\}$ DÅ [28]. The quadrupole moment prediction of the fitted MDCM model for H₂ yields $\{0,0,0.61\}$ DÅ. Both quadrupole moments $Q_{xx}$ and $Q_{yy}$ are zero, as the distributed charges are only displaced along the H₂ bond axis, but $Q_{zz}$ is about 35% larger than the reference one but close to the experimental result.

3. Results

First, the interaction between H₂ and the heme unit is considered. Next, the structural dynamics and vibrational spectroscopy of H₂ in the protein are analyzed from MD simulations.

3.1. H₂ Interaction with Heme

The interaction between H₂ and the heme unit is shown in Figure 3. With respect to H₂ outside the protein, the Fe(II)-H₂ bound state is stabilized by $\sim -$11.5 kcal/mol. This is considerably weaker than the interaction between heme and the physiologically relevant ligands -O₂ and -NO and poisonous -CO and -CN⁻. For -O₂ and -CO, the experimentally determined standard enthalpies of formation with myoglobin are $-18.1\pm 0.4$ kcal/mol and $-21.4\pm 0.3$ kcal/mol respectively [63]. No direct experimental measurements exist for -NO and -CN⁻, but density functional theory calculations indicate comparable or stronger interactions depending on the Fe-oxidation state: $-23.0$ kcal/mol for Fe(II)-NO, more than $-33.0$ kcal/mol for Fe(III)-NO, and stronger than $-50.0$ kcal/mol for binding of -CN⁻ [44,64]. Coordination of H₂ with Fe leads to pronounced frequency shifts by 1000 cm⁻¹ or more to the red [65] but were not further considered in the present work, which focuses on the diffusion and spectroscopy of free H₂ in Mb.

The capability of the heme-Fe(II) to break the H₂ bond was also investigated. The H-Fe-H arrangement was found to be a faint minimum, $13.9$ kcal/mol above the dissociation limit or $25.4$ kcal/mol above the global minimum. The transition state separating the two minima is $52.6$ kcal/mol above the Fe-H₂ state, i.e., $41.1$ kcal/mol above the dissociation limit. Hence, H₂ binds reversibly to heme iron and no “chemistry” is expected to take place.

3.2. Structural Dynamics and H₂ Diffusion

Next, the diffusional dynamics of H₂ within Mb were considered. MD simulations were carried out using two energy functions. The first was the conventional CGenFF energy function, and for the second energy function, the H₂ molecule was described as a Morse oscillator and MDCM electrostatics.

Using the CGenFF setup and MD simulation with five H₂ in Mb, it was found that H₂ can localize in at least nine different locations within Mb. The first four pockets (Xe1 to Xe4) are the those found for xenon in myoglobin [13,15], together with the B-state, which was spectroscopically characterized [5,14,30,31,32]. Pockets 6 to 9 were detected from extended MD simulations. From visual comparison, the newly determined ones were also observed as CO cavities in previous publications [12,66], see e.g., Figure 1 in Ref [66].

Time series for the separation of the center of mass of H₂ to each of the nine pocket centers are reported in Figure 4C,D. The pocket centers were determined from the procedure described in Section 2 MD simulations with one H₂ in Mb. It is found that the H₂ molecule readily migrates between the different pockets, see Figure 4A. For example, the simulation using the CGenFF energy function finds H₂ visiting six out of the nine pockets during 600 ps. Similarly, using the MDCM model for the electrostatic (panel D) also leads to diffusion, but only four different pockets are visited within the same time frame. This indicates that the interaction between H₂ and the protein environment is stronger due to the modified electrostatics. The unassigned spatial location between 150 and 300 ps is likely to be a sub-pocket of Xe4, which was also found for CO diffusing through Mb (called Xe4(2)) [67,68].

3.3. H₂ Vibrational Spectra

Next, the vibrational spectroscopy of H₂ within the protein was analyzed by computing the power spectra from the H₂ bond distance. It is of interest to assess whether the electrostatic interaction between the protein environment and H₂ leads to pocket-specific spectra once MDCM as the electrostatic model for the diatomic is used. Two types of simulations, unconstrained and constrained, were carried out. The unconstrained MD simulations were initialized with H₂ close to the heme-Fe, above the porphyrin plane in the distal site of Mb (B-state). Second, to investigate the impact of the pocket positions on the vibrational spectroscopy of H₂, constrained simulations were performed in which H₂ was weakly harmonically constrained to each center of mass of one of the nine pockets to obtain pocket-specific spectra.

Figure 5 reports the vibrational spectrum from the free dynamics, which featured a single H₂ in Mb. Only the times during which H₂ was within one of the specified pockets were analyzed. Because the residence times of the ligand in each of the pockets differ, the sampling times also differ. In Figure 5A (CGenFF with harmonic H₂), the power spectra are all broad and feature a single maximum, except for that associated with pocket 7. There is a slight shift of the mean peak positions of each spectrum, with that of pocket Xe2 most shifted to the blue by $+16.2$ cm⁻¹ (maximum at $4078.6$ cm⁻¹) and that of pocket 7 least shifted to the blue by $+9.1$ cm⁻¹ (at $4071.5$ cm⁻¹), away from the harmonic gas phase spectrum. The blue shift also indicates that the H₂ bond is slightly strengthened, due to repulsive van-der-Waals interactions between the unbound ligand and the surrounding protein.

For simulations using the Morse potential for the H₂ bond and the MDCM model to include the molecular quadrupole moment, the power spectra are reported in Figure 5B. Here, the ligand samples pockets Xe4, B-state, and 7 but not pocket Xe2. With the refined energy function, the spectra differ in their widths both between the pockets and compared with simulations using the CGenFF energy function (see Figure 5A). For example, the spectrum associated with pocket 7 is rather narrow, whereas that for pocket Xe4 is overall broad but consists of two separate peaks. In addition, the pocket-specific spectra feature both red- and blue-shifts, which indicate slight weakening and strengthening of the H₂ bond due to favorable and unfavorable intermolecular interactions with the protein environment. For H₂ in pocket Xe4, the red shift with respect to the gas phase spectrum amounts to $-9.3$ cm⁻¹, whereas the blue shift for pocket 7 is $+10.8$ cm⁻¹.

Next, the pocket-specific spectra are analyzed, see Figure 6. In all simulations, the H₂ ligand was slightly constrained towards the center of each of the nine pockets using a harmonic constraint. This was to ensure that pocket-specific spectra could be obtained. Three types of simulations were carried out: one using the CGenFF force field, a second one using the MDCM charge model but the harmonic bond potential, and the third combining MDCM with the Morse bond for H₂.

Broadly speaking, the power spectra for H₂ in all nine pockets from simulations using the CGenFF force field are rather similar, not to say largely identical, see Figure 6A. All peaks are shifted to the blue away from the gas phase spectrum, and the peak positions cover a range of only 3 cm⁻¹ (from 4078 cm⁻¹ to 4081 cm⁻¹). Using the harmonic bond together with the MDCM model leads to similar observations, as shown in Figure 6B. However, including a Morse description for the H₂ bond leads to considerable changes, see Figure 6C. Now both blue- and red-shifts of the power spectra appear, and the maxima of the spectra cover a considerably wider range, covering 22 cm⁻¹ (4460 cm⁻¹ to 4482 cm⁻¹). This finding suggests that only the combination of an improved description of the bonded interaction (Morse) together with modeling the electrostatics (quadrupole from MDCM) provides the necessary detail to yield the expected response of H₂ to the inhomogeneous electric field in each of the protein pockets, as was also observed for CO from both simulations and experiments [8,30,31,32,69,70,71].

4. Discussion and Conclusions

In this work the interaction of H₂ with the heme-group of Mb and the diffusional dynamics of the diatomic within the protein were investigated. This was motivated by the observation that, over the past 15 years, H₂ has emerged as a physiologically interesting “small molecule” akin to the well-known diatomics O₂ or NO. The results of the present study demonstrate that H₂ can reside within Mb, can occupy the same spatial regions as those found experimentally and from computations for Xe, CO, and NO, and can even populate less well-characterized internal defects. With H₂ initially localized within the protein, the ligand escaped in approximately half of the cases using the CGenFF energy function but only one in five when using the MDCM model for H₂ on the 5 ns time scale. Examples for such trajectories of escaping H₂ using both models are shown in Figures S2 and S3. However, for a statistically significant result, a considerably larger number of trajectories needs to be run.

The magnitude of the quadrupole moment of H₂ is approximately half that of carbon monoxide (CO), for which a value of $\Theta =-9.47\times {10}^{-40}$ cm² was found experimentally [72], compared with calculations of $\Theta \sim -$2 DÅ, corresponding to $\sim -$6.7 $\times {10}^{-40}$ cm² [29]. Hence, the Stark shifts originating from the interaction between the electrical moments of the ligand and the electrical field of the surrounding protein are expected to be comparable but somewhat smaller for H₂ compared with CO [8,30,31,32,69]. As a validation of the present approach, it is noted that for CO in the B-state of Mb, the experimentally measured splitting between the two conformational substates is 10 cm⁻¹, which is accurately captured from simulations using the same techniques as those used in the present work [30,32]. Furthermore, the Stark shift of the split IR spectrum is −10 cm⁻¹ and −20 cm⁻¹ for the two experimentally measured bands, which correspond to the Fe–CO and Fe–OC orientations in the B-state [8,30,31,32]. For the B-state, the frequency shift for H₂ is ∼4 cm⁻¹ to the blue, see Figure 6C, but for other pockets, the shift can be considerably larger.

A further validation of the approach followed in the present work is afforded by considering previous vibrational spectroscopy studies of H₂ in water clathrates, both from experiments [73] and from classical MD and path integral MD simulations [74]. Similar to the present work, these simulations used a quadrupolar model for H₂, and the frequency ranges over which the H₂ vibrations are distributed from the simulations were found to agree favorably with experiments. This further validates the use of quadrupolar models for H₂ for spectroscopic studies.

Comparing the pocket-specific vibrational spectra, it is found that the H₂ vibrational spectra change as the ligand occupies different regions within the protein, see Figure 6. Because H₂ has a vanishing dipole moment, it is more likely to observe its Raman spectrum [75]. The use of Raman spectroscopy to query myoglobin is well-established [76], and the frequency range for H₂ above 4000 cm⁻¹ is well-removed from other spectroscopic signatures of the protein. Another possibility is to use difference spectroscopy between ligand-occupied and empty Mb. This is also possible from MD simulations, as has been shown recently [77].

In view of the chemical reactivity of H₂ towards heme, Figure 3 establishes the existence of two low-energy states. The first is located where the separation between the H₂ center of mass and Fe is 1.5 Å, and the distance between the H atoms is 0.81 Å. For the second low energy state, the distance between the H₂ center of mass and Fe is ∼1 Å but with a direct Fe-H distance of 1.8 Å. The first state was expected, as the distance between the H two atoms is comparable to the equilibrium separation for H₂ in the gas phase. The second state is more interesting, since the H atoms are separated by ∼2.7 Å, which suggests that the H–H bond can be broken when bound to Fe and remain stable. Nevertheless, reaching the H-Fe-H dissociated state involves a high barrier (∼50 kcal/mol) and is physiologically irrelevant.

If H₂ is found to be associated with physiological function in Mb, it may also be of interest to determine migration paths and barriers between neighboring pockets and between the protein’s interior and the outside. In addition, it will be interesting to map the connectivity of the network sampled by H₂, as has already been performed for Xenon [15]. Such studies, however, require considerably longer simulation times for convergence [68].

In conclusion, the present work finds that the vibrational spectroscopy of H₂ is sensitive to the chemical environment. H₂ as a ligand is well-tolerated within myoglobin, and given the role of H₂ for various physiological processes, it is of interest to further characterize the interaction between H₂ and myoglobin from an experimental perspective. It is hoped that the present work provides an initial stimulus for such studies.