Why Ortho- and Para-Hydroxy Metabolites Can Scavenge Free Radicals That the Parent Atorvastatin Cannot? Important Pharmacologic Insight from Quantum Chemistry

The pharmaceutical success of atorvastatin (ATV), a widely employed drug against the “bad” cholesterol (LDL) and cardiovascular diseases, traces back to its ability to scavenge free radicals. Unfortunately, information on its antioxidant properties is missing or unreliable. Here, we report detailed quantum chemical results for ATV and its ortho- and para-hydroxy metabolites (o-ATV, p-ATV) in the methanolic phase. They comprise global reactivity indices, bond order indices, and spin densities as well as all relevant enthalpies of reaction (bond dissociation BDE, ionization IP and electron attachment EA, proton detachment PDE and proton affinity PA, and electron transfer ETE). With these properties in hand, we can provide the first theoretical explanation of the experimental finding that, due to their free radical scavenging activity, ATV hydroxy metabolites rather than the parent ATV, have substantial inhibitory effect on LDL and the like. Surprisingly (because it is contrary to the most cases currently known), we unambiguously found that HAT (direct hydrogen atom transfer) rather than SPLET (sequential proton loss electron transfer) or SET-PT (stepwise electron transfer proton transfer) is the thermodynamically preferred pathway by which o-ATV and p-ATV in methanolic phase can scavenge DPPH• (1,1-diphenyl-2-picrylhydrazyl) radicals. From a quantum chemical perspective, the ATV’s species investigated are surprising because of the nontrivial correlations between bond dissociation energies, bond lengths, bond order indices and pertaining stretching frequencies, which do not fit the framework of naive chemical intuition.


Introduction
The highly radical scavenging active cholesterol-lowering drug atorvastatin (ATV) [1] is an outstanding success sale story [2]. It was patented in 1985 and approved by the Food and Drug Administration (FDA) in 1996 for medical use. Sold under the name of Lipitor by the world's leading pharmaceutical company Pfizer, it received record high revenues of about 12.8 billion US dollars in 2006, still generated ten billion US dollars in the year of patent loss (2011) and nearly two billion US dollars in 2019. ATV, one of the most prescribed drugs in the US today, is mainly employed to prevent high risk for developing cardiovascular diseases and as treatment for abnormal lipid levels (dyslipidemia). ATV's inhibition of the HMG-CoA (3 hydroxy-3-methylglutaryl coenzyme A) reductase is plausibly related to the high radical scavenging potency against lipoprotein oxidation.
ATV made the object of several theoretical investigations in the past [3,4]. Still, the antioxidant properties of ATV were only recently investigated from the quantum chemical perspective [5]. Unfortunately, as we drew attention recently [6], the only quantum chemical attempt of which we are aware [5] is plagued by severe flaws [6] (e.g., "prediction" of enormous, totally unrealistic O-H bond dissociation energies of ∼ 400 kcal/mol > 17 eV), and this makes mandatory the effort (undertaken in the present paper) of properly reconsidering the antioxidant capacity of ATV and its ortho-and para-hydroxy metabolites in methanol. For the notoriously poor soluble ATV, this solvent is of special interest. ATV is freely soluble in methanol. In addition, antioxidant assays are mostly done in methanolic environment [5,7]. Along with quantities traditionally related to the antioxidant activity, the present study will also reports on the ATV global chemical reactivity indices, relevant bond data as well as spin densities of radical species generated by H-atom abstraction from ATV and related ortho-and para-hydroxylated derivatives (o-ATV, p-ATV, respectively).
Theoretical understanding of the differences between ATV and its ortho-and parahydroxy metabolites, which is missing to date, is of paramount practical importance. A twenty four years old experimental study reported that atorvastatin ortho-and parahydroxy metabolites (o-ATV and p-ATV, respectively) protect, e.g., LDL from oxidation, while the parent ATV does not [8]. Importantly for the results we are going to present in Section 3.5, the free radical scavenging activity of o-ATV and p-ATV was analyzed by the ubiquitous 1,1 diphenyl-2 picryl-hydrazyl (DPPH • ) assay in ref. [8]. Our study is able to provide the first theoretical explanation of this experimental finding.

Computational Details
The results reported below were obtained from quantum chemical calculations wherein all necessary steps (geometry optimizations, frequency calculations, and electronic energies) where conducted at the same DFT level of theory by running GAUSSIAN 16 [9] on the bwHPC platform [10]. In all cases investigated, we convinced ourselves that all frequencies are real. In all calculations we used 6-31+G(d,p) basis sets [11,12] and, unless otherwise specified (see Sections 3.2 and 3.3), the hybrid B3LYP exchange correlation functional [13][14][15][16].
For comparative purposes, we also present results obtained by using the PBE0 [17] functional and Truhlar's M062x [18,19] (see Sections 3.2 and 3.3). Computations for open shell species were carried out using unrestricted spin methods (e.g., UB3LYP and UPBE0). In most radicals, employing the more computationally demanding quadratic convergence SCF methods was unavoidable. We convinced ourselves that spin contamination is not a severe issue. In all these calculations, we invariably found a value S 2 = 0.7501 for the total spin after annihilation of the first spin contaminant, versus the exact value S 2 = 3/4. Still, to better check this aspect, for ATV's cation and anion as well as for the ATV1H and ATV4H radicals (see Section 3.1 for the meaning of these acronyms) we also undertook the rare numerical effort (enormous for molecules with almost 80 atoms) of performing full restricted open shell (ROB3LYP) calculations; that is, not only single point calculations for electronic energy but also geometry optimization and (numerical) vibrational frequency calculations, and all these in solvent. Differences between UB3LYP and ROB3LYP were reasonably small (see Sections 3.2 and 3.3), but they should make it clear that claims (so often formulated in the literature on antioxidation) of chemical accuracy (∼1 kcal/mol) at the B3LYP/6-31+G(d,p) are totally out of place. From experience with much smaller molecules and much simpler chemical structures (e.g., ref. [20]) we had to learn that achieving this accuracy for bond dissociation enthalpies and proton affinity (BDE and PA, quantities entering the discussion that follows) is often illusory even for extremely computationally demanding state-of-the-art compound model chemistries (CBS-QB3, CBS-APNO, G4, W1BD); see, e.g., Figure 10 of ref. [20]. DFT-calculations done by us and by others [21] revealed that, e.g., errors in ionization potential can amount up to 0.7 eV (16 kcal/mol) even when employing the functional B3LYP and the largest Pople basis set 6-311++G(3df,3pd).
Unless otherwise specified, the solvent (methanol) was accounted for within the polarized continuum model (PCM) [22] using the integral equation formalism (IEF) [23]. Although this is the "gold standard" for modeling solvents in the literature on free radical scavenging, one should be aware that this framework ignores specific solvation effects (hydrogen bonds). Because they may play an important role, e.g., in proton transfer reactions, theoretical estimates of PA may not be sufficiently accurate. While this makes comparison with experiment problematic, it should be a less critical issue when comparing among themselves PA values of various antioxidants in a given solvent (e.g., methanol). To better emphasize why we believe that solvent effects in the context of antioxidants deserve a more careful consideration, along with IEFPCM-based results, we also present results obtained in Truhlar's SMD solvation model [24][25][26].
GABEDIT [27] was used to generate molecular geometries and spatial distributions from the GAUSSIAN output (*.log) files. To compute Wiberg bond order indices, we used the package NBO 6.0 [28] interfaced with GAUSSIAN 16. The reason why we use Wiberg bond order indices [29] rather than the heavily advertised Mayer bond order indices [30] was explained elsewhere [31]. All thermodynamic properties were calculated at T = 298.15 K.

Chemical Reactivity Indices
The global chemical reactivity indices investigated in this work are listed below along with their expressions in terms of the ionization potential IP and electroaffinity EA [32][33][34][35][36]: Here, E g ≡ IP − EA is the fundamental (or transport) "HOMO-LUMO" gap [32,37,38]. Noteworthily, the values of IP and EA presented in this paper were calculated as enthalpy differences (cf. Equations (4a) and (6)). Estimating IP and EA using the eigenvalues of the Kohn-Sham (KS) orbitals with reversed sign (Koopmans theorem), is unfortunately a very popular approximation, but it is totally inadequate especially in the presence of a solvent. For clarification, a comment is in order at this point. Although both the approach using Equation (2) and the approach using Equations (4a) and (6) are based on the DFT, there is an important difference between them.
Equations (4a) and (6) rely on total energies computed via DFT. In these computations, the Kohn-Sham (KS) orbitals merely enter as eigenfunctions and eigenvalues of a mathematical (minimization) problem. They are auxiliary mathematical objects useful to compute a quantity with a clear physical meaning (namely, the total electronic energy).
However, it should be well known to any well rounded theoretician that the KS "orbitals" do not have any physical meaning; they are not real molecular orbitals [32,38,39]. What makes Equation (2) problematic is just the fact that it treats the KS-HOMO and KS-LUMO as if they were the true HOMO and LUMO of a real molecule.
To remedy the difficulty related to the KS "energies" (in reality, eigenvalues of a mathematical single-particle problem) in semiconductor physics, which translates into KS-band gaps typically amounting to about 50% of the real band gap, a so-called "scissor" operator procedure is applied [40,41], which consists of empirically shifting the KS eigenvalues. To eliminate this severe difficulty in the case of molecules immersed in solvents, we also proposed a scissor technique [42]. The important difference is that the scissor corrections proposed in ref. [42] are obtained from quantum chemical calculations rather than empirically as done in semiconductor band structure calculations.
Switching back, one may expect that the global chemical reactivity indices can give a flavor of the overall stability of a molecule and are useful in predicting how a certain chemical environment evolves in time [43,44]. In certain situations they turned out to be useful for comparing properties of different molecular species [33,45,46].
The presently calculated global chemical reactivity indices of ATV and its metabolites are collected in Tables 1 and 2, and depicted in Figure 6. Having a chemical hardness η of about 2 eV, ATV, o-ATV, and p-ATV exhibit a good chemical stability. This value lies between the values of the natural antioxidants phenol and trolox, for which our calculations at the same B3LYP/6-31+G(d,p)/IEFPCM level yielded η = 2.56 eV and η = 1.88 eV, respectively. For all three species, the electrophilic indices [33,45,46] are ω ≈ 1.8 eV, a value exceeding the value of 1.50 eV, which is considered the threshold for strong electrophiles [47]. For comparison, let us again mention the values ω = 1.61 eV and ω = 1.85 eV computed by us for phenol and trolox, respectively.
Inspection of Table 1 reveals that, similarly to the quantities η and ω considered above, all global chemical reactivity of ATV, o-ATV, and p-ATV are comparable to those of well known natural antioxidants. Could we then expect that ATV flavors (or other molecules) have indeed good antioxidant potency merely based on global chemical reactivity indices comparable to those of good antioxidants? Table 1. Global chemical reactivity indices (eV) computed via B3LYP/6-31+G(d,p)/IEFPCM for atorvastatin and its ortho-and para-hydroxy metabolites and two natural oxidants in methanol. The analysis in the next section will unravel that, in fact, the global chemical indices have little relevance for assessing the antioxidant activity of a certain molecule. For the time being, let us remark that the values of Table 1 would rather suggest that ATV and o-ATV have similar antioxidant properties, and that ATV (possibly) performs (slightly) better than p-ATV.

Antioxidant Mechanisms and Pertaining Enthalpies of Reaction
As is widely discussed in the literature, an H-atom can be transferred to a free radical in one or two step processes. The three antioxidative mechanisms (HAT, SET-PT, and SPLET) and the corresponding reaction enthalpies (BDE, IP and PDE, PA and ETE, respectively) can be expressed as follows: Direct hydrogen atom transfer (HAT) [48][49][50] Stepwise electron transfer proton transfer (SET-PT) [51,52] Sequential proton loss electron transfer (SPLET) [53,54] In our specific case, X stands for an O or an N atom. Related to the above (albeit not directly entering the aforementioned antioxidation mechanisms), the electron attachment process is quantified by the electroaffinity defined as BDE, IP, PDE, PA, and ETE are enthalpies of reaction which can be obtained as adiabatic properties from standard ∆-DFT prescriptions [38,55,56]. To this aim, along with the enthalpies of the various ATV-based species entering the above reactions, the enthalpies of the H-atom, proton and electron in methanol are also needed [6]. They are presented in Table 3. Table 3. Gas phase enthalpies H 0 and solvation enthalpies ∆H sol in hartree utilized in the present calculations. For the for the gas phase enthalpy of the H-atom we used the value for the B3LYP/6-31+G(d,p) electronic energy (−0.500273 hartree) and the value of 1.4816 kcal/mol for thermal correction to enthalpy common for all compound model chemistries from GAUSSIAN 16.
The presently computed thermodynamic parameters are collected in Tables 4 and 5, and depicted in Figure 7. Table 4. The enthalpies of reaction (in kcal/mol) needed to quantify the antioxidant activity of atorvastatin (ATV) and its ortho-and para-hydroxy metabolites (o-ATV, p-ATV).  The really important effect brought about by the extra OH-group of the hydroxy metabolites is the homolytic bond dissociation at its position (5-OH). Our calculations demonstrate that this process is substantially less expensive energetically than H-atom donation from position 1-OH. The calculated BDE values for both o-ATV and p-ATV at this position are ∼77.5 kcal/mol versus the smallest value ∼91 kcal/mol for ATV at position 1-OH, respectively. Unlike the extremely similar homolytic bond dissociation, there is a certain difference between o-ATV's and p-ATV's heterolytic bond dissociation at position 5-OH, as expressed by the PA values (PA = 34.4 kcal/mol = PA = 37.9 kcal/mol, respectively). However, it is unlikely that this difference in PA's has practical consequences, again because the aforementioned values of PA are both comfortably larger than the lowest PA = 23.8 kcal/mol at position 1-OH, a value that also characterizes the parent ATV molecule.

Molecule
In Section 3.5 we will return to the practical implications of the above finding.

Alternative Approaches to the O-H and N-H Bond Strengths: Vibrational Frequencies and Bond Order Indices
Let us start this section with a short digression. The robustness of a single molecule diode fabricated using the scanning transmission microscopy (STM) break-junction technique [61,62] can be quantified by the maximum force that the junction subject to mechanical stretching can withstand. This rupture (pull-off) force F per molecule, which characterizes the strength of the chemical bond between electrodes and the terminal (anchoring) atom of the embedded molecule, can hardly be directly measured. To circumvent this difficulty, experimentalists use a simple mechanical model which relates F to the vibrational frequency of the pertaining stretching mode ν. The latter quantity can be easily measured by infrared spectroscopy [63]. To exemplify, this is the Au -S stretching mode in benchmark nanojunctions wherein molecules are anchored via thiol groups on gold electrodes.
Transposed to the present context, it is interesting to interrogate the relationship between BDE and the related stretching frequency. In the same vein, a stronger chemical X-Y bond is intuitively expected to have not only a larger BDE and a higher stretching frequency ν(X-Y) but also a shorter length and a larger bond order index.
With these in mind, let us examine the correlation of the aforementioned quantities in the presently considered molecules.
Infrared spectra calculated for ATV, o-OH-AVT, and p-ATV in methanol are depicted in Figure 8.  Table 4) should be done within a multi-reference framework.
The behavior visible in Figure 8b is surprising for several reasons, e.g.,   Let us now switch to bond order indices. Our results are collected in Table 6 and Figures 9 and 10.
To reiterate, based on straightforward chemical intuition, it would be obvious to expect that stronger chemical bonds (larger BDE's) possess larger bond order indices. Figure 9b depicts that for the O-H bonds of ATV, o-ATV, and p-ATV just the opposite holds true: larger BDE's justly correspond to smaller bond order indices. As for their N-H bonds, Figure 9b reveals that the dependence is even nonmonotonic.
To avoid misunderstanding, a clarification is in order before ending this analysis. What chemical intuition in the above example should not overlook is that a pair of atoms X and Y forming an X-Y chemical bond, do not merely interact with each other but also with the neighboring atoms in the molecular surrounding. This is also why a simple (exponential [64]) relationship between bond order indices and bond lengths can hold, e.g., for homologous molecular series [65], but cannot not hold in general; otherwise one arrives at comparing apples with oranges. Figures 9 and 10 illustrate this again using the values of Table 6. BDE values corresponding to different O-H bonds of a given molecule differ from each other depending on the specific chemical environment. These differences can be visualized by inspecting the spin density landscape of the various radicals (Figures 1, 3-5). The stronger the delocalization in a radical, the easier is its formation, and the lower is the corresponding BDE value. Inspection of Figure 1b,c makes it clear, e.g., why ATV's BDE at position 3-OH is higher than that at position 1-OH.

Assessing the Radical Scavenging Activity-A Specific Example
Discussion on free radical scavenging and dominant antioxidant mechanism is very often couched by comparing among themselves values the enthalpies characterizing the HAT, SET-PL, and SPLET of the specific antioxidant(s) under investigation. Every now and then publications conclude, e.g., that SPLET is the dominant pathway because a certain antioxidant has a "small" PA value or a PA substantially smaller than BDE, or that SET-PL prevails because of the small IP value. However, it is worth emphasizing that, along with the antioxidant's properties, a proper evaluation of the antioxidant activity should mandatory consider the specific properties of the radicals to be eliminated (neutralized).
The small value BDE ≈ 77.5 kcal/mol for o-ATV and p-ATV, substantially smaller than the smallest value (BDE = 90.2 kcal/mol) of the parent ATV, is perhaps the most appealing result reported in Section 3.3. Still, the "small" value mentioned above does not demonstrate per se the fact anticipated in Introduction, namely that o-ATV and p-ATV can scavenge can scavenge the ubiquitously employed 1,1-diphenyl-2-picrylhydrazyl (DPPH • ) radical, while the parent ATV cannot.
To demonstrate this, one should mandatory consider the pertaining DPPH • property, namely the enthalpy release in DPPH • 's neutralization (H-atom affinity) Because it amounts to 80 kcal/mol [66], e.g., the reaction is exothermic. H-atom abstraction from position 5-OH of o-ATV (or p-ATV) costs ∼77.5 kcal/mol, a value lower that the enthalpy release of 80 kcal/mol [66] in the neutralization of the DPPH • radical, and this makes the HAT mechanism thermodynamically allowed. Rephrasing, because the BDE of the N -H bond of DPPHH is 80 kcal/mol [66], o-ATV (and p-ATV) can scavenge the DPPH • radical through donating the H-atom at position 5-OH. On the contrary, the parent ATV cannot. The lowest ATV's BDE (at position 1-OH) amounts to 90.2 kcal/mol (Table 4), so the HAT pathway is forbidden.
To conclude, we have presented above the first theoretical explanation of the experimental fact [8] that the antioxidant properties of atorvastatin ortho-and para-hydroxy metabolites differ from those of ATV.
By and large, there is a consensus in the literature that HAT is a possible (or even preferred) antioxidant mechanism in the gases phase but not in polar protic solvents like the presently considered methanol. In this vein, the natural question that arises is: can o-ATV and p-ATV scavenge the DPPH • radical in methanol also via SPLET? Can HAT and SPLET coexist? While the large IP (Table 4) give little chances to an SET-PT pathway, SPLET would a priori be conceivable in view of the "small" value of PA, which is, although not smaller than that of ascorbic acid (as incorrectly [6] claimed in ref. [5]) at least not much larger than the latter (23.8 kcal/mol for ATV's versus 20.5 kcal/mol for ascorbic acid, see ref. [6]).
In fact, Table 4 implicitly gives the negative answer to this question. If o-ATV and p-ATV could scavenge DPPH • via SPLET, then (contrary to experiment [8]) the parent ATV could also do the job; the most favored deprotonation, implying the same enthalpy PA = 23.8 kcal/mol, occurs both for ATV and its metabolites at the same 1-OH position, where furthermore the similar spin density landscapes (compare Figure 1b with Figure 3) indicate a similar chemical reactivity.
Still, let us remain in the realm of theory and demonstrate why neither o-ATV nor p-ATV or ATV can scavenge DPPH • in methanol via SPLET. To this aim suffice it to consider the first step of SPLET where x means "o-", "p-", or "nothing". Straightforward manipulation allows to express the enthalpy of this reaction as follows . (10) Notice that the second brace in Equation (10) corresponds to the proton abstraction from the cation DPPHH •+ of the neutralized free radical DPPHH, or alternatively, the PDE pertaining to the neutralized free radical DPPHH (cf. Equation (4b)).
Our calculations yielded PDE(DPPHH) = 3.9 kcal/mol, a value that is not larger (as the case if the first SPLET step was allowed) but smaller than PA(xATV) = 23.8 kcal/mol. It now becomes clear why neither ATV, nor o-ATV or p-ATV can scavenge the DPPH • radical via SPLET. Their "small" PA is not small enough to fulfill Equation (11).

Conclusions
We believe that the present demonstration that atorvastatin ortho-and para-hydroxy metabolites can scavenge the DPPH • through donating the H-atom at the position of their extra group (5-OH), which is impossible in the parent ATV, is important not only because it theoretically explains for the first time a behavior revealed in experiment [8] but also because, from a general perspective, it provides further insight into the structure-activity relationship (SAR).
By working out a specific example (Section 3.5)-an analysis that can be straightforwardly extended to other cases-, we drew attention that an adequate approach to antioxidant's potency should mandatory account for the thermodynamic properties of the free radicals. Equation (11) expresses a general necessary condition for thermodynamically allowed SPLET, and its application to specific cases may reveal that, even in polar solvents, free radical scavenging via this pathway is forbidden not only for ATV-based species.
In addition, our study emphasize that, while important, e.g., for modeling the temporal evolution of various molecular species interacting among themselves in a given chemical environment [65,67], the global chemical reactivity indices have no direct relevance for antioxidation. Recall that we saw in Section 3.2 that quantitative differences of ATV's o-ATV's, and p-ATV's global chemical reactivity indices are minor. Furthermore, if qualitative differences in these indices were important, then, contrary to Sections 3.3 and 3.5, o-ATV would have antioxidant properties similar to ATV rather than to p-ATV.
Last but not least, from the perspective of fundamental science, we found (Section 3.4) that properties like bond dissociation enthalpy, bond order index, bond length, and bond stretching frequency, expected after all to represent alternatives in quantifying the bond strength, are by no means correlated according to naive intuition. This finding calls for further quantum chemical efforts aiming at comprehensively characterizing ATV's, that inherently remained beyond the scope of this study focused on ATV's antioxidant activity. Finally, the presently reported counter-intutitve relationship between bond stretching frequency and bond strength should also be a word of caution for other communities; for example, for the molecular electronics community, wherein bond stretching frequencies (conveniently obtained via infrared spectroscopy) are used to estimate (pull-off) forces that cause the rupture of a junction subject to mechanical stretching [68].

Data Availability Statement:
The data that support the findings of this study are available from the author upon reasonable request.

Acknowledgments:
The author is much indebted to Ederley Vélez Ortiz for providing valuable details related to her recent work [5].

Conflicts of Interest:
No conflict of interest to declare. Table A1. Z−matrix of ATV.