Antioxidation activity of molecular hydrogen via protoheme catalysis in vivo: an insight from ab initio calculations

Recently, molecular hydrogen has been found to exhibit antioxidation activity through many clinical experiments, but the mechanism has not been fully understandable at atomic level. In this work, we perform systematic ab initio calculations of protoheme-hydrogen complexes to clarify the antioxidation mechanism of molecular hydrogen. We make molecular modeling of iron–protoporphyrin coordinated by imidazole, FeP(Im), and its hydrogen as well as dihydrogen complexes, together with reactive oxygen/nitrogen species (RONS). We carry out structural optimization and Mulliken charge analysis, revealing the two kinds of bonding characteristics between FeP(Im) and H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}: dihydrogen bonding in the end-on asymmetric configuration and Kubas bonding in the side-on symmetric configuration of H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document} molecule. The activation barriers for adsorption and dissociation of H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document} on and further desorption of H atom from FeP(Im) are found to be below 2.78 eV at most, which is remarkably lower than the H–H bond breaking energy of 4.64 eV in free H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document} molecule. We find that the hydrogen bond dissociation energies of FeP(Im)–H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document} and –H complexes are lower than those of RONS–H complexes, indicating the decisive role of protoheme as an effective catalyst in RONS antioxidation by molecular hydrogen in vivo.


Introduction
Molecular hydrogen (or hydrogen gas, H 2 ) has been attracting significant attention as an effective antioxidant against oxidative stress since the first report of its medical efficacy [1]. In fact, H 2 selectively reduces highly cytotoxic reactive oxygen/nitrogen species (ROS/RNS or RONS), being the main cause of oxidative stress in cells, such as hydroxyl radical ( • OH) and peroxynitrite (ONOO − ) [2,3]. Moreover, it has been demonstrated from many cellular and animal experiments that molecular hydrogen exerts clear preventive and therapeutic effects on a variety of pathologies [4][5][6], including metabolic syndrome [7][8][9], diabetes [10], cancer [11], chronic hepatitis B [12], ischemia/ reperfusion (I/R) injury [13,14], radiation injury [15,16], etc. In addition to the antioxidative action, hydrogen gas has more beneficial functions such as anti-inflammation, antiapoptosis, anti-allergy and stimulation of energy metabolism [17][18][19][20]. Therefore, together with its non-toxicity and biological safety, molecular hydrogen can be said to have a bright prospect in future research and clinical applications. In spite of such successes and prospect, mechanisms of H 2 action are in debate.
Among the existing molecules, molecular hydrogen is the smallest molecule in size and mass, being even smaller than molecular oxygen [2,7]. In addition, it is colorless, odorless, non-polar and electrically neutral. These unique physical properties allow for its easy penetration across cell membranes and rapid diffusion into mitochondria, endoplasmic reticulum and nuclei [3], which are the primary sites of RONS generation. The controlled generation of RONS extracellularly by some enzymes such as NADPH-oxidase and nitric oxide synthase was developed evolutionarily to form the innate immune system for killing bacteria, whereas their excessive release or intracelluarly release as a coproduct of metabolic processes was thought to cause oxidative damage to host biomolecules, being referred to oxidative stress [21]. The oxidative stress is believed to largely contribute to the progress of a variety of diseases, especially age-related diseases. When being delivered via various methods, such as inhaling H 2 gas, drinking H 2 -rich water and injecting H 2 -rich saline, dihydrogen was reported to scavenge only the cytotoxic and strong excessive RONS selectively, but not react with advantageous RONS involved in intracellular signaling, cellular metabolic redox reactions, or physiological enzymatic and metabolic reactions [22][23][24][25]. It was also reported that there were no toxic and adverse effects of H 2 even in higher concentrations above the clinical treatment dosages [26,27].
Nowadays, there are several molecular mechanisms of H 2 action, which can be considered mainly in two ways. The first way is direct mechanisms, where H 2 reacts with hydroxyl radicals ( • OH) and peroxynitrite (ONOO − ), resulting in reduction of RONS [1,28,29]. In the second indirect way, the oxidative stress can be reduced by inducing antioxidation systems, such as myeloperoxidase [30], superoxide dismutase [10], catalase [30], and heme oxygenase-1 (HO-1) [6]. These enzymes were found to mediate the H 2 effects through modulation of signal transduction [31,32] and alteration of gene expressions [33,34]. In this line, it has been reported that matrix metalloproteinases (MMPs) can be effective mediators of H 2 effects [28,35]. In fact, inhalation of molecular hydrogen was reported to significantly reduce the expression of MMP-9, which promotes hemorrhagic infarction by bursting cerebral vessels in a rat model [36]. Shi et al. [37] proposed a hypothesis that H 2 could bind to metal in specific metalloproteins as a ligand and regulate the metalloproteins by M-H 2 bonding. There exist several previous works for dihydrogen-metalloprotein complexes [38,39]. However, details of formation and molecular structure of dihydrogen-metalloprotein complexes have yet to be fully understandable in relating with antioxidative effects of H 2 .
In this work, we perform density functional theory (DFT) calculations of protoheme-hydrogen complexes aiming at elucidating the mechanism of antioxidation effects of molecular hydrogen in vivo. We optimize molecular structures of the complexes and evaluate their formation and activation energies for adsorption, dissociation, and desorption of H 2 molecule on protoheme. The antioxidative efficacy is estimated by calculating the bond dissociation energy (BDE) of the complexes in comparison with other antioxidants. We analyze the electronic charge transfer during reaction, frontier molecular orbitals (FMOs) including highest occupied and lowest unoccupied ones (HOMO and LUMO), and charge population.

Molecular models
In this work, we consider hemoglobin (Hb) as matrix metalloprotein for binding H 2 molecule, since in vivo transition metal atoms such as Fe, Co, and Cu are present in heme proteins or non-heme proteins (e.g., ferredoxin). In the structure of Hb, heme is the active site of Hb and is constructed from iron-protoporphyrin (FeP) complex, where Fe can be readily bound with small molecules like O 2 and thus play a key role in transporting oxygen in the vascular systems of animals. In addition, Hb is coordinated by imidazole ligand, which is representative of a histidine amino acid binding to heme. This is referred to protoheme with a formula FeP(Im), as shown in Fig. 1. From the accumulated studies of Fe-O 2 binding, we can expect that H 2 also binds to Fe in FeP(Im), resulting in formation of protoheme-dihydrogen complex, FeP(Im)-H 2 .
For comparison of the antioxidation activity of protoheme-dihydrogen complex, we chose some other molecules that have a structural feature of antioxidant, including phenol ( Fig. S1 for their molecular structures). All of these molecules readily release hydrogen radical for antioxidation effect. Among these, vitamin C is namely ascorbic acid which acts as a hydrogen donor, and vitamin E is namely -tocopherol, a fat-soluble and free radical chain breaking antioxidant, which is also an effective hydrogen donor due to the presence of hydroxyl group in its structure.

Computational details
The DFT calculations were carried out using the pseudopotential -pseudo atomic orbital (PAO) method as implemented in the SIESTA code (version 4.1.b3) [40,41]. For the electrostatic interaction between the ionic core and valence electrons, we constructed the soft normconserving pseudopotentials proposed by Troullier and Martins [42] using the ATOM program included in the SIESTA package, where the valence electron configurations of atoms are C-2s 2 2p 2 3d 0 4f 0 , N-2s 2 2p 3 3d 0 4f 0 , O-2s 2 2p 4 3d 0 4f 0 , H-1s 1 2p 0 3d 0 4f 0 , and Fe-4s 2 4p 0 3d 6 4f 0 . As can be hinted from the configurations, we included the empty states up to l max = 3 (i.e., f state) in constructing the pseudopotentials by applying the generalized approach [43]. The cutoff radii were 1.49, 1.54, 1.14, 0.15, and 2.0 bohr for the s-f components equally in C, N, O, H, and Fe, respectively. To describe the exchangecorrelation interaction among the valence electrons, we used the Perdew-Burke-Ernzerhof (PBE) functional [44] within the generalized gradient approximation (GGA). The dispersive van der Waals (vdW) interactions between molecules were considered by employing the semiempirical Grimme's approach [45] with the parameters provided in the package.
For PAOs, we used a split-valence double-and polarization (DZP) basis set for all atoms, with a split norm of 0.25 and an energy shift of 50 meV. The kinetic energy cutoff for setting the wavelength of the shortest plane wave was fixed to be 150 Ry, which provides a real spacing between the grid points of 0.13 Å for projecting the electron density. We checked that these computational parameters of the real space grid and PAO basis sets yield well converged results, as already proved in the previous studies of heme complexes [46][47][48]. In the atomic relaxations, the atoms were allowed to relax until the forces converged to 0.02 eV/Å. To calculate the energy barriers for reactions, we employed the nudged elastic band (NEB) method [49] as implemented in Python script Pastafarian 1 in conjunction with SIESTA code as applied in our previous work [50,51]. During the NEB run, all the atoms were allowed to relax with the convergence threshold for force of 0.05 eV/Å. The number of images was adopted to be 15 in this work.
We used a cubic supercell with a lattice constant of 60 Å, in which a molecule is placed with a sufficiently long distance from the neighboring periodic image. The basis set superposition error was corrected using the counterpoise method [52,53] The molecular structures and volumetric data were visualized using the VESTA code [54].

Evaluation of antioxidation activity
The direct scavenging processes of RONS free radicals in vivo by antioxidants are considered mainly according to the two mechanisms [55]. The first mechanism is the direct hydrogen radical transfer mode from antioxidant to RONS as follows, where RONS • is the reactive oxygen/nitrogen species radical, R-H is the antioxidant molecule, RONS-H is the reduction product, and R • is the antioxidant radical. Here, the antioxidation activity can be elucidated by how easily the antioxidant (R-H) can release the free hydrogen radical. This can be evaluated by calculating the hydrogen bond dissociation energy (H-BDE) of antioxidant. Supposing that the dissociation of antioxidant is R-H → R • + H • , its H-BDE can be calculated using the total energies as follows, It is clear that the lower the H-BDE is, the stronger the antioxidation activity is.
The second mechanism is the electron-proton transfer mode between antioxidant and RONS like, In the first step of this reaction, the antioxidant cation (R-H + ) is formed. Therefore, this process is mainly governed by the ionization potential (IP) of the antioxidant, which can be obtained directly from the eigen value of the highest occupied molecular orbital ( HOMO ), while the electron affinity (EA) is estimated from the lowest unoccupied molecular orbital ( LUMO ) [56], The lower IP is favorable for stronger antioxidation activity.

Structures of protoheme and protoheme-H 2 complexes
We started with optimization of protoheme, FeP(Im), and protoheme-H 2 complexes, FeP(Im)-H 2 . Here H 2 was suggested to bind to Fe placed at the center of protoporphyrin in the opposite direction to imidazole. In binding H 2 to protoheme, we made two different models of asymmetric end-on and symmetric side-on Fe-H 2 configurations in reference to the previous theoretical and experimental works for Fe-O 2 binding in Hb [57,58]. Figure 2 shows their optimized structures with the isosurface plot of charge density difference upon H 2 binding to protoheme only in the cases of protoheme-H 2 complexes.
The main structural parameters defining the optimized FeP(Im) and FeP(Im)-H 2 models are listed in Table 1 (see atom labeling in Fig. 2). In FeP(Im), the bond length of Fe-N por , where N por are N atoms of the porphyrin ring, was found to range from 1.992 to 2.039 Å with an average length of 2.020 Å, and the bond length between Fe and N ax of ligand imidazole was optimized to be 1.895 Å. These values agree well with the experimental and previous theoretical ones. When binding H 2 to Fe of heme, the Fe-N por bond length was found to be slightly increased to 2.024 and 2.025 Å for asymmetric and symmetric H 2 configurations, while the Fe-N ax was more or less unchanged for the former case but clearly increased to 1.989 Å for the latter case. This indicates that the symmetric H 2 binding mode affects more apparently the host structure than the asymmetric mode.
For the characteristic bonding structures of protoheme-H 2 complexes, we measured the Fe-H bond length and the relevant bond angles. The Fe-H distances in the asymmetric mode are 1.768 and 2.267 Å, which are clearly larger than those of 1.698 and 1.718 Å in the symmetric mode, indicating stronger binding of H 2 to Fe in the latter case. When compared with FeP(Im)-O 2 , these Fe-H distances are somewhat similar to Fe-O distance (1.78 Å) but smaller than Co-O distance (1.94 Å) obtained with the same method to this work [48]. Note that the experimental values are 1.81 and 2.03 Å for Fe-O and Co-O, respectively [59,60]. When compared with 0.726 Å in free H 2 molecule, obtained in this work (the experimental value is 0.74 Å), the H-H bond length in FeP(Im)-H 2 complex is decreased to 0.703 Å for the asymmetric mode while increased to 0.731 Å for the symmetric mode. The Fe-H-H angle in the asymmetric case   [47], and those are 76.1 and 79.2 • in the symmetric case. We note that in the side-on symmetric conformation the Fe-H distance and Fe-H-H angle are slightly deviated from symmetry, which is due to the asymmetry of porphyrin ring with different side groups. The N 3 -Fe-H-H dihedral angles, which define the orientation of the H 2 molecule to the porphyrin plane, are also measured as 42.0 and 35.1 • in the asymmetric and symmetric modes, respectively. For the asymmetric mode, the value is close to 45 • , indicating that the H 2 molecular axis almost bisects a porphyrin quadrant, being similar to that for the case of O 2 binding [48].
In binding H 2 to the protoheme, electron transfer occurs. As shown in Fig. 2(b) and (c) for electronic charge density difference upon binding, the electron transfer is bilateral between protoheme and H 2 . That is, the electron gain (darkred color) and loss (green) are found around both Fe and H 2 sides. To understand this, we remind the general bonding feature between Fe and H 2 in dihydrogen complexes [61]. The H 2 molecule has two molecular orbitals (MOs): twoelectron bonding and empty anti-bonding * . When contacting with Fe atom that has filled and empty d orbitals, the H 2 molecule act as a neutral two-electron donor on one hand, while the retrodonation also occurs from the filled d orbitals of Fe to the anti-bonding * MO of the hydrogen molecule on the other hand, forming the so-called Kubas bonding [61]. Such Kubas bonding was well observed in the symmetric configuration. Meanwhile, for the case of asymmetric end-on binding, the hydrogen atom closer to the host acts as an acceptor by seeing electron accumulation around it, whereas the upper hydrogen atom acts as a donor, forming the dihydrogen bond (DHB) [62]. For quantitative analysis of electron transfer, we show the Mulliken charge of species in Table 1 (see Table S1 for details). The Fe charges were determined to be positive in protoheme and protoheme-H 2 , indicating electron loss, while the porphyrin has negative charge, implying electron gain. Interestingly, for the binding H 2 molecule, one H atom has positive charge of 0.014 (electron donor) and another has negative charge of −0.037 (acceptor) for the case of asymmetric orientation. On contrary, both H atoms in the symmetric mode have similar negative charges, indicating that H 2 molecule acts as acceptor. Moreover, we found that the amount of electron transfer in the symmetric mode is larger than that in the asymmetric one, indicating that the former case has weaker binding of H 2 than the latter case in accordance with the bond length tendency. This is somewhat contradictory to the previous insight that H 2 molecule is easily adsorbed on the top of metal atom of alloy cluster with end-on orientation rather than side-on orientation [63].

Dissociation energetics of H 2 on protoheme
Hydrogen molecule is chemically useful only when the strong two-electron H-H bond is catalytically broken in a controlled way. Transition metals such as Fe and Co supported on Hb or Mb are very effective and widely used catalysts for the splitting of H 2 . Therefore, it is of great importance to clarify energetics of adsorption, dissociation, and desorption of hydrogen molecule catalyzed by Fe atom. To this end, we applied the NEB method to find out the transition states with the activation barrier for these three processes. Figure 3 shows the calculated energy profiles for adsorption, dissociation, and desorption of H 2 on FeP(Im) via the asymmetric end-on and symmetric side-on orientation of H 2 molecule. We set the total energy of the optimized FeP(Im)-H 2 adsorbed complex with the end-on orientation to the reference 0 eV. It was found that the FeP(Im)-H 2 with the side-on mode is energetically lower by 0.773 eV than that with the end-on orientation. This is contrary to the case of O 2 adsorption on heme [57,58], but agrees with the insight into bonding between Fe and H 2 as discussed above. For adsorption of H 2 onto protoheme, we performed the optimization of initial state, where H 2 molecule was far away from the protoheme with a distance of 7.76 Å, and found its total energy higher by 0.076 eV above the reference state. The NEB runs for the adsorptions yielded the activation barriers of 0.129 eV and 0.138 eV with the total energies of 0.205 and 0.214 eV above the reference state for the end-on and side-on configurations, respectively (see Fig. S3 for geometries and energy profiles). The geometries of the transition states were found to be almost identical with the Fe-H distance of ∼5.6 Å in the end-on orientation. Despite the negligible difference in the calculated activation barriers, the H 2 adsorption on protoheme with the end-on orientation can be said to occur more easily than that with the side-on mode.
Then, let us consider the dissociation of adsorbed H 2 on protoheme. In breaking the H-H bond, one H atom is attracted by the nearest C atom because bonding with the orbital of the C-C bond is energetically favorable. The optimization of this H 2 -dissociated complex, denoted as FeP(Im)-H-H, gave the H-H distance to be 3.132 Å and Fe-H bond length to be 1.545 Å, as shown in Fig. 4a. The energy of FeP(Im)-H-H complex was determined to be higher by 0.780 eV above the reference state, giving the H-BDE to be 0.780 eV and 1.553 eV for the end-on and side-on modes, respectively. Therefore, it can be said that the end-on mode is more favorable for H 2 dissociation than the side-on mode. Nevertheless, these are very lower than the dissociation energy of free molecular hydrogen (4.644 eV in this work and 4.510 eV in experiment), confirming the significant advantage of using FeP(Im) as a catalyst. The activation barriers for the H-H bond breaking were determined to be 2.345 eV and 2.784 eV for the end-on and side-on modes, respectively (see Fig. S4 for geometries and energy profiles). At the transition states, one H atom was found at the middle position between another H atom and C atom of porphyrin ring. In accordance with the aforementioned discussion, the end-on mode has lower activation barrier than the side-on mode, because the H-H and Fe-H bonds (Kubas bond) in the side-on mode are stronger than those (dihydrogen bond) in the end-on mode. It is worth noting that our activation barriers are relatively higher than those for the O-O dissociation on FeP, where, more importantly, the side-on mode (1.352 eV) was found to be easier than the end-on mode (1.975 eV) [58].
At the final step, we considered the desorption process of one H atom bound to C atom. For the final state, the dissociated H atom was suggested to be away by about 4 Å from the C atom of the porphyrin ring, and was optimized. The reaction energy defined by the energy difference between the final and initial states was calculated to be 2.918 eV. The activation energy was found to be 2.138 eV at the transition state with the C-H distance of 2.45 Å(see Fig. S5). It should be noted that the total energy keeps almost constant beyond this distance, indicating that the desorption energy for fully removing one H atom (or moving H atom from heme to infinity) is about 2.138 eV. Then, we performed optimization of FeP(Im)-H complex, which will be used as antioxidant for reducing RONS or free radical as discussed below. Figure 4b shows the optimized structure of FeP(Im)-H complex. The Fe-H bond length was slightly decreased from 1.545 Å in FeP(Im)-H-H to 1.537 Å in FeP(Im)-H due to the enhanced interaction between Fe and H atoms. In FeP(Im)-H complex, the H atom originated from the molecular hydrogen is placed on the central position of the porphyrin ring since the bond angle ∠H-Fe-N 3 approaches to almost 90 • . In this state, the Mulliken charges of porphyrin, Fe, and H were found to be −0.649 , 0.860, and −0.211 , indicating that some amount of electron transfers from Fe to H and porphyrin ring (see Table S1). It should be noted that the free hydrogen atom originated from the molecular hydrogen remains in vivo and thus plays a certain role in reducing the free radicals in vivo.

Antioxidation activity of protoheme-H 2 and -H complexes
The protoheme-dihydrogen or -hydrogen complexes, that is, FeP(Im)-H2 and Fep(Im)-H, can act as antioxidant because they can easily release free hydrogen radical H • , which can reduce the RONS free radicals in vivo. In this work, we chose the typical RONS free radicals in vivo, including • OH, • OOH, ONOO − , C 3 H 7 OO • (aliphatic peroxyl radical), and DPPH • (aromatic peroxyl radical). When these free radicals interact with H radical, they should be reduced to water (H 2 O) for hydroxyl radical, peroxides for peroxyl radicals, and nitrous acid (HOONO) for peroxynitrite. To evaluate the antioxidation activity of FeP(Im)-H 2 and FeP(Im)-H, we calculated their H-BDE from the total energies of the optimized molecules by using Eq. (2) and IP from the HOMO energy by using Eq. (4). For comparison, we also considered the conventional antioxidants such as phenol, quercetin, vitamin C, and vitamin E. Table 2 lists the calculated H-BDE and IP of antioxidants under study and furthermore RONS-H adducts (see Fig. S6 for MOs).
In accordance with the aforementioned discussion, the asymmetric end-on configuration in FeP(Im)-H 2 complex has lower H-BDE (2.994 eV) and IP (3.976 eV) than the symmetric side-on configuration (3.767 eV and 4.044 eV), indicating that the former is more favorable for antioxidation than the latter. The protoheme-hydrogen radical complex, FeP(Im)-H, can release H • radical more easily than the protoheme-dihydrogen complexes by seeing the lower H-BDE (2.111 eV), thereby supporting the first mechanism of antioxidation, Eq. (1), but is not more favorable for the second mechanism of electron-proton transfer due to the higher IP value (4.315 eV). When compared with other conventional antioxidants, the end-on FeP(Im)-H 2 and FeP(Im)-H complexes have lower H-BDE values than all others and lower IP values than phenol and vitamin C. Among all the antioxidants considered in this work the FeP(Im)-H is the most favorable for the direct H • transfer mode while vitamin E is the most favorable for the electron-proton transfer mode. It turned out that the protoheme-hydrogen complexes have comparable with or even higher antioxidation activity than the conventional antioxidants. In addition, the energy gaps between HOMO and LUMO vary from 1.13 to 1.36 eV for the protoheme complexes, but from 2.4 to 4.4 eV for the other antioxidants.
In Table 2, we also listed the H-BDE and IP values of RONS-H products. If these values of RONS-H products are higher than those of antioxidants, the reduction of RONS free radicals using the H radical released by antioxidant occurs spontaneously because of stronger interaction between the H • and RONS • radicals. Otherwise, the free H radical should recombine with the antioxidant radical. It should be noted that H 2 has lower H-BDE only than H 2 O but higher than other RONS-H products, meaning that it scavenges only hydroxyl radical. Among the RONS-H products, the lowest H-BDE was found to be 4.152 eV in the ONOO-H product, while the lowest IP to be 5.428 eV in DPPH-H product. These values are clearly higher than those of protoheme-hydrogen complexes as well as other antioxidants except phenol for H-BDE (4.256 eV). This also indicates that using the FeP(Im)-H 2 and -H complexes the antioxidation of RONS is feasible in vivo. We emphasize that the FeP(Im)-H has the lowest H-BDE value, thereby possessing the strongest antioxidation activity among the antioxidants under study in this work.

Reaction between FeP(Im)-H and RONS radicals
In addition to the hydrogen free radical, we suggested that the FeP(Im)-H • radical can directly reduce the RONS free radicals in vivo, resulting in formation of FeP(Im)-H-RONS complexes. To check this suggestion, we further investigated the structural and energetic properties of these reaction products. In Table 3, we summarize the calculation results.
It was found that the Fe-N por bond length is barely affected by combining with RONS, but Fe-N ax and Fe-H bond lengths are clearly increased compared with those in FeP(Im)-H. The bond angle of ∠H-Fe-N 3 is inclined from 89 • in FeP(Im)-H to ∼ 84 • in the reaction complexes. These structural changes are due to the interaction between RONS and protoheme-H complex. We note that among the RONS radicals the hydroxyl radical has the most significant effect on the structural properties of FeP(Im)-H when combining with that. When increasing the molecular weight of RONS, the structural change effect was found to be weakened in general, indicating that the interaction between the RONS and protoheme-H radicals becomes weakening for heavier RONS.
The interaction results in electron transfer between the reactants. To intuitively show the electron transfer, we depicted the electronic charge density difference upon binding of RONS to FeP(Im)-H in Fig. 5. Table 3 presents the Mulliken charges of H, Fe, porphyrin, and RONS, which help give quantitative insight into electron transfer. It is clear that the hydrogen radical bound to and porphyrin of FeP(Im) receive some amount of electron while RONS radicals and Fe lose electron in these five combined complexes.
In order to investigate the energetic properties of these complexes, we calculated the eigen energies of FMOs including  HOMO-1, HOMO, LUMO, and LUMO+1, from which the energy gap ( E gap ) was also obtained, and BDE for dissociating the RONS radical from the FeP(Im)-H radical. Figure 6 shows the isosurface plot of FMOs with their eigen energy values. It was found that HOMO-1 is mainly contributed from Fe 3 d xz , HOMO from 3 d x 2 −y 2 , LUMO from 3 d z 2 , and LUMO+1 from 3 d zy orbitals, together with s and p orbitals of other elements. As listed in Table 3, the highest BDE was found to be 4.434 eV for FeP(Im)-H-OH complex among the five complexes. This is even higher than H-BDE of FeP(Im)-H 2 in the asymmetric configuration (3.994 eV), indicating that FeP(Im)-H can also scavenge the hydroxyl radical. For other RONS radicals, we found the lower BDE than H-BDEs of the antioxidants under study in this work, implying that these RONS free radicals prefer binding with H radical to FeP(Im)-H radical.

Discussion
According to the analysis of H-BDE calculation, the molecular hydrogen H 2 can scavenge only hydroxyl radical but hardly interact with the other RONS radicals in vivo. The dissociation energy of free H 2 molecule was determined to 4.644 eV in agreement with the experimental value of 4.510 eV. According to our calculation, this can be fully compensated by the energy of 5.269 eV released by reaction between H • and • OH radicals in agreement with the experimental value of 5.160 eV. It should be noted that the scavenging reaction of hydroxyl radical with H 2 can be done by the strong reactivity of • OH, not by active role of H 2 . The reactivity of hydroxyl radical is so strong that it can be easily scavenged by any antioxidant. This is in agreement with the previous theoretical studies [64,65]. However, the other radicals are difficult to be scavenged by only H 2 itself due to their relatively lower reactivities.
The clinical experiments confirm that a mild H 2 exhibits various biological effects in vivo [5,6,19,20], implying that either H 2 may be dissociated into H atoms by enzymes with a strong activity or the electronic structure of H 2 may be changed by external effects to break the bond between hydrogen atoms. Such effects can be provided by protoheme including transition metal Fe. The protoheme FeP(Im) catches H 2 molecule with very low activation Fig. 5 Isosurface plot of electronic charge density difference upon binding of RONS radicals to FeP(Im)-H radical at the value of 0.002 |e|/Å 3 . Dark-red (green) color represents electron gain (loss) barriers ( ∼0.139 eV), forming FeP(Im)-H 2 complexes with either Kubas bonding for the symmetric orientation or dihydrogen bonding for the asymmetric orientation. The interaction between FeP(Im) and H 2 bound to it causes easier breaking of H-H bond with the relatively low activation barriers ( ∼2.775 eV) than free H 2 molecule, leading to generation of FeP(Im)-H • and free H • radicals. Not only FeP(Im)-H 2 complexes but also FeP(Im)-H radical can easily release the H • radical due to their lower H-BDE than those of RONS-H complexes, thereby acting as effective antioxidants. Moreover, the FeP(Im)-H • radical itself can react with RONS radicals, which can facilitate the antioxidation effect.

Conclusion
In this work, we have investigated the antioxidation activity of protoheme-hydrogen complexes by using ab initio calculations. We made molecular modelings of FeP(Im)-H 2 Fig. 6 Isosurface view of frontier molecular orbitals including HOMO-1, HOMO, LUMO, and LUMO+1 in FeP(Im)-H-RONS complexes. The corresponding eigen energies are shown in the unit of eV complexes with the end-on asymmetric and side-on symmetric orientations of H 2 and FeP(Im)-H radical as the main antioxidants, together with the four conventional antioxidents, and four kinds of RONS radicals. Through the structural optimizations and Mulliken charge analysis, we revealed that FeP(Im)-H 2 in the end-on asymmetric configuration was characterized by dihydrogen bonding with one H atom as donor and another H atom as acceptor, while the side-on symmetric configuration was featured by Kubas bonding with H 2 molecule acting as donor and acceptor. Our calculations of energetics indicate that the end-on asymmetric mode is more favorable for the H-H bond breaking with the lower activation barriers for adsorption and dissociation of H 2 than the side-on symmetric mode. We found that H-BDEs of FeP(Im)-H 2 and -H complexes were clearly lower than those of RONS-H complexes, indicating the feasibility of using these protoheme-hydrogen complexes as effective antioxidants for all the RONS radicals, whereas molecular hydrogen was revealed to reduce only hydroxyl radical. These findings support the decisive role of protoheme as an effective catalyst in RONS antioxidation by molecular hydrogen in vivo.