A Simple Computational Tool for Accurate, Quantitative Prediction of One–Electron Reduction Potentials of Hypoxia–Activated Tirapazamine Analogues

- The reduction potentials of bioreductively-activated drugs represent an important design parameter to be accommodated in the course of creating lead compounds and improving the efficacy of older generation drugs. Reduction potentials are traditionally reported as single–electron reduction potentials, E (1), measured against reference electrodes under strictly defined experimental conditions. More recently, computational chemists have described redox properties in terms of a molecule’s h ighest o ccupied m olecular o rbital (HOMO) and l owest u noccupied m olecular o rbital (LUMO), in electron volts (eV). The relative accessibility of HOMO/LUMO data through calculation using today’s computer infrastructure and simplified algorithms make the calculated value (LUMO) attractive in comparison to the accepted but rigorous experimental determination of E (1). This paper describes the correlations of eV (LUMO) to E (1) for three series of bioreductively–activated benzotriazine di- N -oxides (BTDOs), ring-substituted BTDOs, ring-added BTDOs and a selection of aromatic nitro compounds. The current computational approach is a closed–shell calculation with a single optimization. Gas phase geometry optimization was followed by a single–point DFT (Density Functional Theory) energy calculation in the gas phase or in the presence of polar solvent. The resulting DFT–derived LUMO energies (eV) calculated for BTDO analogues in gas phase and in presence of polar solvent (water) exhibited very strong linear correlations with high computational efficiency ( r 2 = 0.9925) and a very high predictive ability (MAD = 7 mV and RMSD = 9 mV) when compared to reported experimentally determined single–electron reduction potentials.


INTRODUCTION
Focal hypoxia is caused by microenvironmental deficiency of oxygen as a direct result of the rapid growth and aberrant vasculatures in most solid tumors. As a malicious metabolic aberration, tumor hypoxia enhances the potential for cancer metastasis, treatment failure and hypoxia-based radio-and chemotherapy resistance as compared to oxygenated tumors. Taken together, hypoxia-induced molecular aberrations are major obstacles for effective tumor management [1][2][3][4][5]. The main basis of hypoxic tumor radiosensitization therapy has been the utilization of electron-affinic organic compounds that are bioreductively-activated selectively in hypoxic microenvironments to induce radiation-like molecular damage similar to what happens in presence of oxygen [7][8][9][10]. Of these compounds, tirapazamine (TPZ; Fig.1) is a benzotriazine di-Noxide (BTDO) that has been subjected to extensive preclinical and clinical studies as a radiosensitizer [10]. TPZ requires a single-electron bioreductive activation [11][12][13] to produce free radicals which in turn induce single and double strand breaks in DNA, exploiting the hypoxic microenvironment to be selectively cytotoxic to solid hypoxic tumors ( Fig. 1) [13][14][15][16][17]. _________________________________________ Despite the initial clinical success of TPZ in conjunction with radiation therapy for the management of hypoxic tumor, the poor extravascular transport (EVT) of TPZ to the tumor site was a major limitation. The first-electron reduction potential, E(1), and the partition coefficient (LogP) are the two physico-chemical properties that directly influence the rate of hypoxic metabolism and EVT, respectively, of TPZ and TPZ-analogues, processes that are necessary to achieve effective chemical concentrations in vivo for radiosensitization therapy [18][19][20].
To overcome the poor EVT of TPZ, several TPZ-analogues have been synthesized in attempts to optimize the balance between E(1) and LogP that would improve hypoxia-induced radiosensitization and antiproliferative effects [21][22][23][24][25][26][27][28]. Recently, there is an increasing interest in exploiting the hypoxiaselective properties of TPZ to develop molecular theranostics of hypoxia [26,27,29] and to repurpose TPZ-analogues for treatment of serious anerobic bacterial infection such as E. coli, S. aureus, C. difficile [30] and M. tuberculosis [31]. In drug discovery, the TPZ-analogues are typically synthesized first and then their experimentally-derived E(1) is measured, a protocol that is timeconsuming, expensive and labor-intensive. Ideally, an optimized computational method that can accurately predict the reduction potentials of the bioreductive drugs, in the current context the proposed TPZ-analogues, would guide medicinal chemists in the process of decision making by selecting only those compounds that show favourable theoretically-predicted reduction potentials. Based on their E(1), potentially active molecules can be selected for full synthesis development and complete biological evaluation, which may substantially enhance an otherwise ineffective drug development approach.
A good linear correlation between the Hückelcalculated LUMO (lowest unoccupied molecular orbital) energies and the experimentally-derived reduction potentials of hydrocarbons was first demonstrated by Maccoll in 1949 [32]. This finding prompted extensive research in the field of computational chemistry, with the purpose of finding better computational methods that include both the solvent stabilisation and molecular reorganizational energies in the calculation of LUMO energies, which Maccoll`s method excluded. Interested researchers are referred to selected papers for more detailed information [33][34][35][36][37][38][39]. It is noteworthy to mention that both Fry et al [37] and Gillmore et al [38] have provided good computational methods for accurately predicting the redox potentials of mainly polycyclic aromatic hydrocarbons (PAHs), expanded to few other classes in case of Gillmore`s method [39], with good correlation to the experimentally-derived reduction potentials. Despite both methods having been proven to provide accurate values, they are computationally expensive and require an experienced computational chemist, as they rely on three [38] to six [37] calculations per molecule. The high computational cost and complexity of both methods were resolved by Mujica et al [40]. These authors reported a simple computational method that accurately correlates DFT (Density Functional Theory)-calculated HOMO (highest occupied molecular orbital) /LUMO energies (HLE) to the experimentally-derived redox potentials through a single point calculation per molecule using B3LYP/6-31G(d) functional. This method provided good correlation between the DFT-calculated HOMO/LUMO energies and the redox potentials of 51 PAHs, as well as to other organic compounds of different structural families [41]. Despite the accuracy and efficiency of this method, a linear correlation must be established and calibrated before any prediction can be made. Therefore, Mujica et al suggested that new correlations should be established for other structural families of compounds, functional groups or solvents to maintain accurate predictability, especially because most compounds used in their methodology were rigid organic compounds [40].
Our special interest is to design new TPZanalogues as theranostic tools for management of hypoxia and treatment of anaerobic bacterial infections. The aim of this work is to exploit a quantitative tool to predict the reduction potentials of newly proposed molecules that is accurate, fast, cost-efficient and simple enough to be conducted in medicinal chemistry laboratories. With this purpose in mind, the implementation, calibration and application of a linear correlation between the DFTcalculated LUMO energies and experimental first electron reduction potentials using the method reported by Mujica et al are now reported for a series of TPZ and nitroimidazole analogues.

Computational Method
All calculations were done on an HP Envy Phoenix 810-160 with 4 th generation Intel® Core i7-4770 processor, 16GB of DDR3 system memory and 16 GB SATA 6G SSD acceleration cache. Spartan 16 Parallel Suit (Wavefunction Inc., Irvine CA) was used for optimization of the calibrant molecules [42,43]. Molecules 1-64 were constructed using the 3D model kit of the program and preliminarily minimized by the program's comprehensive molecular mechanics tool. These structures were subsequently submitted to geometry optimization in gas phase and followed by single-point energy calculation in gas phase and in polar solvent (water). The Conductor-like Polarizable Continuum Model (CPCM) [44,45] was employed in calculations involving polar solvent (water). All calculations were carried out at DF level of theory with the B3LYP [46,47] hybrid density functional and 6-31G(d) [48][49][50] basis set. Geometry optimization and single-point energy calculations were performed in a single step using the intuitive graphical user interface (GUI) of Spartan 16 Parallel Suit.

RESULTS AND DISCUSSION
The experimental E(1) values for tirapazamine analogues (1 -64, Tables 1 -4) were obtained from the literature [21][22][23][24]. The E(1) experiments are usually conducted in anaerobic aqueous solution at pH 7 by measuring the one-electron transfer equilibrium constants between the radical anions of the test compounds and a reference standard (viologen or quinone), with compensation for the effect of ionic strength on the equilibrium constant by collecting the data at three concentration ratios at room temperature (22 ± 2 °C) and using them to calculate the ∆E between the test compounds and the appropriate reference [24,51,52].
The DFT-calculated LUMO energies of compounds 1-64, calculated for the gas phase and for the presence of polar solvent (water), were plotted against their experimental E(1) values to construct linear correlations, Figure 2. A summary of the parameters of the linear correlations between the computed LUMO energies and the experimental E(1) values established in the current study is provided in Table 5.      Four linear correlations between the DFTcomputed LUMO energies and the experimental E(1) were established in the current study. Correlation #1, between the DFT-calculated LUMO energies of all calibrant compounds (1-64) and their experimental E(1) values in gas phase, resulted in low r 2 (0.737), indicative of low computational efficiency. Correlation #3 included the solvent in the calculation, which improved the r 2 (0.951). This variance is due to the absence of the solvent stabilization effect on compounds 1-64 in gas phase calculations, which span a wide range of structural substitutions in the tirapazamine parent structure that in turn induce various levels of structural flexibility and therefore different solvent reorganization energies [39]. It is noteworthy that Correlation #2 for compounds 1-32 yielded an higher r 2 (0.915) for the gas phase calculation, an effect primarily due to the similar solvent reorganization energies of compounds 1-32, as they differ only in the substitution pattern of the aromatic ring of the tirapazamine parent structure. The substitution pattern imparts only a limited effect on the structural flexibility of these molecules, and hence the error associated with not including the solvent in the calculation is very small. Furthermore, when the solvent was included in the calculation for compounds 1-32 (correlation #4), the r 2 was greatly improved, approaching the optimal value of 1 (0.9925), demonstrating high computational efficiency and high predictive ability as reflected in the very small values of MAD and RMSD. To further test and establish the predictive ability of the current method, Correlation #4, which demonstrated the strongest correlation in the current study, was used to predict the reduction potentials of five test compounds, 65-69 [53] that were not included in the data set used to establish these correlations ( Figure  3). Furthermore, test compounds are not structurally related to TPZ; 65-67 are 5-nitroimidazoles, 68 is a 2-nitrofuran and 69 is a 2-nitroimidazole. The experimental one-electron reduction potentials E(1) of these test compounds were determined at pH 7 in aqueous solution using pulse radiolysis [7]. The predicted E(1) values calculated from Correlation #4 for the test compounds were very close to the experimental values, providing a further demonstration of the robust predictive ability of this correlation (Table 6).

CONCLUSION
These analytical data show that a very strong linear correlation between the DFT-calculated LUMO energies and the experimental E(1)s of TPZ analogues can be established using a polar solvent model. The established correlation has a high computational efficiency (r 2 = 0.9925) and a very high predictive ability (MAD = 7 mV and RMSD = 9 mV), which can be extended to other hypoxiaselective compounds, e.g., 65-69, that are structurally unrelated to TPZ.
The current methodology used a closed-shell calculation and a single optimization.
Gas phase geometry optimization followed by a single-point energy calculation for the presence of polar solvent reduced both the number and the complexity of the computations.
In conclusion, this approach provides a simple, quick and accurate quantitative method to predict reduction potential without requiring the skills, computing resources and expense normally associated with molecular computation. By removing the necessity for a strong background in computational chemistry and the knowledge of complex computational protocols, our objective is to encourage medicinal chemists to exploit computational tools through easy-to-use graphical user interfaces. This can accelerate research in drug discovery and optimization, and as applied in this case, lead to better management of hypoxia.

ACKNOWLEDGEMENTS
This work was supported through grants 201201164 and 20131051 from Alberta Innovates through the CRIO program and CRIO Project-Cancer.