Density Functional Theory Transition-State Modeling for the Prediction of Ames Mutagenicity in 1,4 Michael Acceptors

Assessing the safety of new chemicals, without introducing the need for animal testing, is a task of great importance. The Ames test, a widely used bioassay to assess mutagenicity, can be an expensive, wasteful process with animal-derived reagents. Existing in silico methods for the prediction of Ames test results are traditionally based on chemical category formation and can lead to false positive predictions. Category formation also neglects the intrinsic chemistry associated with DNA reactivity. Activation energies and HOMO/LUMO energies for thirty 1,4 Michael acceptors were calculated using a model nucleobase and were further used to predict the Ames test result of these compounds. The proposed model builds upon existing work and examines the fundamental toxicant−target interactions using density functional theory transition-state modeling. The results show that Michael acceptors with activation energies <20.7 kcal/mol and LUMO energies < −1.85 eV are likely to act as direct mutagens upon exposure to DNA.


■ INTRODUCTION
Global frameworks such as green chemistry should sit at the forefront of chemical research in the modern era. To help achieve this ideology, Anastas and Warner proposed twelve principles of green chemistry that acts as fundamental guidelines for the development of greener chemical processes and products. 1 The field of computational toxicology is concerned with using qualitative and quantitative data to develop models that improve our understanding about the toxicological risk of chemicals. 2,3 This directly aligns with the fourth principle of green chemistry; to design safer chemicals and reduce toxicity wherever possible. Furthermore, computational toxicology plays a pivotal role in attempting to reduce the widespread use of animal testingan ethical issue of great concern to many people around the world. 4 Despite continuous improvements in the treatment of cancer, carcinogenicity continues to be a toxicological endpoint of great concern when developing new chemicals. Mutagenesis, the process by which genetic mutations are created, is the underlying process in the development of many cancers. Understanding the chemical mechanism of mutagenic activity is therefore vital to producing less toxic chemicals. 5 The most widely used bioassay for testing mutagenic activity is the Ames test. This test uses Salmonella typhimurium bacteria with pre-existing mutations that prevent the synthesis of histidine. Without histidine, the bacteria cannot grow. When a chemical shows mutagenic activity, DNA reactivity is apparent and this can cause a mutation that allows the bacteria to revert to a state where it can produce histidine and therefore start to grow. 6 The extent of bacterial growth therefore gives quantitative and qualitative indication to the mutagenic potential of a test chemical. The prediction of Ames test results using computational toxicology is possible, however, many currently used methods often rely on expert knowledge of chemical systems. 7 These in silico prediction systems, such as Derek Nexus, group chemicals into categories (e.g., DNA reactive, see Figure 1) based around so-called "structural alerts". 8,9 Examples of structural alerts within the category of DNA reactivity are α,β-unsaturated carbonyls, acyl halides, and isocyanates to name but a few. However, category formation based upon the presence of a functional group (or "structural alert") does not provide an adequate explanation of the fundamental chemistry associated with mutagenicity. The initial reaction between an exogenous chemical and DNA can be described as a "molecular initiating event". 10 This reaction can trigger the "Adverse Outcome Pathway", which is defined as a series of events between initial exposure and the adverse outcome, for example, cancer, respiratory failure, and so forth. 9 Acquiring a thorough understanding of the fundamental chemistry behind the molecular-initiating event would allow us to better predict the outcome of experimental methods in toxicology such as the Ames test. This could then lead to a high throughput in silico screening process for drug candidates, making predictions both faster and cheaper. Quantum chemistry provides the necessary toolset to probe the molecular-initiating event. In particular, density functional theory (DFT) transitionstate modeling has previously been shown to allow a thorough assessment of the steric and energetic details of chemical reaction mechanisms. 11 This is typically achieved through calculating activation energies to assess which reaction pathway is likely to proceed. Goodman et al. have previously shown that DFT-derived activation barriers can be used to predict Ames test results for a variety of Michael acceptors which form covalent bonds with DNA. 12 This work improves upon the published model by Goodman et al. and demonstrates that the energy of the lowest unoccupied molecular orbital (LUMO) shows significant predictivity for the assessment of mutagenic potential in Michael acceptor-type compounds.

■ MATERIALS AND METHODS
Ames test results for 30 archetypical α,β-unsaturated carbonyl compounds were obtained from two resources: the OECD QSAR Toolbox and a dataset published by Peŕez-Garrido and co-workers. 13,14 19 of these compounds were obtained from the study by Goodman and co-workers. 12 At the time of this study, 16 of the total compounds were experimentally classified as Ames negative whilst 14 of the compounds were classified as Ames positive. All experimental Ames test data were acquired without the use of an S9 enzyme activation system. DFT was used to study the reaction between the Michael acceptor dataset and a model nucleophile (see Figure 2). Thus, the activation energies were derived for the molecular initiating events. To ensure a reasonable calculation time, methylamine was chosen as the nucleophile to model a nitrogenous DNA nucleobase (see Figure 3). Compounds 16 and 23 were truncated from citral and pulegone, respectively. 15 Where conformational flexibility was evident, conformational searches were performed using Schrodinger's Macromodel (ver. 11.3) with the MMFF force field. 16−18 A low-mode sampling approach was used for both individual Michael acceptors and transition states. 19 Reactant and transition-state conformations obtained from Macromodel were optimized using DFT calculations within Gaussian 16 (Revision A.03). 20 Optimizations were performed using Grimme's D3 dispersion-corrected B3LYP functional with Becke−Johnson damping and the polarized triple-zeta valence (def2-tzvpp) basis set to minimize basis set superposition errors. 21,22 Throughout this work, water was used as a solvent. The implicit integral equation formalism− polarizable continuum model was used, which previously, has been extensively utilized for modeling organic chemical reactions. 23,24 Temperature and concentration-corrected quasiharmonic free energies (under the Grimme approximation 25 ) were obtained using GoodVibes. 26 This included a vibrational scaling factor of 1, a temperature of 310.15 K, and a concentration of 1 mol dm −3 . 26 Electronic Supporting Information was generated for the lowest energy ground and transition state conformer for each compound using ESIGen. 27 ■ RESULTS AND DISCUSSION Activation Energies. In total, the activation energies for 30 compounds were obtained (see Table 1). The results showed a wide range of activation energies, ranging from 11.4 to 33.5 kcal/ mol (see Figure 4). Ames positive compounds ranged between 11.4 and 20.7 kcal/mol. The Ames negative compounds ranged between 22.0 and 33.5 kcal/mol with two outliers at 13.7 and 22.8 kcal/mol (see discussion below). The average free energy of activation for Ames positive compounds was 16.8 kcal/mol, whilst the average for Ames negative compounds was 24.9 kcal/ mol. This leads to a significant energy gap of 8.1 kcal/mol between the average Ames positive and Ames negative compounds.
For the 19 compounds published in the dataset by Goodman et al., the energy gap separating the highest Ames positive compound and the lowest Ames negative compound was improved from 2.2 to 2.9 kcal/mol for these 19 compounds. This was achieved by using a dispersion-corrected B3LYP functional and a larger basis set. 12,28 Upon extension of the model, the Ames positive compound with the highest activation energy of 22.8 kcal/mol was compound 9, and an energy barrier of this magnitude sits within the expected region for Ames negative compounds. However, Haddon et al. previously showed evidence that compound 9 is indeed Ames positive with the TA100 strain of S. typhimurium. Despite this classification, compound 9 was shown to exhibit an extremely low molar mutagenicity of 0.24 rev/nmol, explaining why its activation energy sits within the negative region. Further to this, in the same study, an Ames negative, nonmutagenic compound was quoted to have a similar molar mutagenicity of <0.3 rev/ nmol. Therefore, compound 9 was removed from the model, and the authors recommend complete re-evaluation of Ames test activity for compound 9. A further outlier was compound 1, initially an Ames negative compound obtained from the OECD QSAR toolbox with an activation energy of 13.7 kcal/mol. Upon inspection of the compound activation energy and LUMO energy, it became evident that compound 1 shows prototypical values expected from an Ames positive compound. To explore this result, the existing literature for compound 1 was examined.
A study was found that shows positive activity in the Ames test, thereby validating the prediction of our model. 29 The Ames test result for this compound was changed accordingly.
The chosen compounds in this study included a carboxylic acid, compound 30. The standard acidic strength of carboxylic acids sits between 2 and 5 pKa, and at physiological pH, the predominant form of compound 30 would therefore be the ionized carboxylic acid species. 30 For this reason, calculations were performed on the conjugate base of compound 30. The experimental Ames test result for this compound was negative,     HOMO/LUMO Energies. The highest unoccupied molecular orbital (HOMO)/LUMO energies for 30 compounds were calculated. Despite previous successful attempts in literature, the energy of the HOMO showed no predictive potential for mutagenicity as a toxicological descriptor (see Table 1). 31,32 However, the LUMO energy of Michael acceptors proves to be an extremely promising descriptor. The LUMO energies ranged from −3.65 eV (for the lowest lying LUMO) to −0.81 eV for the highest calculated LUMO energy (see Table 1). The LUMO energies for Ames positive compounds ranged from −3.65 to −1.85 eV, with an average Ames positive LUMO energy of −2.57 eV. The Ames negative compounds ranged from −1.83 to −0.81 eV, with an average Ames negative LUMO of −1.50 eV. It is evident that Ames positive compounds have lower LUMO energies when compared to Ames negative compounds. Our results show that as a toxicological descriptor, LUMO energy shows great promise for accurately predicting Ames test results in Michael acceptors. Further to this, LUMO energies are easily and routinely obtained from quantum chemical calculations. 33 Furthermore, the LUMO energies showed statistically significant correlation with the activation energies being calculated (see Electronic Supporting Information). Regression analysis was performed resulting in an r 2 = 0.75, suggesting good correlation between activation energy and LUMO energy. Also, the coefficient of correlation "r" was equal to 0.89, suggesting that as the reaction barrier increases, there is a simultaneous increase in LUMO energy.

Journal of Chemical Information and Modeling
Comparison to the Existing Model. In order to assess the performance of our model, a pre-existing (Q)SAR software was used to predict the in vitro Ames test result for the compound dataset. The "Toxicity Estimation Software Tool" (TEST) was used to predict Ames test results through a "nearest neighbor" approach. 34 In this approach, the mutagenicity is predicted by comparing the test chemical to the three most similar compounds in the QSAR training set. TEST correctly predicted 9 of the 14 Ames positive compounds and 11 of the 15 Ames negative compounds. Despite TEST showing relatively good performance, the model, we present, predicts Ames test results with greater consistency. It also provides a deeper insight into the fundamental chemistry associated with toxicant−target interactions. Overall, both LUMO energies and activation energies show great promise for use as primary (Q)SAR descriptors and readily meet guidelines within toxicological policy. For example, the International Council for Harmonisation declares that two complimentary in silico methods should be used for the assessment of the mutagenicity of impurities in pharmaceutical products. 35

■ CONCLUSIONS
Activation free energies and HOMO/LUMO energies were calculated for 30 Michael acceptors. For Ames positive compounds, the activation energies ranged from 11.4 to 20.7 kcal/mol, whilst the Ames negative compounds ranged from 22.0 to 33.5 kcal/mol. LUMO energies ranged from −3.65 to −1.85 eV for Ames positive compounds and from −1.83 to −0.81 eV for Ames negative compounds. Compound 1 was originally classified as Ames negative; however, it was corrected to Ames positive upon application of our model and assessment of the literature. We can predict with confidence that Michael acceptors with a LUMO energy < −1.85 eV and an activation energy <20.7 kcal/mol will be Ames positive, whilst those with activation energies >22.0 kcal/mol and LUMO energies > −1.83 eV will be Ames negative. These results have shown that DFT transition state modeling shows great potential for understanding the molecular initiating event and therefore the intrinsic chemical origin of mutagenicity. When activation energy and LUMO energy are combined as descriptors, they show excellent predictive potential for Ames test classification in Michael acceptor-type compounds. By considering the information presented above, the authors recommend that the activation energy and LUMO energy are thoroughly and simultaneously investigated when probing the mutagenic potential of pre-existing and newly synthesized 1,4 Michael acceptors. In future, this framework can be applied to different groups of compounds that are known to react with DNA, for example, Schiff base formers. 36 However, the range of activation barriers and LUMO energies is likely to differ between separate groups of compounds. Therefore, DFT transition-state modeling should be applied independently to each category. In conclusion, we have developed a fast, low-cost framework, with advantages over current methods (e.g. TEST) for assessing the mutagenic risk of pharmaceutically interesting Michael acceptors. Activation energy and LUMO energy show significant predictivity for mutagenic risk in Michael acceptors, and these descriptors should be incorporated into both pre-existing and future QSAR models in predictive toxicology.

Author Contributions
The manuscript was written through contributions of all the authors. All the authors have given approval to the final version of the manuscript.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
This work was supported by the Engineering and Physical Sciences Research Council (EP/L016354/1) and the University of Bath. Part of this work was completed using the Balena HPC service at the University of Bath (https://www.bath.ac.uk/ corporate-information/balena-hpc-cluster/). Dr. Marianne Ellis is thanked for helpful discussions regarding this work.