Accelerated discovery of superoxide-dismutase nanozymes via high-throughput computational screening

The activity of nanomaterials (NMs) in catalytically scavenging superoxide anions mimics that of superoxide dismutase (SOD). Although dozens of NMs have been demonstrated to possess such activity, the underlying principles are unclear, hindering the discovery of NMs as the novel SOD mimics. In this work, we use density functional theory calculations to study the thermodynamics and kinetics of the catalytic processes, and we develop two principles, namely, an energy level principle and an adsorption energy principle, for the activity. The first principle quantitatively describes the role of the intermediate frontier molecular orbital in transferring electrons for catalysis. The second one quantitatively describes the competition between the desired catalytic reaction and undesired side reactions. The ability of the principles to predict the SOD-like activities of metal-organic frameworks were verified by experiments. Both principles can be easily implemented in computer programs to computationally screen NMs with the intrinsic SOD-like activity.

Reactions S1a−S1e can be expressed with the Hess cycles as illustrated in Fig. 3b. The reactions involved in this Hess cycles and the corresponding changes of G are defined as follow, * (aq), ∆ sol * ° (S7) In eqs. S2 and S3a−S3e, ∆ r i° (i = 6−11) are the changes of standard-state G for the reactions; in eqs. S4b, S5b, and S6b, ∆ sol R * ° (R = HO, H, O) are the changes of G for the solvation of the adsorbates; in eqs. S4a, S5a, and S6a, ∆ ads R° (R = HO • , H • , O •• ) are the changes of G for the surface adsorptions of radical R in gas phase; in eq 7, ∆ sol * ° is the change of G for the solvation of the nanomaterial.

Calculation of ∆ r G º 6
The value of ∆ r 6° was calculated based on the following reactions: In eq. S12b, HO 2 •°(g) is the standard G of HO 2 • in gas phase; H 2 O°( aq), O 2 •−°(aq), and H 3 O +°(aq ) are the standard G for H2O, O 2 •− , and H3O + in water, respectively. These four values were obtained by frequency analysis calculations based on the optimized molecular geometries. The B3LYP method (1,2) in conjugation with the 6-311++G(d,p) basis set (3,4) as implemented in the Gaussian 09 program (5) was used for these calculation. The solvation effect of water was considered using the SMD model.

Calculation of ∆ r G º i (i = 7−11)
The values of ∆ r i° (i = 7−11) were calculated using the thermodynamic data taken from the chemistry handbook (6). The thermodynamic data that were used to calculate the ∆ r i° (i = 7−11) was shown in Supplementary Table 2. According to eqs. S1a−S1e, ∆ r 7°= 0. Estimation of ∆ sol,2 G º HO* , ∆ sol,2 G º O* , and ∆ sol,2 G º H* , ∆ sol,2 R * ° are defined in eqs. S9c-S9e. Such solvation energies are hard to determine either experimentally or computationally. For the purpose of estimation, the values of ∆ sol,2 HO * ° and ∆ sol,2 O * ° were supposed to be equal and their values on difference material surfaces were also supposed to be equal. Then, the values of ∆ sol,2 HO * ° and ∆ sol,2 H * ° could be estimated by the combination of Fig. 3d and available experimental results. According to the experimental results of Table 1, noble metal Pd and Pt had the SOD-like activities while defect-free CeO2 did not have the activity. Such information can be used to determine the border of the line corresponding to x1 = 0.5 in Fig. 3d, which in turn could be used to derive the values of ∆ sol,2 HO * ° and ∆ sol,2 H * °. Using this method, the values of ∆ sol,2 HO * °, ∆ sol,2 O * °, and ∆ sol,2 H * ° were estimated to be 1.2 eV, 1.2 eV, and 0.81 eV, respectively.

Calculation of electronic density of states
The geometry optimization of all nano-surfaces of Supplementary Fig.2 were performed using the projector augmented wave (PAW) method (7) implemented in the Vienna Ab initio Simulation Package (VASP) (8). The Perdew-Burke-Ernzerhof (PBE)(9) exchange-correlation functional with the generalized gradient approach (GGA) was applied. The cutoff energy of plane-wave was 400 eV. The vacuum values of 15 Å along z orientation were set to avoid interaction between units. The electronic and geometry optimization convergence criteria were set to 10 -5 eV and 0.02 eV Å -1 , respectively. For metallic nano-surfaces, such as Au(111), Ag(111), Pt(111), Pt(111), graphene, Nb2C and V2C, gaussian smearing with a width of 0.2 eV was used for the Fermi level. However, for semiconductive nano-surfaces, such as NiO(100), MoS2-x, Mn3O4(001), Co3O4(001), δ-MnO2, MnO(001) and CeO2(111), gaussian smearing with a width of 0.05 eV was used for the Fermi level. Electronic density of state (DOS) of all nano-surfaces were calculated in the level of Heyd-Scuseria-Ernzerhof (HSE06)(10) hybrid functional using the optimized structure in PBE + U GGA function. More calculated parameters, such as thickness and supercell size of slabs, K-points and Ueff values, one could see Supplementary Table 1.

Calculation of adsorption energies
Because only the most thermodynamically favorable adsorption sites are relevant, the adsorption energies of metals, metal oxides, and metal fluorides were calculated only for the sites that are most likely to have the strongest affinities for the corresponding adsorbates. The surface adsorption chemistry of these materials have already been well studied before. According to the known knowledge, the hydroxyl (HO·) group slightly prefers to adsorb on the "hollow" or "bridge" sites on metal surfaces, where the oxygen can coordinate with more metals (11). So does the hydrogen atom (11). This might be understood by that both hydroxyl and hydrogen adsorbates accept electrons from the surrounding metals, which subsequently have electrostatic attractions with the surrounding metals. So, both adsorbates prefer to form the closely-packed configurations with their surrounding metals to maximize the nonmetal-metal coordination numbers. As for metal oxides and fluorides, the hydroxyl group prefers to adsorb on the surface metals because these metals have less coordination numbers than metals in the bulk and are coordinatively unsaturated. In contrast, the hydrogen atom prefers to adsorb on the O or F sites on metal oxides and fluorides, because atomic hydrogen is highly reductive and prefers sites with strong electron affinity (e.g., atoms with strong electronegativity). For these reasons, the hydroxyl and hydrogen adsorption energies for all metal surfaces were calculated based on the computationally relaxed adsorption configurations where the adsorbates were located near the hollow sites of the surfaces. The hydroxyl adsorption energies for all metal oxides and fluorides were calculated based on the configurations where the hydroxyls were near the surface metals and the hydrogen adsorption energies were calculated based on the configurations where the hydrogen atoms were near the O or F sites in the surfaces. All these relaxed adsorption configurations could be readily obtained by geometry optimizations not using symmetry constraints and using the initial structures where the groups were placed near the desired sites.
MIL-53(Fe)-H contains many possible addition sites for hydroxyl and hydrogen radicals. However, only several are probable for the side reactions illustrated in Fig. 3a. In the MIL-53(Fe)-H structure, each iron is coordinated with five oxygen atoms; it has a free d-orbital for further covalent bonding. Therefore, irons in the MIL-53(Fe)-H structure can well model irons in the surfaces of the MOF materials. Except for irons, all the other atoms in the MIL-53(Fe)-H structure are covalently saturated. So, irons are the only probable sites for the addition of hydroxyls on MIL-53(Fe)-H. Atomic hydrogen is a highly reductive racial, which prefers to adsorb on atoms with strong electronegativity. So, oxygen atoms are the most probable sites for the additions of hydroxyls on MIL-53(Fe)-H. Although the additions of hydroxyls and hydrogens to other sites of MIL-53(Fe)-H are possible to yield products with even larger thermodynamic stabilities, these additions severely destructed the structure of MIL-53(Fe)-H and did not need to be considered because they were unlikely the side reactions of O2 ·− . The one-dimensional nanowire structure was used to model the MOF structure (e.g., see Supplementary Fig.6). Supplementary Fig. 7 comparably shows the adsorption energies of H on different sites. As can be seen from Supplementary Fig.7, oxygen bridging two irons is indeed the most favorable site for hydrogen.

Energies vs. hydrogen electrode (HE) potential
The energy levels of a material surface with respect to HE potential (Evs HE) were calculated using the following procedures: 1) Calculating the electronic band structure for the material using the slab model, which obtains energy bands (Evs fermi) with respect to Fermi energy (Efermi) for the material; 2) Calculating vacuum energy (Evac) for the NM; 3) Calculating Evs HE using this equation, Evs HE = Evac -Evs fermi -Efermi -4.5 in which 4.5 eV is the scaling factor relating the normal hydrogen electrode scale (NHE) to absolute vacuum scale (AVS);

Characterization of ceria particles
Diffuse reflectance ultraviolet-visible (UV-vis) spectra were recorded using an Agilent Cary 5000 ultraviolet-visible-infrared spectrometer (BaSO4 as a reference). Valence-band X-ray photoelectron spectra were given by an ESCALab220i-XL spectrometer equipped with a twinanode Al Kα (1486.6 eV) X-ray source.

Calculation of optical energy band-gap for ceria particles
The optical band-gap (Eg) of samples was determined according to the following equation, where hν is the photon energy (J), A is a proportionality constant, and n is an index that indicates the nature of transition (n = 1/2 for indirect transition; n = 2 for direct transition). Considering that CeO2 nanoparticles are claimed to be an indirect band-gap material (n = 2), the recorded UV-vis spectra were first transformed to the curve of (αhν) 2 versus ℎ and the variations of (αhν) 2 versus ℎ were plotted. The straight line range of these plots was extrapolated to the x-axis (y = 0) to obtain the values of Eg.

Supplementary Note 1. Work function (W f ) version of the energy level principle for metals.
For metal surfaces, the Fermi energies with respect to the vacuum energy are known as Wf. So, the Wf-version of energy level principle can be obtained: φ1 + 4.5 < Wf < φ2 + 4.5 In the above equation, 4.5 V is the difference between vacuum energy and the hydrogen potential. The Wf of the elemental substances available in the Lange's Handbook of Chemistry(6) are plotted in Supplementary Fig. 3b. As can be seen, Au, Pt, and Pd have their Wf located in the range. Because Wf is usually measured in an ultrahigh vacuum condition, which is quite different from the aqueous condition for the substances to function as SOD mimics. Unlike noble metals, which are chemically inert, many substances of Supplementary Fig. 3b easily react with O2 or water at room temperature, which would dramatically change their Wf. So, the result of Supplementary Fig.  3b cannot always be directly used to predict whether the elemental substances are potential SOD mimics.

Supplementary Note 2. Verification of the energy level principle with ceria particles.
Recent experimental study has demonstrated that ceria particles synthesized under the same condition except with at different temperatures, 0, 30, 60, and 90 ℃ (denoted as Ceria_0, Ceria_30, Ceria_60, and Ceria_90, hereafter) have markedly different SOD-like activities: Ceria_0 and Ceria_30 have much stronger SOD-like activities than Ceria_60 and Ceria_90 [see Fig. 3A of Ref. (12)]. To check whether these ceria particles obey the energy level principle, the frontier molecular orbitals of these ceria samples were measured by the combination of experiments and calculations. To this end, the UV−vis diffuse reflectance spectra (DRS) of these ceria particles were measured ( Supplementary Fig. 8). These DRS spectra were then transformed to the (αhν) 2 versus hν plots to obtain the gap between valence band and conduction band and the gap between valence band and defect levels ( Supplementary Fig. 9). The energies from the valence bands to the Fermi levels were measured using the XPS valence band spectra (Supplementary Fig. 10). Because the energy levels (versus hydrogen electrode, HE) of conduction band (ECB, V) and valence band (EVB, V) of ceria particles cannot be measured using the ultraviolet photoelectron spectroscopy because of poor electron conductivity of the samples, these energies were calculated using the method described in the above Supplementary Method section. The results are shown in Supplementary Fig. 11. As can been seen, ceria_0 and ceria_30 rather than ceria-60 or ceria_90 have the iFMO, which agrees with their SOD-like activity order and verifies the energy level principle.

Supplementary Note 3. Computational support for the predicted SOD-like activities of the 2D materials.
To our knowledge, 14 of the 121 2D materials have been experimentally realized. These 14 materials include nine transition metal dichalcogenide (TMDC) structures. All these nine TMDCs have been chosen for computational investigation to check their SOD-like activities. The CBM and VBM energies of these nine TMDC structures are shown in Supplementary Table 3. As can be seen from this table, six structures including ZrS2, PtS2, ZrSe2, ZnS2, HfS2, and VS2 have their CBMs located in the range (−0.16 eV, 0.94 eV). According to our energy-energy principle, they should follow the LUMO-mediated mechanism to mimic the activity of SOD. Similarly, the remaining three structures including MoTe2, WSe2, and WTe2 have their VBMs in the range and thus they should follow the HOMO-mediated mechanism. To verify this prediction, we have located the intermediates and transition-state structures involved in both the HOMO-and LUMOmediated pathways for ZrS2 and MoTe2 and plotted the reaction energy profiles in Supplementary  Fig. 12. As shown in Supplementary Fig. 12a, the energy barrier of the rate-determining step (RDS) of the LUMO-path is much lower than that of the RDS of the HOMO-path (0.83 eV vs 2.38 eV), suggesting that ZrS2 indeed prefers the LUMO-mediated mechanism, in good agreement with the prediction. On the other hand, Supplementary Fig. 12b suggests that MoTe2 prefers the HOMOmediated mechanism, which also agrees with the prediction. For the remaining seven TMDC structures, their kinetically favorable pathways were investigated by locating only the intermediate structures. As shown in Supplementary Fig. 13, the LUMO-mediated paths for HfS2, SnS2, PtS2, ZrSe2, and VS2 do not contain any intermediate structures with relatively high energies, in agreement with the prediction that these structures prefer the LUMO-mediated mechanism. Similarly, the results of Supplementary Fig. 14 suggest that the HOMO-mediated mechanisms for feasible for WTe2 and WSe2, in agreement with the prediction. Therefore, these computational results have supported the SOD-like activities predicted by the high throughput calculations.
Supplementary Fig. 1. Reaction energy profiles indicating the different mechanisms of the NMs. a The surface of NiO(100) prefers the mechanism (red) in which the surface is first oxidized (forming the H2O2) and then reduced back (forming the O2); the mechanism (black) is kinetically disfavored. b The surface of MoS2−x prefers the mechanism (red) while the mechanism (black) is less kinetically competitive. In A, structure labeled with SA was used as the starting point for the reaction energy profile of Fig. 2f.  In the left panel, the H atom was adsorbed on the oxygen bridging two irons and made a hydrogen bond with the neighboring bridging oxygen. In the middle panel, the H atom was also adsorbed on the oxygen bridging two irons but made a hydrogen bond with the ligand oxygen. In the right panel, the H atom was adsorbed on the iron. In all these panels, the iron and adsorbed hydrogen atoms were shown in the ball-and-stick mode and all the other atoms in the stick mode. Relative energies of these adsorption configurations were marked.
Supplementary Fig. 8. UV−vis diffuse reflectance spectra (DRS) of CeO 2 _0, CeO 2 _30, CeO 2 _60, and CeO 2 _90. These DRS spectra were used to obtain the gap between valence band and conduction band and the gap between valence band and defect levels for the samples, and the details can be found in the Supplementary Notes.