First-principles study of molecular hydrogen binding to heme in competition with O₂, NO and CO

Abstract

Molecular hydrogen shows antioxidant activity and distinct efficacy towards vascular diseases, but the understanding of this is not yet satisfactory at the atomic level. In this work, we study the binding properties of H₂ to the heme group in relation with other diatomic molecules (DMs), including O₂, NO and CO, and their displacement reactions, using first-principles calculations. We carry out molecular modeling of the heme group, using iron-porphyrin with the imidazole ligand, i.e., FePIm, and smaller models of Fe(C_nH_n+2N₂)₂NH₃ with n = 3 and 1, and of molecular complexes of heme–DM and –H. Through analysis of optimized geometries and energetics, it is found that the order of binding strength of DMs or H to the Fe of heme is NO > O₂ > CO > H > H₂ for FePIm-based systems, while it is H > O₂ > NO > CO > H₂ for model-based systems. We calculate the activation energies for displacement reactions of H₂ and H by other DMs, revealing that the H₂ displacements occur spontaneously while the H displacements require a large amount of energy. Finally, our calculations corroborate that the rate constants increase with increasing temperature according to the Arrhenius relation.

---

1 Introduction

Recently, molecular hydrogen (or hydrogen gas, H₂) has shown great potential for treating various fatal diseases including cancer and angiopathy.^1–3 H₂ is the smallest and lightest molecule, coupled with its non-polar nature, thereby readily accessing any organ through facile penetration across cell membranes and fast diffusion into cellular organelles.^4,5 Moreover, hydrogen gas is odorless, colorless and non-toxic, and hydrogen is abundant on Earth as part of water molecules. At the early stage of therapeutic applications, hyperbaric H₂ was applied to treat skin squamous carcinoma in mice,⁶ and later a clinically viable dose of H₂ was found to significantly reduce cytotoxic reactive oxygen/nitrogen species (RONS),⁷ indicating that hydrogen may be an effective antioxidant in the clinic. The clinical effects of molecular hydrogen have been steadily reported in humans and animals over the past two decades.^8–10

In particular, various vascular diseases, such as diabetic vasculopathy, atherosclerosis and hypertension, which are mostly caused by cellular redox dysregulation, can be effectively treated by inhaling hydrogen gas or drinking H₂-rich water.^11–13 When inhaling hydrogen, the H₂ concentration in venous and arterial blood was found to increase in proportion to the applied dose.¹⁴ According to Henry’s law, a 2% H₂ concentration leads to a 15.6 μM H₂ concentration in the blood. Considering that cellular redox dysregulation occurs as a result of homeostasis depression between RONS and self-defense antioxidants, H₂ may react with RONS, typically the hydroxyl radical (˙OH) and peroxynitrite (ONOO⁻), in the presence of catalytic metal cations such as Fe²⁺ and Cu⁺.^15,16 Although some enzymes, such as catalase and superoxide dismutase, may also alleviate RONS-induced cellular injury,¹⁷ their clinical efficacy is little demonstrated because of their inability to permeate cell-membrane barriers. In addition, hydrogen has also been applied to treat cerebral ischemia by buffering the destructive effects of oxidative stress in the brain, proving its ability to cross the blood–brain barrier.^18,19

To explain the efficacy of H₂ towards vascular diseases, it was suggested that the effects of H₂ could be mediated by heme groups, such as those in catalase,^20,21 peroxidases²² and matrix metalloproteins.²³ Taking account of the fact that heme is the active site of hemoglobin (Hb), the main component of blood cells, we recently carried out the systematic study of hydrogen–protoheme (Fe-protoporphyrin with an imidazole ligand: FePIm) complexes using density functional theory (DFT) calculations.¹⁵ We suggested two types of bonding between H₂ and the Fe(II) of heme (dihydrogen and Kubas bonding), and revealed that H₂ can be dissociated with a relatively low activation barrier to form an FePIm–H complex, clarifying the key role of heme as an effective catalyst for antioxidant activity towards RONS by hydrogen gas. Some previous works also reported dihydrogen–metalloprotein complexes.^24–28 Meanwhile, other diatomic molecules (DMs), such as O₂, NO and CO, can be strongly bonded with heme, playing important biological roles in living organisms. The oxygen molecule O₂ reacts reversibly with Hb, allowing the transport and storage of dioxygen in the respiratory chain. Having a high affinity for the iron of heme, both NO and CO exhibit medicinal effects, acting as intercellular signaling molecules on one hand, while being poisons on the other hand. There is considerable previous theoretical work in the literature on the binding of DMs to heme groups^29–32 and on the geometrical and electronic properties of heme groups themselves.^33–37 Putting all accounts together, a natural question is raised: if hydrogen is bound to heme by inhaling molecular hydrogen, can DMs such as O₂, NO and CO displace the H₂ molecule or H atom bound to heme? If so, how high are the activation barriers for such displacement reactions?

The aim of this work is to answer these questions. We perform DFT calculations and ab initio molecular dynamics (AIMD) simulations of FePIm complexes with DMs to estimate the binding energies and activation energies for substitution. We organized the paper as follows. After describing the computational methods, we briefly explain the binding features of the complexes together with their molecular vibrations and electronic structures. Then we present the results regarding the energetics of substitution of O₂, NO and CO for a H₂ molecule or H atom bound to heme.

2 Computational methods

DFT calculations were carried out with the NWChem (version 6.6) package,³⁸ using the hybrid B3LYP exchange–correlation (XC) functional^39,40 and the basis set combination described below. For the systems with odd numbers of electrons, a spin-unrestricted open-shell formalism was utilized throughout. In all the DFT calculations, the 6-311G basis was used for the non-metallic elements C, O, N, and H, while the standard effective core potential (ECP) was used for Fe(II) with the standard LANL2DZ_ECP basis.⁴¹ After the geometry optimizations, we performed frequency calculations on all the molecular systems. We considered the thermal correction (TC) and the zero-point correction (ZPC) to the total DFT energy. The basis set superposition error (BSSE) was corrected for binding energy calculations using the counterpoise method.⁴² For the geometry optimizations, the DRIVER module was used based in the quasi-Newton algorithm with the convergence criteria of “tight”, being the highest accuracy for geometry optimization among the possible selections in the NWChem package. The nudged elastic band (NEB) method was applied to estimate the activation energy for substitution reactions using 10 NEB beads, followed by a zero-temperature string calculation with 10 beads.⁴³ For the AIMD simulations, the cutoff radius was set to 2.8 nm, the constant-temperature ensemble using Berendsen’s thermostat at T = 300 K was used, with a temperature relaxation time of 0.05 ps, and the time step and duration were set to 1 fs and 1 ps, respectively.

Our calculations were performed on three molecular systems to mimic the heme group, as shown in Fig. 1. The first system is Fe(II)-porphyrin-imidazole ligand (FePIm), being the largest and most realistic system in the present work (Fig. 1a), and which has been widely used to study the reactions of small molecules with heme groups in DFT calculations.^15,34,35,44 Since FePIm is too large to be used for calculation of the rate constant with direct dynamics for variational transition state theory (DIRDYVTST), the smaller models 1 and 2 were also considered. The chemical formulae of models 1 and 2 can be expressed as [Fe(C_nH_n+2N₂)₂NH₃], where n = 3 and 1, respectively. Small models similar to 1 and 2 but with different ligands have been proved previously to reproduce some features of the FePIm model.^29,35 In these models, the two ligands of vinylogous amidine C₃H₅N₂⁻ or amidine CH₃N₂⁻ are located symmetrically around the Fe centre, having the same bonding geometry and charge as the full porphyrin ring. They are distinct from the FePIm system in their conjugation properties because there is no cyclic ring. The larger model1 is more similar to the FePIm system than model2, due to it having the same number of C atoms between the N atoms and thus a more similar bond angle to the porphyrin of FePIm. The ammonia molecule in the two models, having a small size, high symmetry and single lone-pair donated to Fe, was selected to mimic the imidazole group.

Fig. 1 Molecular structures of (a) Fe-porphyrin with the imidazole ligand (FePIm), (b) model1, Fe(C₃H₅N₂)₂NH₃, and (c) model2, Fe(CH₃N₂)₂NH₃.

The chemical reaction rates were calculated for only model2 because of the quite heavy computational load for the FePIm system and model1. The calculations were carried out by connecting the NWChem and Polyrate packages.^45–47 Using the optimized geometries of reactants, products and saddle points determined by the NEB calculations, we applied the DIRDYVTST module included in the NWChem package at B3LYP/6-31G and LANL2DZ level to carry out the electronic structure calculations. We selected the Page–McIver corrected local quadratic approximation (CLQA) for the integration method to follow the reaction path. The resultant energies, gradients, and Hessians were merged into file30, which is the input file for Polyrate. In using the Polyrate package, we selected reaction-path variational transition-state theory (RP-VTST),⁴⁶ in which the canonical variational theory (CVT) forward rate constant was obtained by minimizing the generalized transition (GT) state theory rate constant with respect to the reaction coordinate s as follows:

where s is defined as the distance between two pivot points along the minimum energy path (MEP). In Polyrate, the reverse reaction rates and equilibrium rate constants were also calculated.

3 Results and discussion

3.1 Geometrical and binding features

The aim of this work is to study displacement reactions of H₂ or H bound to heme with DMs such as O₂, NO and CO. Therefore, the molecular systems to be considered were heme–DM and –H, where heme was modeled using FePIm, model1 and model2, and the DMs were H₂, O₂, NO and CO, resulting in 15 molecular systems in total. To begin with, the geometries of these 15 molecular systems were optimized at the B3LYP/6-311G and LANL2DZ_ECP level. Fig. 2 shows their optimized geometries, together with the frontier molecular orbitals including the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). As already clarified both in experiment⁴⁸ and theory,³¹ the DMs O₂, NO and CO were bound to ferrous iron, forming a new chemical bond. As such, molecular hydrogen and the hydrogen atom were also bound to the iron of heme, which has been confirmed in our previous work.¹⁵ Table 1 lists the typical bond lengths and bond angles, mostly around the Fe atom, for all the molecular systems.

Fig. 2 Optimized geometries of heme–DM and –H compounds, where heme is modeled by FePIm, model1 and model2, and the DMs are H₂, O₂, NO and CO. The highest occupied molecular orbitals (HOMOs) and the lowest unoccupied molecular orbitals (LUMOs) are presented with the MO levels in eV units. The DFT method with B3LYP/6-311G and LANL2DZ_ECP was utilized.

Table 1 Typical bond lengths d (unit: Å) and bond angles α (unit: degrees), and binding energies E_b (unit: kcal mol⁻¹) of heme–H₂, –H, –O₂, –NO and –CO compounds, where heme is modeled with FePIm, model1 and model2. Geometries were optimized with the B3LYP/6-311G and LANL2DZ_ECP methods. N_p represents the nitrogen atom of porphyrin or heme, and N_O for the cases of NO indicates the nitrogen atom of the NO molecule

		FePIm		Model1	Model2
		This	Prev.
a PBE/DZP DFT supercell calculation.^15b Experiment.^48c B3LYP/lacv3p**+ DFT cluster calculation.^31d Experiment.^49e B3LYP/CCSD(T) DFT cluster calculations.^29f B3LYP/CCSD(T) DFT cluster calculations.^35
H2	dFe–H	1.952	1.71^a	1.839	1.827
	dH–H	0.755	0.73^a	0.754	0.767
	dFe–Np	2.023	2.03^a	1.978	2.064
	αFe–H–H	78.8	77.7^a	79.1	77.9
	Eb	4.8		5.8	8.2
H	dFe–H	1.694	1.54^a	1.500	1.524
	dFe–Np	2.014	2.02^a	1.938	1.999
	Eb	−10.6		−44.9	−44.1
O2	dFe–O	1.794	1.81^b, 1.89^c	1.784	1.782
	dO–O	1.317	1.24^b, 1.35^c	1.328	1.334
	dFe–Np	2.019	2.02^b, 2.09^c	1.970	2.015
	αFe–O–O	122.3	122^b, 118.1^c	123.0	121.8
	Eb	−27.3		−32.6	−38.4
NO	dFe–NO	1.821	1.82^c	1.818	1.815
	dNO–O	1.203	1.20^c	1.217	1.225
	dFe–Np	2.026	2.06^c	1.965	2.033
	αFe–NO–O	141.1	142^c	136.4	139.8
	Eb	−28.2	−25^d, −20^e	−24.9	−35.6
CO	dFe–C	1.810	1.80^c	1.786	1.792
	dC–O	1.163	1.17^c	1.169	1.168
	dFe–Np	2.027	2.05^c	1.986	2.056
	αFe–C–O	180.0	179.9^c	180.0	180.0
	Eb	−16.7	−23.3^f	−20.2	−31.3

---

The binding geometries of O₂, NO and CO molecules to heme were confirmed to reproduce those clarified previously with DFT calculations.^29,31,44 That is, these DMs were bound to heme in the asymmetric mode by forming a chemical bond between one atom of the DM and the Fe of heme, such as Fe–O, Fe–N and Fe–C, and the other atom of the DM was not bonded, resulting in the formation of a bond angle between Fe and DM. The typical bond lengths and bond angles around Fe, calculated in the present work, were in overall agreement with the available experimental and theoretical data, as shown in Table 1. In contrast, the H₂ molecule preferred the symmetric mode to the asymmetric mode when binding to the heme, forming two Fe–H bonds, as found in our previous work with a different DFT method.¹⁵ In the heme–H compounds, which could be formed by dissociation of an H₂ molecule adsorbed on heme with a relatively low activation barrier, to act as an antioxidation catalyst, the hydrogen atom was bonded to Fe in a vertical line to the porphyrin ring.¹⁵ When comparing the different heme groups, the newly formed bond lengths of Fe–H, Fe–O, Fe–N and Fe–C were slightly decreased, while the intramolecular bond lengths of H–H, O–O, N–O and C–O were increased going from FePIm to the models. These indicate that the interaction between Fe and the DM can be enhanced while the interaction between the atoms within the DM can be reduced by using smaller models for the heme group. Compared with those in FePIm, the Fe–N_p bond lengths (N_p; nitrogen atom of porphyrin or heme) were decreased in model1 and increased in model2. The Fe–DM bond angles were more or less similar among the three heme models.

To evaluate the binding strength, the binding energy E_b was calculated as E_b = E_heme–DM(H) − (E_heme + E_DM(H)), where E_b is the DFT total energy of the s species. For the heme–DM(H) compounds, the BSSE correction was considered. For all the species, ZPC and TC at T = 298.15 K were considered for the BSSE-corrected DFT total energy by performing calculations of vibrational frequencies (see Table S1†). For the three heme–H₂ models, E_b was found to be positive, indicating an endothermic reaction of heme and H₂ binding. However, the heme–H systems showed negative values of E_b, namely −10.6, −44.9 and −44.1 kcal mol⁻¹ for FePIm, model1 and model2, respectively. This means that once the H₂ molecule is dissociated on heme, the product of the heme–H compound exists in a stable state. For other DMs of O₂, NO and CO, the values of E_b were negative for the three heme models. For the case of FePIm–NO, the calculated value of −28.2 kcal mol⁻¹ agreed well with the experimental value of −25 kcal mol⁻¹⁴⁹ and the previous theoretical value of −20 kcal mol⁻¹.²⁹ For the case of FePIm–CO, our value of E_b was 6.7 kcal mol⁻¹ higher than the previous calculation using the B3LYP/CCSD(T) method.³⁵ It should be noted that the binding structures and energies in these molecular systems depend strongly on the choice of XC functional and basis sets, thus giving a very wide spread of values, and the hybrid XC functionals like B3LYP predict much more accurate binding energies than the pure XC functionals for O₂ and CO.³⁵ Harvey and co-workers^29,35 found that the B3LYP functional underestimated the bond strength of NO to heme and thus extensive large basis set CCSD(T) calculations should be performed. Alternatively, a multi-configurational method like the complete active space SCF (CASSCF) method can be applied to get more accurate electronic structures with much heavier computational load, as performed by Radón et al.⁴⁴ In general, the models 1 and 2 showed E_b values that were larger in magnitude than those of FePIm, except for the heme–NO systems, but the differences were not so great.

The order of binding strength was NO > O₂ > CO > H > H₂ for the FePIm-based systems, while it was H > O₂ > NO > CO > H₂ for the model-based systems. Such binding characteristics can be explained by performing the analysis of the frontier molecular orbitals (FMOs) and the atomic charges. Fig. 2 shows the FMOs, such as the HOMO and LUMO, proving the quite strong covalent interaction between the Fe d orbitals and DM π* orbitals, especially for the O₂ and NO cases. In fact, the π* orbital of O₂ is closer in energy to the Fe d_z² orbital than the Fe π* orbitals, thereby obtaining more energy by bending and increasing the overlap.³¹ In previous work,²⁹ it was found via B3LYP calculations that FePIm has a triplet ground state with a +2 oxidation state of iron. When binding of O₂ to Fe(II)-heme occurs, the resultant [Fe^IIO₂] system is isoelectronic with a total of 8 electrons in the d orbitals of iron and π* orbital of O₂, leading to a singlet ground state. When binding of NO to Fe(II)-heme occurs, the interaction between three Fe d orbitals (d_xz, d_yz, and d_z²) and two NO orbitals with a total of 7 electrons induces two bonding orbitals , two antibonding orbitals and one nonbonding orbital ,⁴⁴ leading to a doublet ground state and larger bending angle (141°) than that of O₂ (122°). The energy difference between the Fe d_z² orbital and the CO π* orbital is larger than those for O₂ and NO, and stabilization is not gained by bending, leading to the almost straight geometry in the three models.³¹ The FMOs of the heme–H₂ and –H compounds looked similar to those of heme–CO compounds, implying similar binding features. From the Mulliken charge analysis (see Table S2†), it was revealed that the ligands, i.e., the side ligand (porphyrin, vinylogous amidine, and amidine) plus the Im or NH₃ ligand, accept electrons, and iron donates electrons, when forming the heme–DM or –H molecular systems. Among the adducts, H₂, H and CO act as electron donors, whereas O₂ and NO act as electron acceptors. Although the absolute amount of transferred electrons upon binding tends to decrease going from FePIm to model1 and to model2, the variation tendencies of the Mulliken charges according to the adducts look quite similar.

Using the HOMO levels (E_HOMO) and LUMO levels (E_LUMO), we further estimated the chemical stability and reactivity of the molecules, such as the ionization potential (I = −E_HOMO), electron affinity (A = −E_LUMO), electronegativity , global hardness (η = I − A), and electrophilicity index . Note that the electronegativity is the negative of the chemical potential (χ = −μ) and the electrophilicity index (ω) is a measure of the energy stabilization of the molecule caused by maximal electron flow from the environment.⁵⁰ As listed in Table 2, the heme–O₂ complexes have the lowest hardness (1.760, 1.699, and 2.332 eV for FePIm, model1 and model2, respectively) and the highest electrophilicity (5.092, 4.530, and 3.504 eV), indicating that they are the most soft and reactive. The second-lowest hardness and highest electrophilicity were found in NO-based compounds. Meanwhile, the heme–H complexes showed high values of hardness, especially in model1 and model2. Such analysis of global stability parameters is in accord with the findings from the analysis of binding features. It should be noted that the general tendencies of these parameters according to adducts in models 1 and 2 resemble the FePIm systems.

Table 2 The HOMO level (E_HOMO) assigned to the negative of the ionization potential (I = −E_HOMO), LUMO level (E_LUMO) assigned to the negative of the electron affinity (A = −E_LUMO), electronegativity , global hardness (η = I − A), and electrophilicity index in the heme–DM and –H compounds (units: eV)

Compound	EHOMO	ELUMO	χ	η	ω
FePIm	−5.039	−1.965	3.502	3.074	1.995
FePIm–H2	−4.941	−1.957	3.449	2.983	1.994
FePIm–H	−4.747	−2.075	3.411	2.672	2.178
FePIm–O2	−5.114	−3.354	4.234	1.760	5.092
FePIm–NO	−5.013	−2.172	3.592	2.842	2.271
FePIm–CO	−4.998	−2.085	3.542	2.913	2.153
Model1	−3.897	−0.212	2.054	3.684	0.573
Model1–H2	−3.959	−0.206	2.082	3.752	0.578
Model1–H	−4.680	−0.365	2.522	4.315	0.737
Model1–O2	−4.773	−3.074	3.924	1.699	4.530
Model1–NO	−4.507	−2.397	3.452	2.110	2.823
Model1–CO	−4.432	−0.471	2.451	3.961	0.759
Model2	−4.126	−0.095	2.110	4.031	0.552
Model2–H2	−4.193	0.183	2.005	4.376	0.459
Model2–H	−5.156	−0.145	2.650	5.011	0.701
Model2–O2	−5.208	−2.876	4.042	2.332	3.504
Model2–NO	−4.879	−1.912	3.396	2.968	1.943
Model2–CO	−4.748	0.012	2.368	4.760	0.589

---

3.2 Energetics for displacement reactions

Then, the work proceeded to the next stage of studying displacement reactions between H₂ or H and other DMs, namely O₂, NO and CO, binding to the heme groups. Firstly, the proper initial configurations were prepared by putting one DM on the FePIm–H₂ and –H complexes, and were optimized to produce the initial states (IS) of FePIm–H₂⋯DM and FePIm–H⋯DM complexes, where the distance between the DM and H₂ or H was controlled to not form a chemical bond. For these optimized structures, we performed AIMD simulations at T = 300 K, confirming that the temperature, starting from 0 K, arrived at ∼300 K within about 0.3 ps and then fluctuated around 300 K after that (see Fig. S1†). For the FePIm–H₂⋯DM (DM = O₂, NO) complexes, it was observed that H₂ was detached while O₂ or NO was attached to the FePIm group after about 0.5 ps. In contrast, the CO molecule was observed to be hardly attached to the FePIm group, in accordance with the tendency of the binding energies as discussed above. These results indicate that O₂ and NO molecules can easily displace a H₂ molecule bound on the heme, but CO cannot do so at room temperature. Meanwhile, it was found from the AIMD simulations of the FePIm–H⋯DM complexes that none of the molecules of O₂, NO and CO could displace the H atom attached to FePIm at room temperature. Therefore, there should be activation barriers for the displacement reactions of H by DMs.

We determined the activation barriers for these reactions by applying the NEB method, which yielded the transition state (TS) structures lying on the MEP. To do NEB runs, the proper final states (FS) were also prepared by putting a H₂ molecule or H atom on the heme–DM complexes and were optimized, resulting in heme–DM⋯H₂ or ⋯H complexes. Since the NEB method constructs an approximate MEP between the initial (reactant) and final (product) structures, the TS structures were further refined by performing a saddle-point search using the Driver module, so that the TS had only one imaginary frequency. The displacement reactions of heme–H₂ by O₂, NO and CO were first considered.

Fig. 3(a)–(c) show the calculated energy profiles and the optimized geometries corresponding to the IS, TS and FS for displacement reactions of H₂ by O₂, NO and CO on FePIm, model1 and model2. The final states, being the complexes composed of the heme group with an attached DM and detached H₂ molecule, were found to be energetically lower than the initial states of the heme–DM⋯H₂ complexes for all the cases. When replacing H₂ with O₂, the final states were 26–32 kcal mol⁻¹ lower in energy than the initial states, and the reaction barriers were calculated to be 3, 7, and 8 kcal mol⁻¹ for FePIm, model1, and model2, respectively (Fig. 3a). For those reactions with NO, the energy reductions of the final states were found to be 5–8 kcal mol⁻¹, being much smaller than those for the O₂ and CO cases, and the activation barriers were 2, 3, and 6 kcal mol⁻¹ for FePIm, model1, and model2, respectively (Fig. 3b). The final states of heme–CO⋯H₂ were also found to have lowered in energy by 25–28 kcal mol⁻¹ compared with the initial states of heme–H₂⋯CO. For these displacement reactions, the activation barriers were calculated to be 3, 1, and 11 kcal mol⁻¹ for FePIm, model1, and model2, respectively (Fig. 3c).

Fig. 3 Energy profiles (bottom panels) for displacement reactions of a H₂ molecule (a–c) or H atom (d–f) bound to heme modeled using FePIm, model1 and model2. The optimized geometries (top panels) of molecular complexes composed of FePIm, H₂ or H, and a DM, namely O₂ (a and d), NO (b and e), and CO (c and f) are shown, corresponding to the initial state (IS), transition state (TS) and final state (FS).

From the optimized structures of the molecular complexes corresponding to the IS and FS for the FePIm cases (Fig. 3a–c), it was found that the adsorbed structures of H₂ and DMs onto FePIm were almost preserved, even in the presence of another molecule, for the initial and final states. The shortest distances between the adsorbed H₂ and DM in the initial states were found to be 2.27, 3.34, and 3.22 Å for O₂, NO, and CO, respectively. In the final states, the distances between the adsorbed DM and H₂ were 2.83, 3.99, and 3.52 Å for O₂, NO, and CO. Meanwhile, the transition states were observed at a distance over FePIm, and the distances between the H₂ and DM were 2.03, 2.79, and 2.42 Å for O₂, NO, and CO. Similar findings were obtained for models 1 and 2 (see Fig. S2†).

As found in this work, the activation barriers for displacement reactions of H₂ by O₂, NO, and CO on heme were quite small, ranging from 1 to 11 kcal mol⁻¹, indicating that these reactions can occur spontaneously at finite temperature, as verified through AIMD simulations. Therefore, one can suppose that molecular hydrogen absorbed onto Hb in a blood cell can be readily displaced by other DMs, especially O₂, thereby with no expected harmful influence on the storage and transport of oxygen in the respiratory chain. The same reasoning was also applied to other DMs, such as NO and CO, which exhibit a medicinal effect of propagating intercellular signals. When compared between the three DMs, the energy lowering of the final state was found to be the smallest for the case of NO, while those of O₂ and CO were more or less comparable to each other and much larger than that of NO, but the activation barriers were found to be in a similar order. This indicates that O₂ and CO can be more favourable for replacing H₂ than NO. Interestingly, although they are very simplified structures, models 1 and 2 resembled the displacement reaction paths on FePIm for the three DM cases, in accordance with the binding tendencies.

Then, the displacement reactions of a H atom adsorbed onto heme by other DMs were considered. Fig. 3(d–f) present the energy profiles and optimized geometries of the molecular complexes obtained via NEB runs. In contrast to the H₂ cases, the final states of heme–DM⋯H were found to be energetically higher than the initial states of heme–H⋯DM, except in the cases of FePIm–NO and –CO. Furthermore, the activation barriers were calculated to be much higher compared to the cases of H₂, indicating that the displacement reaction of a H atom adsorbed on heme could require energy from external sources. This might be related to the stronger binding of H to the Fe of heme compared with that of H₂. For the cases of O₂, the final states were 35, 27, and 23 kcal mol⁻¹ higher in energy than the initial states, and the activation barriers were calculated to be 55, 35, and 40 kcal mol⁻¹ for FePIm, model1 and model2, respectively. When displacing H with NO, the energy differences between the final and initial states were −6, 16, and 15 kcal mol⁻¹, and the activation barriers were 5, 46, and 38 kcal mol⁻¹ for FePIm, model1 and model2. For the CO cases, the energy differences of the final states from the initial states were −15, 61, and 58 kcal mol⁻¹, and the activation energies were calculated to be 11, 35, and 45 kcal mol⁻¹ for FePIm, model1 and model2. It is distinctive for models 1 and 2 that the relative energies of the transition states are lower than those of the final states and relatively flat plateaus are observed in the energy profiles.

The optimized geometries corresponding to IS and FS were also shown to be well-reproduced compared with the adsorbed structures of H and DM without another molecule or atom. The distances between O₂ and H were found to be the smallest (2.20, 1.82 and 2.79 Å in the IS, TS and FS; Fig. 3d) among the three DMs, indicating a strong interaction between O₂ and H. Meanwhile, the distances between CO and H were observed to be the largest (3.37, 3.31 and 4.29 Å in the IS, TS and FS; Fig. 3f), proving the weak interaction between CO and H. Similar findings were observed for models 1 and 2 (see Fig. S3†). Such strong and weak interactions between the adsorbates may cause higher and lower activation barriers for the displacement reactions.

In most of the cases, the activation barriers were higher, over 30 kcal mol⁻¹, compared with the H₂ displacement reactions, indicating that the H displacement could not occur readily but H and a DM could coexist on heme for a long time. Such coexistence of H and O₂ or CO was verified through AIMD simulations carried out at room temperature. Therefore, one can surmise that a H atom adsorbed on heme may disturb the regular respiration of cells with hemoglobin or myoglobin. However, the activation energy for H₂ dissociation on protoheme was found to be as high as about 64 kcal mol⁻¹ in our previous work.¹⁵ Exceptionally, the activation energy for the reaction with NO on FePIm was found to be very low (5 kcal mol⁻¹), indicating the spontaneous displacement of H by NO, as confirmed by AIMD simulation. This might be associated with the unique interaction between NO and heme, which have unpaired electrons and thus make lone-pair electrons forming a dative chemical bond.

3.3 Kinetics for reactions based on model2

Finally, we calculated the rate constants (k) of H₂ or H displacement reactions with DMs on model2 by applying the CVT and RP-VTST formalisms, as implemented in Polyrate using the output of a DIRDYVTST run of NWChem. The equilibrium constants (K) were derived from the calculated rate constants. The backward as well as forward reactions were also considered. The temperature T was changed, ranging from 250 to 1010 K with a step of 30 K.

Fig. 4 illustrates the calculated k(T) and K(T) for the H₂ displacement reactions. It was confirmed that the calculated rate constants, both for the forward and backward reactions, satisfied the general trend of an exponential increase with increasing temperature according to the Arrhenius equation k(T) = A exp(−E_a/RT), where R = 8.31 J mol⁻¹ K is the gas constant and E_a is the activation energy. Table 3 lists the determined fitting parameters. For the forward reactions of model2–H₂ + DM → model2–DM + H₂, the activation energies were determined to be 11, 3, and 8 kcal mol⁻¹ for O₂, NO, and CO respectively, which are in reasonable agreement with those (8, 6, and 11 kcal mol⁻¹) from the NEB calculations. The rate constants for the O₂ reaction show the highest values in magnitude while those for CO exhibit the lowest values among the three DMs, indicating that the O₂ displacement reaction has the slowest rate and the NO reaction has the fastest rate (Fig. 4a). The equilibrium constants exhibit a decreasing tendency with increasing temperature. In accordance with the tendency of the rate constants, the equilibrium constants of the O₂ reaction have the lowest values while those of the CO reaction show the highest values in the whole range of temperatures. The rate and equilibrium constants of the NO reaction have moderate values, but the lowest value of the activation energy. For the backward reactions, the activation energies were calculated to be higher than those for the forward reactions, 41, 18, and 38 kcal mol⁻¹ for the O₂, NO, and CO cases, respectively, in agreement with the NEB calculations. In contrast to the forward reactions, the equilibrium constants for the backward reactions show an increasing trend with increasing temperature.

Fig. 4 Rate constants (k; black colour) and equilibrium constants (K; blue colour) of H₂ displacement reactions with O₂, NO and CO on model2 in the (a) forward and (b) backward directions as functions of temperature. Red-coloured solid lines indicate the regressions into the Arrhenius equation.

Table 3 Parameters A (units: cm³ per molecule s) and B (units: K) obtained by fitting the rate constants into the Arrhenius equation k(T) = A exp (−B/T) and the determined activation energy E_a = RB (units: kcal mol⁻¹) for displacement reactions on model2

Displacement reaction	Forward			Backward
	A	B	Ea	A	B	Ea
O2–H2	2.5 × 10^−9	5545	11	5.0	20 774	41
NO–H2	2.7 × 10^−10	1537	3	4.5 × 10^2	9254	18
CO–H2	3.8 × 10^−6	3761	8	3.3 × 10^19	19 234	38
O2–H	0.2	19 753	39	1.2 × 10^−17	9618	19
NO–H	6.0 × 10^−11	19 918	40	1.0 × 10^−8	11 720	23
CO–H	3 × 10^8	21 763	43	0.9	11 488	23

---

The calculated rate and equilibrium constants for the H displacement reactions of model2–H + DM → model2–DM + H are presented in Fig. 5. As listed in Table 3, the activation energies for the forward reactions were determined to be 39, 40, and 43 kcal mol⁻¹ for the O₂, NO, and CO cases, respectively, in agreement with the NEB calculations (40, 38, and 45 kcal mol⁻¹). For the backward reactions, the activation energies were found to be lower at 19, 23, and 23 kcal mol⁻¹ for the O₂, NO, and CO cases, respectively. The H displacement reaction with CO might occur most rapidly among the DMs, whereas the reaction with NO might occur most slowly, due to them having the highest and lowest values of rate constants. For the backward reactions, the CO reaction also occurs most rapidly, while the O₂ reaction progresses most slowly among the DMs. The equilibrium constants for the forward reactions show an increasing trend as temperature increases, but those for the backward reactions show a decreasing tendency.

Fig. 5 Rate constants (k; black colour) and equilibrium constants (K; blue colour) of H displacement reactions with O₂, NO and CO on model2 in the (a) forward and (b) backward directions as functions of temperature. Red-coloured solid lines indicate the regressions into the Arrhenius equation.

4 Conclusions

In this work, we have studied the binding properties of molecular hydrogen to heme in competition with DMs, and their displacement reactions, using DFT and VTST calculations. The heme group was modeled using FePIm, model1 and model2, while various molecular complexes were constructed by combining DMs, namely H₂, O₂, NO and CO, or a H atom with the modeled heme groups. From the analysis of the structures and energetics of the molecular complexes, heme–DM and –H, which were optimized using the B3LYP/6-311G and LANL2DZ_ECP method, it was found that the Fe–DM bond angles were similar among the three heme models and the O₂ complexes have the lowest hardness and the highest electrophilicity for the three models. Moreover, the general tendencies of the structural properties and reactivity in models 1 and 2 resemble those of FePIm, indicating that the simplified models can be used for further calculations. The difference was in the order of the binding strength of DMs or H to the Fe of heme, being NO > O₂ > CO > H > H₂ for the FePIm-based complexes but H > O₂ > NO > CO > H₂ for the model-based complexes. Through AIMD simulations of FePIm–H₂⋯DM and FePIm–H2⋯DM complexes, it was observed that on FePIm, H₂ could be readily replaced by O₂ and NO but not by CO, whereas H could be displaced by NO but not by O₂ and CO at room temperature. The activation energies for H₂ or H displacement reactions with other DMs were determined by applying the NEB method, revealing that the H₂ displacement by O₂, NO and CO could occur almost spontaneously due to the very low activation energies of below 11 kcal mol⁻¹, but H-atom displacement should require a large amount of energy, due to the high activation energies of over 30 kcal mol⁻¹ for the three heme models. Finally, we calculated the rate and equilibrium constants as functions of temperature for the displacement reactions based on model2, confirming that the rate constants satisfy the general tendency of increasing with increasing temperature and the activation energies determined by fitting into the Arrhenius equation agreed well with those from the NEB calculations. With these findings, we believe that the present work can contribute to the understanding of the efficacy of molecular hydrogen towards vascular diseases and the development of enhanced medicine.

Author contributions

Chol-Jun Yu developed the original project. Yun-Kyong Ri, Song-Ae Kim and Chol-Jun Yu performed the calculations and drafted the first manuscript. Yun-Hyok Kye, Yu-Chol Jong and Myong-Su Kang assisted with the post-processing of calculation results, and contributed to useful discussions. Chol-Jun Yu supervised the work. All authors reviewed the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is supported as part of the basic research project “Design of New Energy Materials” (No. 2021-12) funded by the State Commission of Science and Technology, DPR Korea. Computations have been performed on the HP Blade System C7000 managed by the Faculty of Materials Science, Kim Il Sung University.

References

1. T. W. LeBaron, B. Kura, B. Kalocayova, N. Tribulová and J. Slezák, Molecules, 2019, 24, 2076.
2. D. Wang, L. Wang, Y. Zhang, Y. Zhao and G. Chen, Biomed. Pharmacother., 2018, 104, 788–797.
3. T. Hasegawa, M. Ito, S. Hasegawa, M. Teranishi, K. Takeda, S. Negishi, H. Nishiwaki, J. Takeda, T. W. LeBaron and K. Ohno, Int. J. Mol. Sci., 2022, 23, 2888.
4. R. Yamamoto, K. Homma, S. Suzuki, M. Sano and J. Sasaki, Sci. Rep., 2019, 9, 1255.
5. S. Ohta, Biochim. Biophys. Acta, Gen. Subj., 2012, 1820, 586–594.
6. M. Dole, F. R. Wilson and W. P. Fife, Science, 1975, 190, 152–154.
7. I. Ohsawa, M. Ishikawa, K. Takahashi, M. Watanabe, K. Nishimaki, K. Yamagata, K. ichiro Katsura, Y. Katayama, S. Asoh and S. Ohta, Nat. Med., 2007, 13, 688–694.
8. Q. Hu, Y. Zhou, S. Wu, W. Wu, Y. Deng and A. Shao, Biomed. Pharmacother., 2020, 130, 110589.
9. B. Kura, A. K. Bagchi, P. K. Singal, M. Barancik, T. W. LeBaron, K. Valachova, L. Šoltés and J. Slezák, Can. J. Physiol. Pharmacol., 2019, 97, 287–292.
10. Y. Ming, Q. Ma, X. Han and H. Li, Exp. Ther. Med., 2020, 20, 359–366.
11. M. Barancik, B. Kura, T. W. LeBaron, R. Bolli, J. Buday and J. Slezak, Antioxidants, 2020, 9, 1281.
12. Y. Zhang, S. Tan, J. Xu and T. Wang, Cell. Physiol. Biochem., 2018, 47, 1–10.
13. T. Kiyoi, S. Liu, E. Takemasa, H. Nakaoka, N. Hato and M. Mogi, PLoS One, 2020, 15, e0227582.
14. A. Tamasawa, K. Mochizuki, N. Hariya, M. Saito, H. Ishida, S. Doguchi, S. Yanagiya and T. Osonoi, Eur. J. Pharmacol., 2015, 762, 96–101.
15. S.-A. Kim, Y.-C. Jong, M.-S. Kang and C.-J. Yu, J. Mol. Model., 2022, 28, 287.
16. D. Moris, M. Spartalis, E. Spartalis, G.-S. Karachaliou, G. I. Karaolanis, G. Tsourouflis, D. I. Tsilimigras, E. Tzatzaki and S. Theocharis, Ann. Transl. Med., 2017, 5, 326.
17. E. Hood, E. Simone, T. Wattamwar, P. Dziubla and V. Muzykantov, Nanomedicine, 2011, 6, 1257–1272.
18. H. Li, Y. Luo, P. Yang and J. Liu, J. Neurol. Sci., 2019, 396, 240–246.
19. S. Takeuchi, K. Nagatani, N. Otani, H. Nawashiro, T. Sugawara, K. Wada and K. Mori, BMC Neurosci., 2015, 16, 22.
20. C. S. Huang, T. Kawamura, Y. Toyoda and A. Nakao, Free Radical Res., 2010, 44, 971–982.
21. Y.-S. Ge, Q.-Z. Zhang, H. Li, G. Bai, Z.-H. Jiao and H.-B. Wang, Hepatobiliary Pancreatic Dis. Int., 2019, 18, 48–61.
22. J. Yu, Q. Yu, Y. Liu, R. Zhang and L. Xue, PLoS One, 2017, 12, e0173645.
23. C. H. Chen, A. Manaenko, Y. Zhan, W. W. Liu, R. P. Ostrowki, J. Tang and J. H. Zhang, Neuroscience, 2010, 169, 402–414.
24. G. J. Kubas, R. R. Ryan, B. I. Swanson, P. J. Vergamini and H. J. Wasserman, J. Am. Chem. Soc., 1984, 106, 451–452.
25. Y. Jean, O. Eisenstein, F. Volatron, B. Maouche and F. Sefta, J. Am. Chem. Soc., 1986, 108, 6587–6592.
26. P. J. Hay, J. Am. Chem. Soc., 1987, 109, 705–710.
27. P. J. Hay, Chem. Phys. Lett., 1984, 103, 466.
28. G. J. Kubas, Acc. Chem. Res., 1988, 21, 120–128.
29. J. Oláh and J. N. Harvey, J. Phys. Chem. A, 2009, 113, 7338–7345.
30. I. Degtyarenko, R. M. Nieminen and C. Rovira, Biophys. J., 2006, 91, 2024–2034.
31. L. M. Blomberg, M. R. A. Blomberg and P. E. M. Siegbahn, J. Inorg. Biochem., 2005, 99, 949–958.
32. M. Tsuda, E. S. Dy and H. Kasai, Eur. Phys. J. D, 2006, 38, 139–141.
33. J. P. Arcon, P. Rosi, A. A. Petruk, M. A. Marti and D. A. Estrin, J. Phys. Chem. B, 2015, 119, 1802–1813.
34. H. Chen, W. Lai and S. Shaik, J. Phys. Chem. B, 2011, 115, 1727–1742.
35. N. Strickland and J. N. Harvey, J. Phys. Chem. B, 2007, 111, 841–852.
36. I. Degtyarenko, X. Biarnés, R. M. Nieminen and C. Rovira, Coord. Chem. Rev., 2008, 252, 1497–1513.
37. M. Tsuda, W. A. Dino, H. Nakanishi and H. Kasai, Chem. Phys. Lett., 2005, 402, 71–74.
38. M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus and W. A. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489.
39. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
40. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
41. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 299–310.
42. F. B. van Duijneveldt, J. G. C. M. Van Duijneveldt-Van De Rijdt and J. H. Van Lenthe, Chem. Rev., 1994, 94, 1873–1885.
43. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
44. M. Radón, E. Broclawik and K. Pierloot, J. Phys. Chem. B, 2010, 114, 1518–1528.
45. B. Long, Y. Wang, Y. Xia, X. He, J. L. Bao and D. G. Truhlar, J. Am. Chem. Soc., 2021, 143, 8402–8413.
46. J. L. Bao and D. G. Truhlar, Chem. Soc. Rev., 2017, 46, 7548–7596.
47. R. Steckler, W.-P. Hu, Y.-P. Liu, G. C. Lynch, B. C. Garrett, A. D. Isaacson, V. S. Melissas, D. h. Lu, T. N. Truong, S. N. Rai, G. C. Hancock, J. G. Lauderdale, T. Joseph and D. G. Truhlar, Comput. Phys. Commun., 1995, 88, 341–343.
48. J. Vojtechovsky, K. Chu, J. Berendzen, R. M. Sweet and I. Schlichting, Biophys. J., 1999, 77, 2153–2174.
49. O. N. Chen, S. Groh, A. Liechty and D. P. Ridge, J. Am. Chem. Soc., 1999, 121, 11910–11911.
50. R. G. Parr, L. V. Szentpaly and S. Liu, J. Am. Chem. Soc., 1999, 121, 1922–1924.

---

Footnote

† Electronic supplementary information (ESI) available: Tables for binding energetics and Mulliken charges, and figures for fluctuations in temperature and energy during AIMD simulations and optimized geometries corresponding to the IS, TS and FS during NEB runs for models. See DOI: https://doi.org/10.1039/d4ra02091j

---