Influence of Varying Functionalization on the Peroxidase Activity of Nickel ( II ) –Pyridine Macrocycle Catalysts: Mechanistic Insights from Density Functional Theory

: Nickel(II) complexes of mono-functionalized pyridine-tetraazamacrocycles (PyMACs) are a new class of catalysts that possess promising activity similar to biological peroxidases. Experimental studies with ABTS (2,2 ′ -azino-bis(3-ethylbenzothiazoline-6-sulfonic acid), substrate) and H 2 O 2 (oxidant) proposed that hydrogen-bonding and proton-transfer reactions facilitated by their pendant arm were responsible for their catalytic activity. In this work, density functional theory calculations were performed to unravel the influence of pendant arm functionalization on the catalytic performance of Ni(II)–PyMACs. Generated frontier orbitals suggested that Ni(II)– PyMACs activate H 2 O 2 by satisfying two requirements: (1) the deprotonation of H 2 O 2 to form the highly nucleophilic HOO − , and (2) the generation of low-spin, singlet state Ni(II)–PyMACs to allow the binding of HOO − . COSMO solvation-based energies revealed that the O–O Ni(II)–hydroperoxo bond, regardless of pendant arm type, ruptures favorably via heterolysis to produce high-spin ( S = 1) [(L)Ni 3+ –O · ] 2+ and HO − . Aqueous solvation was found crucial in the stabilization of charged species, thereby favoring the heterolytic process over homolytic. The redox reaction of [(L)Ni 3+ – O · ] 2+ with ABTS obeyed a 1:2 stoichiometric ratio, followed by proton transfer to produce the final intermediate. The regeneration of Ni(II)–PyMACs at the final step involved the liberation of HO − , which was highly favorable when protons were readily available or when the p K a of the pendant arm was low. abstract a proton either from H 2 O or H 2 O 2 . The positive charge of Schiff base rendered it susceptible to by free basic species, leading to new intermediates that were identical to INTMD1 of 2B and 3B . Such a case would turn out different if calculations were made in aqueous solution. Water molecules will stabilize the iminium group by electrostatic attraction and hydrogen bonding.


Introduction
Oxidation reactions involving electron transfer are customary in the study of biological processes, in organic synthesis, and in the production of novel energy sources, among others [1]. While reactions of this type are essentially fast, the rate of electron transfer between two redox species may be hastened if a more favorable reaction pathway is conceived in the presence of a catalyst. For reactions involving hydrogen peroxide (H2O2) as an oxidant, typical catalysts usually involve transition metals (e.g., Fenton reagent: Fe 2+ + H2O2) [2][3][4][5]. From the studies reported to date, two general types of mechanisms for H2O2 activation have been proposed: (1) the radical mechanism that involves HO· and HOO· that both create a redox cycle with the metal species [6]; and (2) the peroxide complex mechanism, which causes a two-electron oxidation of the metal, leading to the heterolytic cleavage of the peroxide O-O bond [7]. The activation of H2O2 commences when the oxidant coordinates to the metal center of the catalyst via nucleophilic addition [8,9]. Its rate of binding to the complex is affected by the nature of the metal, the ligands, and the solution matrix [10]. Once the oxidant is bound, those complexes that contain metals with a single stable non-zero oxidation state (e.g., Al, Be, Mg, Ca, Zn, Cd, etc.) will not participate in the electron transfer reaction, although they function as Lewis acids. In contrast, complexes of transition metals that have at least two relatively stable oxidation states (e.g., Fe, Ni, Co, Cu, Mn, etc.) are typical for the Fenton-type chemistry and will serve as a source and sink of electrons during the reaction [11].
Peroxidases are among the natural enzymes that have high efficiency and specificity for H2O2-driven oxidation of substrates [12][13][14]. However, similar to other enzymes, peroxidases are expensive and are prone to denaturation in response to extreme conditions [15][16][17]. With such a disadvantage comes the initiative to design cheap metal-based biomimetics with improved catalytic sustainability over a wide range of pH and temperatures [18]. The most prominent of which are iron-based complexes, yet many other compounds containing metals such as Mn, Cu, Co, Ru, V, and Ti have also been reported [9,[19][20][21][22][23]. In 2009, Organo and colleagues synthesized novel Ni(II) complexes of monofunctionalized pyridine-tetraazamacrocycles (PyMACs) (Figure 1), which exhibited catalytic activity toward the one-electron oxidation of 2,2′-azino-bis(3-ethylbenzo thiazoline-6-sulfonic acid) (ABTS) with H2O2 as the oxidant [24]. These catalysts, bearing different pendant arms, displayed remarkable reaction kinetics where they accelerated the production of ABTS +· at varying rates. Ni(II)-PyMACs of acidic (pKa < 7) pendant arms were observed to possess higher activity than Ni(II)-PyMACs with a basic functional group. Parallel findings were also published for copper complexes of analogous macrocycle structures, where the pentadentate carboxylic acid complex demonstrated mild peroxidase-like activity, while the unfunctionalized one showed no activity at all [25]. So far, among the assumptions highlighted in the mechanism were the possible interactions between the coordinated peroxide and the pendant arm via hydrogen bonding and intramolecular proton transfer. This scenario has already been implicated in carboxyl-appended "hangman" metalloporphyrins in which H2O2 is activated through proton-coupled electron transfer that facilitates the cleavage of its O-O bond [26]. While the study on Ni(II)-PyMACs provided preliminary evidence of the influence of functionalization on their catalytic activity, no succeeding reports have ever investigated the reaction mechanism by which the pendant arm group facilitates the oxidation of the substrate in the presence of H2O2.  In the present work, insights into the role of Ni(II)-PyMAC pendant arm in the oxidation of ABTS by H2O2 are presented. Density functional theory (DFT) calculations were employed to model the reaction systems [27][28][29][30][31][32]. The nature, energy, and reactivity of Ni(II)-PyMACs were examined together with the potential intermediates and transition states suggested by the collective information reported in the literature [8,13,[33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. The selection of the most plausible reaction mechanism was based primarily on the frontier orbitals of the involved structures, their relative stabilities, the estimated O-O bond dissociation energies, and the ease of electron transfer from ABTS to active oxidative species using energy levels of highest occupied molecular orbitals (HOMO) and lowest unoccupied molecular orbitals (LUMO). The general pendant arm (L) of all hypothetical structures presented in this paper was treated neutral in charge unless stated otherwise.

Proposed Reaction Pathways
The mediated one-electron oxidation of ABTS by Ni(II)-PyMAC catalysts in the presence of H2O2 was initially hypothesized to follow either of the two reaction pathways (A and B) shown in Schemes 1 and 2, respectively. The pendant arm of the catalyst is proposed to undergo a reversible proton dissociation to form the acid-base pair. The key primary step of the catalytic cycle, regardless of the form of the catalyst, is the binding of H2O2 to the nickel complex either as is or in the form of the highly nucleophilic HOO − species made possible by the presence of an abstracting base. Either way, the O-O bond of the first intermediate undergoes cleavage (either homolytic or heterolytic) to generate the active nickel-oxo species. The four spin multiplicities (i.e., S = 0, 1/2, 1, and 3/2) of nickel-oxo are all accepted as candidate forms of the active intermediate as necessary. Regardless of the spin state, ABTS undergoes one-electron oxidation via electron transfer to the nickel-oxo species, which may or may not be coupled with a proton, depending on the pathway followed. This step subsequently leads to the formation of ABTS + and nickel(II)-hydroxo complexes.
To capture the intrinsic properties and the possible electronic variations of structures in response to solvation, separate calculations were performed in the gas phase and in solution. The Conductor-like Screening Model (COSMO) [65][66][67], as implemented in ADF, was utilized to calculate energies in solution where each solute structure was embedded in a molecule-shaped cavity of the chosen dielectric medium. The dielectric constant and radius of the rigid sphere were set to 78.39 and 1.93 Å, respectively to describe water as the solvent.
The initial geometries of all molecular models, except for the transition states, were fully optimized using the default integration accuracy and convergence criteria of the software. For systems with hard-to-converge SCF, corresponding occupation numbers were specified in the input to make sure Aufbau solutions were obtained. SCF convergence criterion was also tightened to 1 × 10 −5 from the default of 1 × 10 −3 . The optimized spin states of the structures with the lower electronic energies were taken as their electronic ground states, although all spin multiplicities were considered for the succeeding vibrational frequency calculations. To verify predicted ground states of structures with closely lying states, single-point calculations using BLYP-related hybrid functionals (B3LYP and B3LYP*) were performed on the optimized geometries.
For the transition states, linear transit (LT) calculations were first carried out, where one or a combination of the three LT parameters (i.e., distance, angle, dihedral angle) were set and defined in the input, getting initial and final values. All coordinates were frozen to run a sequence of single-point calculations at equidistant steps along the LT path. The one with the highest energy among all geometries was selected to undergo the final and true transition state search. The same convergence criteria for both the geometry and SCF cycles were applied in the transition state search as those defined in geometry optimization runs. The nature of all the stationary points generated was verified through vibrational frequency calculation.
After obtaining all the optimized structures from geometry optimization and transition state search, the generated Hessian matrices were calculated analytically by doing vibrational frequency. One-point numerical differentiation, symmetric displacements, and medium integration accuracy were applied purposely to minimize the computational time. Imaginary frequencies were checked for all optimized structures to confirm whether the stationary points are local minima (equilibrium geometries) or first-order saddle points (transition states). Statistical thermal analysis (i.e., to calculate for nuclear internal energy and entropy) was carried out at a constant temperature and pressure of 298.15 K and 1 atm, respectively. The corresponding zero-point vibrational energies (ZPE) were derived automatically from the measured frequencies and were included in all quoted energies.

Structural Characterization of Ni(II)-PyMAC Complexes
Fully optimized structures of mono-functionalized Ni(II)-PyMAC complexes are shown in Figure 2, where both the acid and base forms of each catalyst were studied separately, except complex 4. The decision to investigate both forms of the catalysts was influenced by previous findings that the protonation of 5B to 5A reduced the catalytic strength of the original complex during the first few seconds of the reaction but made it acquire a much longer shelf-life [24]. These notable differences between this acid-base pair prompted us to consider that the functional state of the pendant arm affects the catalytic action of these complexes. Both forms must be treated separately to be able to elucidate several possible mechanisms of action for these catalysts.
The optimized structures of complexes 1 to 6 supported the idea that pendant arms have different binding modes toward the nickel metal. The lowest energy conformation of the three basic complexes, 1B, 2B, and 3B, adopted a five-coordinate square pyramidal geometry, while their protonated counterparts favored the slightly distorted, four-coordinate square planar conformation. Deprotonation of the pendant arms increased their nucleophilic character and allowed them to coordinate readily to the Lewis acidic metal at the axial position. In fact, with the strong sigma-donor interaction in these three complexes, specifically between the orbital of the axial ligand and the orbital of the nickel metal, the separation gap between their respective HOMOs and LUMOs was small, which allowed them to achieve a high-spin (S = 1) electronic ground state. In contrast, optimization of the structures without available axial ligands, such as 4, all the protonated forms, and 6B, only led to their low-spin (S = 0) states, as expected. These results regarding the ground state of the complexes were in agreement to what was previously reported [24], except for 5B, which after optimization in both spin states emerged to be more stable in the singlet rather than in the triplet form by a small free energy difference of 0.19 kcal mol −1 . We can account for such nonconformity to the fairly weak Ni-Naxial interaction in 5B brought by a lesser nucleophilic character of neutral amine relative to other basic pendant arms. As a result, the energy gap between the low-spin and high-spin states is small, to which the BLYP functional incorrectly predicted the ground state. This incongruity with experimental findings is not peculiar in most GGA functionals where the low-spin multiplicity is often over-stabilized when calculating weak-field complexes [68][69][70]. Fouqueau et al. (2004) believed that such a phenomenon is caused by the underestimation of the Fermi correlation and is normally resolvable by the inclusion of the Hartree-Fock (HF) exchange parameter into hybrid functionals [71]. The energy splitting was shown to depend linearly on the coefficient of exact exchange admixture and that hybrid functionals provided much more reasonable values [72,73]. Thus, at this point, we were prompted to calculate the bonding energy of 5B at multiple spin states using BLYP and two hybrid functionals -B3LYP and B3LYP* to check for possible discrepancies in the resulting ground state. Based on the obtained values (Table 1), high-spin/low-spin state separation is more pronounced when employing functionals with a higher percentage of HF. B3LYP (with 20% exact exchange) resulted in a triplet state that is 11.5 kcal mol −1 more stable than its corresponding singlet pair. Similarly, B3LYP* also correctly predicted the ground state of 5B with a 7.7 kcal mol −1 energy gap between the two spin multiplicities. Recognizing this as a limitation of GGA BLYP, we eventually accepted the high-spin (S = 1) to be the ground state of 5B consistent with the previous experiment. All calculated energies from the statistical thermal analysis of fully optimized Ni(II)-PyMAC catalysts (in both spin states) are summarized in Table 2.  Table 1. High-spin/Low-spin bond energy splitting (kcal mol −1 ) in 5B complex as derived from single-point calculations at selected exchange-correlation (XC) functionals. A positive value in energy splitting indicates that the predicted ground state is high-spin.

XC Functional HF Exchange Admixture Bond Energy Splitting
0.20 11.5 Table 2. Calculated energies in the gas phase of optimized Ni(II)-PyMAC structures in their low-spin (S = 0) and high-spin (S = 1) states at 298.15 K and 1 atm. U, H, S, and G are the internal energy, enthalpy, entropy, and Gibbs free energy, respectively, corrected with ZPE (zero-point vibrational energy).  (Table 3) were not much different [24,74] (refer to Figure 3 for the atomic labels). There were only slight elongations in Ni-Ni bond lengths ranging from 0.039 to 0.083 Å for 1A, 0.050 to 0.068 Å for 6B, 0.032 to 0.138 Å for 6A, and 0.035 to 0.086 Å for 5BHS, and bond length shortening between 0.047 and 0.101 Å for 2A. The Ni metals in the optimized models were also slightly undisplaced from the ring except only for 2A, which showed a 0.90° constriction of the N1-Ni-N3 bond angle due to the considerable protrusion of Ni above the equatorial plane. Nonetheless, the results of the optimization were consistent throughout all types of Ni(II)-PyMAC complex, where the four Ni-Nequatorial bonds were longer in the deprotonated forms than in their conjugate acid counterparts. The pendant arms in these forms could favorably act as axial ligands and slightly draw the nickel metal out of the plane. Table 3. Bond lengths (Å) and angles (°) for the optimized electronic ground states of protonated and deprotonated forms of Ni(II)-PyMAC complex cations. Experimental pKas are also shown.

H2O2 Activation and Reactivity with Ni(II)-PyMACs
The binding of H2O2 to the complex is a crucial step in the catalytic process of Ni(II)-PyMACs.  2+ , the coordinate covalent bond between Oaxial and Ni elongated to the extent that it eventually disappeared during optimization. The presence of the unabstracted proton made the terminal oxygen of H2O2 weak to coordinate axially to Ni. In fact, its energy level diagram consistently showed an unoccupied orbital below HOMO, which suggested that the structure was not in the ground state when calculated in high-spin. In contrast, geometry optimization of high-spin [(LMe)Ni 2+ -OOH] + resulted in a distinct binding of hydroperoxo with Ni and followed Aufbau occupations in its energy-level diagram. While Ni(II)-PyMACs can take on a slightly distorted square planar geometry, it can be assumed that HOO − can coordinate axially to Ni either above (b1) or below (b2) the equatorial plane of the macrocycle. To determine which side is more susceptible to binding, we optimized [(LMe)Ni 2+ -OOH] + with HOO − positioned below and compared its energy with that of the overhead-coordinated HOO − . Based on the bonding energies in the gas phase, the nucleophilic attack of hydroperoxo ion to Ni is more favorable above than below the plane by approximately 2.6 kcal mol −1 . These close energies between the two forms suggest that HOO − can attack the nickel metal either way. However, the top-axial coordination is preferred, since Ni tends to protrude above the plane due to the presence of the pendant arm. As Ni gets more exposed above than below, the preference of HOO − to bind at this location also increases. With this result, we took the "overhead" [(L)Ni 2+ -OOH] + as the most plausible form of INTMD1 across all catalysts, and the attempt to bind HOO − anti to the axially bound pendant arm at this point was no longer made.
HOMO orbitals of optimized Ni(II)-PyMACs in the gas phase (Figure 4) display the portions of the complexes that are likely to interact with H2O2 before it binds to Ni. Note that for complexes with coordinated pendant arms, HOMO for both the high-spin and low-spin structures are shown. Based on the shape and location of these orbitals, we can group the 11 catalysts into four: (G1) 1A, 4, 5A, and 6A; (G2) 2A, 3A, and 6B; (G3) 1B and 5B; and (G4) 2B and 3B. G1 catalysts are those that have HOMO localized on the nickel center, specifically of the orbital type , which is in contrast to G2 catalysts that have lobes mainly confined in the pendant arm. On the other hand, low-spin G3 and G4 catalysts have their highest-energy electrons occupy the * − orbitals, but they differ in their high-spin states where G3 HOMO is of the type * − / , while that of G4 is * − .
Using this classification (Table 4), the catalysts of the same group will follow the same mechanism in activating H2O2 toward binding to the complex. There are two requirements for this step to proceed: (1) H2O2 should first be deprotonated by a base to produce the more nucleophilic HOO − , as discussed earlier; and (2) the top-axial position of Ni in the complex should first be free and available for the consequent binding of HOO − . The first requirement renders G2, G3, and G4 catalysts the advantage over G1 in forming INTMD1, since they can abstract a proton from H2O2 using the electrons (HOMO) of their pendant arms. In terms of this ability, G4 catalysts are the strongest since they are the most basic, followed by G3 and then G2. G1 catalysts, in contrast, would require the presence of any other base just to deprotonate H2O2 before binding to them. Meanwhile, the second requirement implies that high-spin complexes specifically G3 and G4 should have their pendant arms "off" of Ni before the oxidant can bind either below or above. G3 catalysts in their high-spin cannot favorably react with H2O2 due to the absence of proper symmetry and shape of HOMO, which requires them to turn to low-spin states before they can deprotonate H2O2 by a significant overlap of their orbitals. High-spin G4 catalysts in comparison to G3 can meet the first condition but not the second. The pendant arms of G4 hardly dissociate from Ni and thus, HOO − cannot bind easily on top.
In fact, high-spin G4 HOMO can overlap sufficiently with LUMO of H2O2 so that they do not need to detach their amide pendant arms from Ni to deprotonate the oxidant. This makes them the least effective catalysts among the four groups being classified. The prerequisite to change from high-spin to low-spin before HOO − can coordinate to the complex is further supported by the HOMO ( ) orbitals of the catalysts, in which all of them provide the same molecular platform only when they are in low-spin. Although coordination to Ni below the plane is also possible, it cannot render any of the catalysts better than the rest as they are all capable of the same. Thus, it can be inferred that at the first reaction step of the mechanism, G3 complexes may be the most effective, followed by G2, then G1, and lastly, G4.  In effect, four-coordinate complexes rather achieve HOMO similar to G1 of the previous classification, and five-coordinate complexes acquire * − / similar to G3. 6B is a special case given that hydrogen bonding with N is forbidden due to the steric hindrance and electron-donating abilities of two ethyl groups. The HOMO of 6B is the same for both cases. Nonetheless, the addition of a water cavity generally masks the ability of some pendant arms to possibly react with H2O2. For this step, gas-phase calculations were better at differentiating the pendant arms in terms of their reactivity with H2O2.

Structure and Energetics of Nickel-Hydroperoxo Species
Nickel-hydroperoxo intermediates were found identical for complexes that form an acid-base pair, such as in 1A/1B, 5A/5B, and 6A/6B. However, for those with amide pendant arms (e.g., 2A/2B and 3A/3B), the pairs did not immediately produce identical nickel-hydroperoxo species right after the first reaction step. The LUMO of the first intermediate of 2A and 3A displayed an orbital type different from those of other catalysts' INTMD1s. Orbital lobes of their LUMO were strictly confined in the carbonyl-amide functional group of their pendant arms (Scheme 4). By shape and symmetry, these orbitals were * in nature. Their pendant arms rather formed a protonated Schiff base as a result of nitrogen's lone pair forming a pi-bond with carbonyl carbon. This can shuttle the oxygen to abstract a proton either from H2O or H2O2. The positive charge of Schiff base rendered it susceptible to deprotonation by free basic species, leading to new intermediates that were identical to INTMD1 of 2B and 3B. Such a case would turn out different if calculations were made in aqueous solution. Water molecules will stabilize the iminium group by electrostatic attraction and hydrogen bonding.  (1A, 4, 5A, and 6A), whose reactivity also requires the assistance of H2O or any free base. Nonetheless, the 2A and 3A reaction mechanism still passes through the formation of INTMD1, which is a prerequisite step in forming the active intermediate. For

O-O Bond Cleavage: Heterolytic vs. Homolytic
Following the coordination of HOO − to the complex, the O-O bond of the generated nickel-hydroperoxo complex is proposed to undergo cleavage via either of the two pathways shown in Scheme 5. Many studies about oxidation catalysts similar to Ni(II)-PyMACs have reported contrasting evidence regarding which process is more favored [8,38,[40][41][42][43]. In this study, we first identified the possible intermediates arising from each of the two pathways.  [75], ruling out the formation of [(L)Ni 3+ =O] + . The energy of the high-spin (S = 3/2) structure was also calculated to be approximately 8 kcal mol −1 less than that of the low-spin (S = 1/2) structure. This prediction was consistent across all types of pendant arm and irrespective of whether calculations were made in the gas phase or in solution.  The gas-phase energies of nickel-oxo species (with consideration of HO▪/HO − energies) showed that the heterolytic pathway leading to [(L)Ni 3+ -O·] 2+ is thermodynamically forbidden by over 171 kcal mol −1 relative to the homolytic pathway that produces high-spin [(L)Ni 2+ -O·] + . This inference is synonymous with the findings of Sankaralingam et al. (2014), who also investigated the Ni(II)-catalyzed alkane hydroxylation with m-CPBA using similar DFT methods [76]. They reported [Ni 2+ -O·] + species to be the active intermediate of their catalysts, being more energetically stable than [Ni 4+ =O] 2+ by at least 900 kJ mol −1 (or approximately 215 kcal mol −1 ). This looked acceptable at first and supportive of our initial findings, but our calculations with solvation effects suggested otherwise. Due to solvation, the heterolytic pathway emerged to be more favorable with energy splitting of the resulting species (considering HO·/HO − energies) falling within 3 to 9 kcal mol −1 depending on the pendant arm. Hydration stabilized both homolytic and heterolytic intermediates, but the change in energy was more significant in the latter (e.g., ∆EHeterolytic = −131 kcal mol . This implies that after O-O bond heterolytic cleavage, Oaxial received an electron from Ni, causing the metal to achieve a +3 oxidation state. Specifically, the unpaired electron of Ni occupying the orbital * − / (Ni-Nequatiorial interactions) was utilized for this process. This is supported by the decrease in the total spin density of Ni + 4Nequatorial atoms by 0.94. Regardless of the pendant arm, [(L)Ni 3+ -O·] 2+ showed the same type and order of SOMOs in terms of energy, none of which is the highest lying orbital. Therefore, further oxidation by electron removal is highly unfavorable in these intermediates. Instead, [(L)Ni 3+ -O·] 2+ appear to be good oxidizing agents. They can readily accept two electrons from substrates given that they have two lowest-lying unoccupied orbitals that are nearly degenerate. The two SOMOs with corresponding energies (in eV) of the active (S = 1) [(L)Ni 3+ -O·] 2+ intermediates are shown in Figure 5.  Table 5) and their interaction. The x, y, and z axes point along N(1)-N(3), N(4)-N(2), and Ni-Oaxial directions, respectively (numbering of nitrogen atoms is similar to that in Figure 3).

Dissociation Energy of Nickel-Hydroxo O-O Bond
The ease in breaking the O-O bond is a factor that affects the degree of reactivity of H2O2 as oxidant. In the absence of any catalytic influence, the rate of peroxide bond dissociation may be very slow under mild conditions. From our calculations (Table 6), it would require 53.0 kcal mol −1 of energy in aqueous solution to produce two HO· radicals from H2O2. The transition state for this process verified by a single imaginary frequency (vi) of −254 cm −1 was also determined to lie 11.2 kcal mol −1 downhill of the resulting radical pair. These results agree well with the published report of Sandhiya and Zipse in 2017 [77]. Compared to that of [(L)Ni 2+ -OOH] intermediates, the H2O2 peroxo bond appeared to be the strongest. The energy barrier associated with the heterolytic breaking of the [(L)Ni 2+ O-OH] peroxo bond was only 26 to 27 kcal mol −1 . The energies of the resulting intermediates were also 21.6 to 28.3 kcal mol −1 lower than their corresponding INTMD1s. Thus, we can ascribe the catalytic activity of Ni(II)-PyMACs to their potential role of lowering the bond strength of the attached peroxo group to generate the active intermediate that will serve as the actual oxidant of ABTS. Based on calculated energy barriers, the ease of O-O bond breaking in [(L)Ni 2+ -OOH] follows the trend: 4 > 3A/3B > 2A/2B > 5A/5B > 6A/6B > 1A/1B. Unexpectedly, 1A/1B appears to be the poorest pair of catalysts in this step, while in the previous experimental study, it was shown to be the most effective [24]. Apparently, the energy difference between 1A/1B and 4 is not significantly large (approximately 1.8 kcal mol −1 ), and to treat the catalytic strength of the complexes this way may not be the most ideal. Therefore, we also looked into each catalyst's degree of partitioning between heterolysis and homolysis. Although homolytic intermediates are less stable in aqueous solution, it does not guarantee that the [(L)Ni 2+ -O-OH] + bond may not break evenly, forming a hydroxyl radical and [(L)Ni 2+ -O·] + . Nam et al. (2000) had previously pointed out that both pathways occur in iron porphyrin complexes, and the preference of the complex is affected by its electronic nature [42]. It was emphasized that electron-deficient iron porphyrin showed a tendency to cleave the hydroperoxide O-O bond heterolytically, whereas electron-rich ones cleave the same bond homolytically. As shown in Figure 6, all catalysts will favor the homolytic pathway if reactions are carried out in the absence of any solvent. However, the opposite is true for catalysis occurring in aqueous solution. Heterolysis is more preferred since water will contribute to the stabilization of charged species that are produced after bond cleavage. Among the six types, 5A/5B and 6A/6B showed the smallest partitioning between the two pathways in aqueous solution, followed by 4, 1A/1B, 2A/2B, and 3A/3B of relatively the same disparity. The pendant arms in 1A/1B, 2A/2B, and 3A/3B are all electron-withdrawing due to the presence of a carbonyl group. As a result, the nickel metal in these complexes is more electron-deficient compared to that in 5A/5B and 6A/6B. Then, this greatly weakens the O-O bond in [(L)Ni 2+ -OOH] + , thus increasing its chance to break heterolytically. Although for all catalysts the homolytic cleavage may still proceed regardless of pendant arm type, a partition gap between the two processes would indicate that preferably more of [(L)Ni 3+ -O·] 2+ and the HO − pair are produced after O-O bond cleavage rather than [(L)Ni 2+ -O·] + and HO·.

ABTS Oxidation by [(L)Ni 3+ -O·] 2+ Active Intermediates
ABTS in its fully deprotonated di-anionic form (HABTS + pKa = 2.2) [78] was optimized in solution and was verified to have a singlet electronic ground state (E = −6865.0 kcal mol −1 ). Its highest energy electrons were found occupying the antibonding * molecular orbital delocalized around the aromatic rings up to the bridge nitrogens. From this, it can be expected that in the presence of an oxidizing species (e.g., [(L)Ni 3+ -O·] 2+ or HO·), the first electron that will be removed from ABTS will come from this region (Figure 7). The computed ground state (S = 1/2) of ABTS +· (E = −6762.5 kcal mol −1 ), the one-electron oxidized form of the substrate supports this assumption. It was noted that the HOMO and LUMO obtained from this structure were the same in shape, although each has a different energy ( HOMO = −4.979 eV and LUMO = −4.650 eV). In addition, the obtained spin density values indicate that the unpaired electron of ABTS +· is distributed along the two central nitrogens ( Ncenter = 0.198) and the two other ring nitrogens ( Nring = 0.111).  The energy difference between the HOMO of ABTS and the LUMO of each oxidant was calculated in solution to compare the oxidizing power of [(L)Ni 3+ -O·] 2+ intermediates in inducing electron transfer from ABTS. Since in general, electron transfer reactions do not involve the formation of transition states, it is not possible using the proposed method to estimate the energy barrier associated with this step of the reaction. However, we may rationalize the electron transfer by treating the reaction system as a Lewis acid-base pair wherein, based on frontier orbital combinations, an electron transfer from HOMO to LUMO is highly permissible when the energies of the involved orbitals are more disparate [79]. The electron is shuttled completely to a more electronegative atom rather than being shared as in the case of adduct formation that involves covalent bonding. Applying the same concept, it is safe to assume that if the oxidant provides an unoccupied orbital that is significantly lower in energy than that of the source, electron transfer between the two redox species would be favorable. In Figure 8, all the intermediates have two low-lying unoccupied orbitals that are nearly degenerate and are both significantly lower in energy than HOMO of ABTS. They can accommodate two electrons, which we think will come either from ABTS or ABTS +· . Apparently, HOMO of the radical form is also high at −4.98 eV considering the energy of the LUMOs in subject. However, the further oxidation of ABTS +· to ABTS 2+ with λmax at 513-520 nm has been reported to be very slow, roughly three times that of the first oxidation [80]. Hence, we considered the parent ABTS as the principal electron donor for this redox process.
Among the six types, 5A/5B and 6A/6B INTMD2s are the two with the lowest LUMO energies each at 1.7 to 1.8 eV lower than the HOMO energy of ABTS. We believe they are the strongest oxidants in the group, since they offer empty orbitals where the two electrons will be the most stable. The low energies of their LUMO and LUMO+1 can be attributed to the very weak antibonding interactions along their Ni-Oaxial bond and Ni-Nequatorial bonds. Their oxidizing power is followed by those of other catalysts in the order: 3A/3B > 1A/1B > 2A/2B > 4.   Nickel-hydroxo species were believed to be the final intermediates of the catalytic cycle before the whole process gets back to the resting state of the catalyst. Based on the generated frontier orbitals, the spin-conversion of INTMD3 from high-spin (S = 1) to low-spin (S = 0) is expected by the time the terminal HO − dissociates from the nickel metal, similar to the reversible coordination of HOO − at the initial step. The ease of the axial HO − to leave the complex is influenced by the nature of the surrounding ligands, which in the case of Ni(II)-PyMACs is determined by the functional group of the pendant arm being the only one unique for each catalyst. It can be assumed that HO − dissociation is slower in complexes containing an electron-withdrawing pendant arm than those that have an electron-donating pendant arm. In the presence of an electron-withdrawing group such as that in 1A/1B, 2A/2B, and 3A/3B, the metal is made more electron-deficient, thereby increasing both its Lewis acid character and its propensity toward nucleophilic binding. However, 1A/1B may also appear as having the greatest tendency to liberate HO − as it can easily stabilize the free ion with its very low pKa. 1A/1B INTMD3 has a pendant arm that can provide ionizable H + to neutralize HO − ion to form H2O. Taking all these into consideration, basic complexes, 2A/2B and 3A/3B, tend to have the lowest turnover number compared to the rest of the catalysts. In this particular step, HO − dissociation from the complex likely follows the trend: 1A/1B > 5A/5B > 4 > 6A/6B > 2A/2B > 3A/3B.

Reaction Energy Profile of Ni(II)-PyMAC Catalysis
The functional state of the pendant arm dictates the reaction pathway of the catalytic action. In the first part of the discussion, we introduced the concept of acid-base equilibria of the pendant arm and how it can affect the reaction rate. Synonymous with experimental findings, the basic form B of each catalyst will boost the reaction kinetics of ABTS oxidation during the initial reaction timeframe. Based on the energy diagram (Figure 9), B will proceed on a two-step, minimal energy process, leading to the formation of INTMD1. B undergoes spin-transition (S = 0→S = 1), and the originally coordinated pendant arm will serve as hydrogen-abstracting moiety to transform H2O2 to HOO − . Such a process is more thermodynamically favorable compared to that of its conjugate counterpart A (e.g., +36.9 kcal mol −1 from 1A to INTMD1). A depends solely on the availability of free basic species to deprotonate the oxidant. However, since A and B coexist in solution, it is likely that apart from H2O and HO − , B may also assist A in generating HOO − prior to INTMD1 formation. Note that due to computational limitation and the complexity of the modeled system, we were unable to locate the first-order saddle point (transition state) associated with this reaction step. We merely relied on the thermodynamic stabilities, frontier orbitals, and spin density distribution of equilibrium structures as a fair start to orchestrate the overall mechanism of the reaction. Remarkably, the results obtained so far were suggestive of the most plausible pathway undertaken by each form of the catalyst. Starting from INTMD1 toward catalyst regeneration in the form of triplet BHS, both forms proceed along the same reaction scheme. In fact, the relative energies that correspond to INTMD1, TS2, and INTMD2 are equal for both functional states. They only deviate when INTMD2 is converted to INTMD3, where it is thermodynamically forbidden for B while thermodynamically favorable for A. This result is surprisingly noteworthy. The high energy requirement for complex B can be attributed to the limited protons that can stabilize the reduced form of INTMD2 after accepting an electron from ABTS. The formation of INTMD3 will only proceed when the system can sustain the proton requirement of the proton-coupled electron transfer without disrupting the equilibrium among participating species. Such a criterion is fulfilled by the reaction environment of A, being more acidic. For complex 1A, the stabilization energy associated with the formation of INTMD3 measures −18.6 kcal mol −1 . While this is not the case for complex B, we still believe that regardless of whether protons are sufficient, electron transfer from ABTS should still irreversibly proceed given the much lower energy of the accepting orbital of INTMD2. However, the major implication of the high-energy formation of INTMD3 in B is the likely low turnover of the complex after going through successful catalysis. B is practically consumed after INTMD2 formation, and little of it will undergo another round of catalysis after serving as an electron repository for ABTS oxidation. Similarly, A may not be regenerated also because it needs to overcome a considerably high energy to convert from BHS (e.g., 22.0 kcal mol −1 from 1BHS to 1A). Such a conversion is thermodynamically restricted, as the liberated HO − from INTMD3 will instead favor the formation of B to maintain equilibrium. Yet, this allows complex A to go another round of reaction catalysis, but this time following the pathway of the basic form, B.
Overall, the calculated reaction energies provide a theoretical explanation as to why in the previous experiment [24], 5B reaction rate was observed to be superior during the first few seconds of the reaction but had gone constant toward the end. Practically, 5B was only able to catalyze up to a certain extent, and that was limited both by the starting concentration of the complex and the actual pH of the solution. In contrast, 5A catalysis was fair throughout as it acted twice in the redox reaction −first as 5A and second as 5B.

Conclusions and Recommendation
The reaction mechanism of ABTS oxidation with H2O2 over Ni(II)-PyMAC catalysts has been investigated using density functional theory. By frontier orbital inspection, the catalysts were grouped in terms of their ability to deprotonate H2O2 and bind to HOO − . Some catalysts would require a free base (e.g., H2O, HO − ) in forming the first intermediate, [(L)Ni 2+ -OOH] + . The calculated energies in solution suggested that the heterolytic cleavage of the O-O bond, which produces [(L)Ni 3+ -O·] 2+ and HO − , is more favorable thermodynamically compared to the homolytic cleavage. The resulting charged species in heterolysis were stabilized significantly and preferentially by the surrounding aqueous solvent. High-spin (S = 1) [(L)Ni 3+ -O·] 2+ was the active intermediate in the catalytic process, oxidizing ABTS in a 1:2 stoichiometric ratio followed by proton transfer. The acid dissociation constant of the pendant arm influenced the liberation of the terminal OH from Ni. The more ionizable the protons, the easier it was for the dissociated hydroxyl ion to be stabilized by neutralization. Lastly, Ni(II)-PyMAC catalysts in their basic form were more difficult to regenerate compared to their conjugate acid counterparts because of the high energy formation of (S = 1) [(L)Ni 2+ -OH] + due to proton requirement. Overall, deprotonated Ni(II)-PyMACs were verified to be

Progress of Reaction
LCOO LCOOH superior catalysts during the first few seconds of the reaction, and their protonated pairs were better toward the end.
The present study will serve as a theoretical framework for future use of Ni(II)-PyMACs as peroxidase-like catalysts. Given the experimental studies available in the literature, prospective computational works may focus on studying the variations in bonding MOs of mononuclear PyMAC complexes with different types of metal center (e.g., Ni, Fe, Co, Cu). Doing so will shed light on the distinctive catalytic strength displayed by each catalyst. In addition, the influence of other solvents on the preferred degradation pathway of the oxidant may have to be studied as well in order to expand the application of these novel catalysts to other reaction systems.