Coupled molecular design diagrams to guide safer chemical design with reduced likelihood of perturbing the NRF2-ARE antioxidant pathway and inducing cytotoxicity†

The NRF2-ARE antioxidant pathway is an important biological sensing and regulating system that responds to xenochemicals. NRF2 senses chemically-caused production of reactive oxygen species (ROS) and electrophilic interactions with chemical species. Upon NRF2 activation, the expression of a wide array of genes will be upregulated to counteract oxidative or electrophilic insults. However, when the external disruption exceeds the inherent resilience of the biological system, cellular damage can occur, eventually leading to cytotoxicity. Induced NRF2 activity in in vitro assays is therefore a signal that a man-made chemical may cause unwanted cellular activity. This was the motivation to derive a chemical design strategy to minimize the risk that new chemicals would perturb this pathway. We constructed a logistic regression model using design variables derived from density functional theory (DFT) calculations and physical properties. The model showed excellent predictive power to distinguish between NRF2-active and inactive chemicals based on the EPA ToxCast high throughput screen (HTS) assay data (tested in the concentration range of 10–10 μM). External evaluation showed that the area under the curve (AUC) for the receiver operating characteristic (ROC) of the model is 0.81 and the precision is 0.90. Combining this model with a previously developed cytotoxicity model, we developed a probabilistic design diagram to guide chemical design with the twin goals of minimizing NRF2 antioxidant pathway activity and cytotoxicity. This work initiated a simultaneous design strategy against two toxicity pathways of interest in molecular design research.


Introduction
The fourth principle of green chemistry states that "chemical products should be designed to preserve their desired function while minimizing their toxicity". 1 This idea leads to the desire among chemists to seek design solutions to reduce the likelihood of chemicals perturbing various toxicity pathways. In this report, we chose the NRF2-ARE antioxidant pathway as a target for study. NRF2 (Nuclear factor E2-Related Factor 2) is a bZip transcription factor and a member of Cap 'n' Collar family of regulatory proteins. 2,3 At the basal level, NRF2 is repressed by Kelch-like ECH-associated protein 1 (Keap1) and constantly undergoes ubiquitin-assisted degradation. 4,5 Upon stimulation initiated by redox active or electrophilic chemicals, signaled by the presence of reactive oxygen species (ROS), the Keap1-NRF2 interaction is destabilized or the degradation of NRF2 is hindered. In either case, free NRF2 forms a transcription factor that causes upregulation of anti-oxidation genes via the antioxidant response element (ARE). [6][7][8] The resulting detoxifying enzymes are upregulated to both reduce ROS levels and clear xenobiotics. 9,10 At low levels of perturbation (low xenobiotic concentration, short exposure durations), NRF2 activation can effectively protect the cell from the effects of ROS or electrophilic attacks. This is considered to be a protective mechanism to help an organism adaptively respond to this type of insult. The NRF2-ARE pathway was thus proposed as a potential therapeutic target for anti-oxidation intervention. [10][11][12][13] However, at greater levels of exposure, the protection mechanism can be overwhelmed and adverse consequences may follow. Oxidative stress provides a well-known example of tipping over the inherent redox balance in the biological system due to the pressure from ROS. NRF2-Keap1 is responsible for sensing this event. 5,[14][15][16] The adverse effects of oxidative stress have been linked to multiple clinical syndromes including neurodegenerative diseases, cancer, arthritis, atherosclerosis and inflammatory diseases. 3,[17][18][19][20][21] Therefore, man-made chemicals should be designed to minimize the likelihood of generating ROS. The NRF2-Keap1 is used as a biosensor in this study to gauge the potential of a chemical inducing ROS generation and posing an electrophilic threat. It should be noted that certain chemicals may bypass the NRF2-Keap1 monitor but still be involved in ROS generation. Those cases deserve separate treatment and are not considered this research endeavor.
Three critical components are needed to build a predictive model of NRF2 activity to be used in green molecular design. First, we need a reliable biological assay that probes the chemically induced activation of the NRF2-ARE pathway. In vitro high throughput screening (HTS) methods answer this need. HTS is a cost-effective technology to examine how chemicals disrupt biological pathways and hence potentially lead to adverse health outcomes. It embodies a paradigm shift in toxicology from using in vivo to using in vitro and in silico testing approaches, as described by the U.S. National Research Council (NRC) in their report on Toxicity Testing in the 21st Century. 22,23 In order to evaluate practical approaches to implement the "Tox21" vision, the U.S. National Toxicological Program (NTP), the U.S. Environmental Protection Agency (EPA) and the NIH National Center for Advancing Translational Sciences (NCATS) collaboratively forged a collaborative research partnership. This Tox21 partnership is using HTS methods to test thousands of chemicals in a wide variety of cells, pathways and technologies, relevant to many aspects of chemical toxicity. [24][25][26][27][28][29] The NRF2-ARE antioxidant pathway has been probed using the ToxCast/Tox21 assay battery. Second, we need to construct a design variable space with physically meaningful and accurately calculable parameters. Physical observables and chemically intuitive variables are preferred because of their accessibility from experimental measurements. Additionally, these design variables should also encapsulate mechanistic information regarding the molecular initiating events that lead to the toxicity endpoint of interest. Thirdly, the predictive model needs to be built so that the design variable space can be made searchable and the inverse function of the model can be presented in a probabilistic and complete manner. This requires a careful design of the modeling algorithm.
Previously, acute aquatic toxicity, 30,31 mutagenicity/carcinogenicity 32 and cytotoxity 33 have been investigated as single target toxicity endpoints to derive green molecular design strategies. In this report, we present for the first time in molecular design research a coupled design diagram that simultaneously considers two design goals, in this case the minimization of both cytotoxity and the NRF2-ARE pathway perturbation. This attempt is motivated by the significance that certain chemicals can induce one toxic response without the other. With the progress of molecular design research, it is anticipated that more and more models can and will be linked together to provide a comprehensive design strategy.

Toxicity data resource
The U.S. EPA Toxicity ForeCaster (ToxCast™) program release 34 contains two assays that were designed to monitor chemicals perturbing the NRF2-ARE toxicity pathway: ATG_NRF2_ARE_CIS_up (subsequently abbreviated as ATG_NRF2) and Tox21_ARE_BLA_agonist_ratio (subsequently abbreviated as Tox21_ARE).
ATG_NRF2 is part of a multiplexed assay using an inducible reporter that measures mRNA for gain-of-signal activity at the transcription factor-level as they relate to the gene NFE2L2. This is one of 52 assay component measured or calculated from the ATG_CIS assay. Activity is detected with fluorescence intensity signals by reverse transcription polymerase chain reaction (RT-PCR) and capillary electrophoresis technology. Activity is quantified by the level of mRNA reporter sequence unique to the cis-acting reporter gene response element ARE, which is responsive to the endogenous human nuclear factor, erythroid 2-like 2 [GeneSymbol:NFE2L2|GeneID:4780| Uniprot_SwissProt_Accession:Q16236]. The assay is cell-based and uses HepG2, a human liver cell line, with measurements taken at 24 hours after chemical dosing in a 24-well plate. 35 Tox21_ARE uses an inducible reporter to measure gainof-signal activity related to the gene NFE2L2. The assay measures beta-lactamase induction, as detected with fluorescence intensity signals by GAL4 beta-lactamase reporter gene technology. Changes to fluorescence intensity signals produced from an enzymatic reaction involving the key substrate [CCF4-AM] are indicative of changes in transcriptional gene expression due to agonist activity regulated by NRF2. The assay is cell-based, producing a single-readout that uses HepG2, with measurements taken at 24 hours after chemical dosing in a 1536-well plate. 29 The two assays were tested on chemicals with similar concentration ranges, approximately 10 −3 -10 2 μM (see Fig. S1 in ESI †). Consequently, our predictive model is strictly applied to the likelihood for chemicals incurring biological responses over this particular concentration range. The testing outcome concordance between these two assays is 0.73 (see Fig. S2 in ESI †). Cytotoxicity was also examined in the ToxCast data. Assay data related to cytotoxicity were used to construct a cytotoxicity model that was linked to the NRF2-ARE toxicity pathway model. Detailed description of the cytotoxicity model is provided in our recent publication. 33

Chemical selection
Four criteria were set for chemical selection in the following analyses: 1. Molecular weight < 1000. The bioavailability for chemicals with molecular weight greater than 1000 is negligible. 36 2. Single compound with a definite structure, excluding geometrical and optical isomers and mixtures. The molecular mechanisms for mixtures can be complex, and are not considered here.
3. Containing no metal elements. The calculation of electronic properties of metals is significantly more challenging than for organics, so the focus here is on the latter. 4. Concordant testing outcomes between ATG_NRF2 and Tox21_ARE assays (both positive or both negative). Concordance is required because both assays have false positives and false negatives, so chemicals with concordant behavior are more likely to be correctly assigned. 37 A total of 1261 chemicals satisfied these four criteria. Binary hit calls provided in the ToxCast release 34 were used to categorize chemicals into two classes: chemicals detected to be not perturbing the NRF2-ARE toxicity pathway (inactive) and otherwise (active). The chemical membership distribution is biased towards the inactive chemicals, which account for ∼72% of the total sample population.
All chemical data were evenly split into training and test sets with the same ratio between active and inactive chemicals for cross-validation and external evaluation. The first treatment of the chemical structures is to remove counter-ions from the molecules using the open-source OpenBabel chemistry toolbox. 38 Subsequently, 3D structures with the lowest energy conformation were generated using the ChemAxon Marvin calculator plugin. 39

Generation and selection of design variables
Seven design variables were generated and examined in this study: molecular softness (SOF), electrophilicity index (EPH), ionization potential (IP), electron affinity in the aqueous phase (EA.aq), polarizability (PLRZ), radical stability energy (RSE) and log P (log water-octanol partition coefficient). log P was calculated using the ChemAxon Marvin calculator plugin. 39 The remaining variables were computed based on DFT implemented in Gaussian 09 rev.D.01. 40 Boese and Martin's τ-dependent hybrid functional 41 and basis set 6-31+G(d) were used for full geometry optimization. Vertical IP and electron affinity (EA) were calculated in vacuum. The species with an extra electron were then transferred into the implicit aqueous environment based on a universal solvation model 42 to obtain vertical EA.aq. SOF and EPH were derived according the following formulas. 43 Radical stability energy (RSE) is defined as the reaction energy for the reaction: where MH is a generic presentation of a chemical with H being a hydrogen atom and M • is the leftover radical after one hydrogen atom is removed from MH.
The numerical distribution of each design variable is plotted in ESI Fig. S3. †

Mechanistic rationales for design variables
The molecular initiating event for NRF2 activation is the chemical modification of thiol groups on Keap1. 7,8 Such modification can be effected through a covalent interaction or a redox mechanism. SOF is a quantitative variable that describes the likelihood for covalent interaction. 44 Covalent interactions account for a number of toxicological events including the activation of NRF2-ARE pathway and cytotoxicity. 5,45 Therefore, it is reasonably included in the predictive model design. For the redox event, EA.aq describes the energy gain of a chemical after acquiring one electron. Here, it can viewed as a measure of the potential of a chemical stripping one electron from critical biological molecules in the aqueous phase. The higher the EA.aq value, the higher is the reactivity of a chemical. It was within the original design variable consideration, but discarded later because of its inferior predictive power compared to other preserved variables (see Table S1 in ESI †). In order for xenobiotics to interact with Keap1 inside a cell, they must permeate the membrane. log P has been frequently used to estimate such likelihood and has shown acceptable predictivity in quantitative structure-activity relationship (QSARs) of acute aquatic toxicity. 46 Subsequent to molecular transport into a cell, xenobiotics proceed to interact with biological targets. This process is driven at least partially by induced electric dipole moments and dispersion forces. 47 Polarizability (PLRZ) is a physical quantity that describes the relative tendency towards molecular electron cloud distortion under the influence of an external electric field. This reflects the energy alteration during a molecular recognition process. This provides the physical basis for including PLRZ in the predictive model.
The last variable to mention is radical stability energy (RSE). It approximates the degree to which a chemical is liable to form reactive radicals, which is the trigger for releasing the NRF2-Keap1 interaction. However, the low predictive power of this variable in the model suggests that direct radical activation is not the dominant mechanism of NRF2 activation in chemicals included in this study. It is consistent with the fact that HepG2 cell lines do not possess full metabolic activities [48][49][50] and that biotransformation is an important avenue for generating reactive radicals. 51,52

Predictive model construction and design probability derivation
The predictive model for this work is built on the logistic regression algorithm. 53 where π(X) denotes the probability for a chemical not perturbing the NRF2-ARE toxicity pathway (benign probability) and X denotes the design variables. The predictive power of the seven design variables were assessed using ROC AUC analysis (see ESI Table S1 †). The three most predictive variables (PLRZ, SOF and log P) were selected based on LASSO regularization 54 for the predictive model construction. The logistic regression model was trained on ∼600 chemicals with 10-fold cross validation, i.e. chemicals were randomly partitioned into 10 equal sized groups and a single group is retained as validation set for the model built upon the other nine groups in rotation until every group has been tested. Then the established model was externally tested on the remaining ∼600 chemicals. Model construction, data analysis, graphics and design diagrams in this study were coded using the Python programming language and the R statistical computing environment. 55-67

Performance evaluation for the predictive model
We used the following measures to evaluate the performance of the logistic model: Because we are interested in producing 'safer' chemicals in the design process, we define chemicals that do not perturb the NRF2-ARE antioxidant pathway as positive in this study. Therefore, type I error is the focus of interest from the green design perspective. Precision measures the probability that a chemical will be a true positive relative to the fraction of chemicals that are predicted to be positive. This probability is of high interest to green chemists. It should be maximized for the purpose of safer chemical design. Accuracy is a measure of the balanced performance of the model. Given that the overall rate for a chemical to be inactive against the NRF2-ARE pathway is imbalanced, the ROC AUC (equivalent to the balanced accuracy) is a more robust measure for the overall classification accuracy of the predictive model. The model performance is summarized in Table 1.
From Table 1, one can see that three performance indicators are in reasonable agreement. Additionally, the crossvalidation and external evaluation results agree well with one other. This suggests that the classifier has reasonable forward predictive performance and can be used for deriving green design probabilities for novel chemicals.

NRF2-ARE pathway and cytotoxicity
Phase II enzymatic detoxification is triggered by electrophilic and redox-active species, and among other actions serves to maintain a stable intracellular environment. 68,69 When chemical perturbations exceed the cellular recovery capacity, molecular damage can result that could eventually lead to cell death. Therefore, a correlation is expected to exist between NRF2-ARE pathway perturbation and cytotoxicity. For the purpose of comparison, the NRF2-ARE and cytotoxicity assay results 33 are summarized in Table 2. There is an 82% concordance between the two assays, among which 30% of chemicals cause positive responses in both assays. 99 chemicals (∼14%) specifically invoke cytotoxicity without perturbing the NRF2-ARE pathway and 23 chemicals (∼3%) solely activate the NRF2-ARE system.
Mechanistically, the NRF2-ARE pathway activation can be triggered by either covalent or redox modifications of the thiol groups on Keap1. 15,16 Cytotoxicity can be induced by multiple causes including cell membrane disruption, activation of specific signaling pathways or DNA modification. 70,71 Given the multitude of mechanisms of cytotoxicity, it is of interest to investigate the differences between chemicals that induce only one or both of NRF2-ARE activity and cytotoxicity. The following statistical discussion will focus on the chemicals that only cause cytotoxicity (cytotox group, 99 chemicals) or incur both cytotoxicity and NRF2-ARE response (cytotox + ARE group, 209 chemicals). Table 3 shows the calculated ROC AUC 72,73 and p-values from a Mann-Whitney U test 74 for the three design variables (SOF, PLRZ and log P) used in the logistic model. In accordance with the Mann-Whitney U test, SOF strongly suggests that the two groups of chemicals (cytotox and cytotox + ARE group) are drawn from two unequivalent distributions. This is consistent with the ROC AUC result, where SOF exhibits the highest discriminative power among the three design variables. The overlaid SOF histogram (Fig. 1) further illustrates that the distributions of the two chemical groups are stag-   gered, with the cytotox group slightly shifted to the left. Based on the analysis, it follows that the two groups of chemicals have differential selectivity toward invoking the NRF2-ARE pathway and cytotoxicity. The distributions are not totally separated, which is an indication that SOF is not solely responsible for separating the two groups. Many processes leading to toxicity are driven by chemical reactivity toward critical biological targets including DNA, enzymes, lipids and others, and thus should be able to be described by chemical design variables. In the current case, SOF was found to be helpful in discriminating between the cytotox and cytotox + ARE chemical groups. This can be understood from Pearson's HSAB theory 75 -acid-base reactions preferentially occur between chemicals with a similar level of softness/hardness. Since covalent bond formation plays a crucial role in chemical initiation of cytotoxicity and the NRF2-ARE pathways, it is thus plausible to expect certain discriminative power from SOF in distinguishing the two classes of interest. In a biochemical interaction, active functional groups on proteins and DNA normally act as nucleophiles, while xenobiotics often act as electrophiles. 76,77 The sulfhydryl groups on protein cysteine residues are generally softer than the amino and hydroxy groups on nucleic acids. 78 In rapidly dividing cells, cytotoxicity has been attributed to chemical modification of DNA. 79,80 Interactions between chemicals and critical proteins can also lead to cytotoxicity. 45 The thiol group is an important functional unit on many proteins, and is responsible for sensing a variety of the chemical signals. The NRF2 activation mechanism is an example of the thiol group functioning as a molecular switch. 81 Therefore, it is reasonable that chemicals with a preference to induce only cytotoxicity may initiate this event primarily through reacting with DNA, which is a 'harderharder' interaction. On the other hand, chemicals that cause both cytotoxicity and NRF2-ARE activity are more likely to achieve this effect through modifying thiol groups, 81 which is Fig. 1 Histogram overlay for SOF. Color coding: yellow-chemicals causing both cytotoxicity and perturbing the NRF2-ARE pathway, green -chemicals only causing cytotoxicity, vertical blue bar-median value for chemicals causing cytotoxicity only, vertical red bar-median value for chemicals perturbing the NRF-ARE pathway and incurring cytotoxicty. Fig. 2 Coupled nomograms for NRF2-ARE activity and cytotoxicity. 33 P(non-NRF2-ARE): probability of not activating NRF2-ARE ( purple). P(non-cytotox): probability of not inducing cytotoxicity (blue). The benign probability axis is colored while other design variables axes are printed black. R1 indicates an auxiliary axis. The numerical range of each design variable axis are calculated using 5 to 95 percentiles of numerical data to produce probabilities corresponding to the range of the benign probability axis. Benign probability above 50% is considered for the green molecular design purpose. a 'softer-softer' interaction. Thus, it is observed that the SOF distribution for the cytotox chemical group is shifted to the left (harder) compared to the cytotox + ARE group.

Probabilistic design diagram
The full design diagram is composed of two parts (Fig. 2). The left panel depicts the probability of chemicals not activating NRF2-ARE and the solution to any probability in the design variable space. The right panel is the design diagram for reducing the likelihood of causing cytotoxicity. 33 Given the high concordance between cytotoxicity and NRF2-ARE activation, one would expect to observe a degree of similarity between the two diagrams. However, its should be kept in mind that the concordance matrix ( Table 2) possesses non-zero off-diagonal elements. This means that the predicted probabilities may vary between the two toxicity endpoints. To obtain a design solution with a target probability of not stimulating the NRF2-ARE pathway and not causing cytotoxicity, one can simply start from desired probabilities for each endpoint and work out values for design variables on the diagrams. A sample solution is given in Fig. 2, using dotted lines. In this example, suppose we have planned probabilities for a chemical not perturbing the NRF2-ARE pathway (73%) and the probability of not incurring cytotoxicity (58%). We start by identifying the positions of these two probabilities on their corresponding probability axes and link them with a dotted line. Next, we put dots on the desired values for softness (0.12) and log P (2.9) on the design variable axes because those values fit into the feasible design space. Then we connect those dots on the design variable axes with the dots on the probability axes and we arrive at the values for polarizability. This case study shows the different probabilities for non-cytotoxic and non-NRF2-ARE activity for the same set of values of the design variables. Of course, users have the freedom to choose any three fixed points on either diagram to determine the last unknown parameter (either a design variable or a probability). One can also work on one side of the diagram first to determine the benign probability and then to check its probability of not invoking the other endpoint by filling those numbers on the other half of the diagram.

Conclusions
We have trained a classifier based on the logistic regression algorithm to identify chemicals that do not disrupt the NRF2-ARE antioxidant pathway in the concentration range of 10 −3 -10 2 μM. The overall performance of this model is 0.81 for ROC AUC and 0.90 for precision, as assessed using an external validation set. Based on this model, we derived a probabilistic diagram to guide molecular design to reduce the chance for chemicals perturbing the NRF2-ARE pathway. In addition, we statistically demonstrated that molecular softness (SOF) possess information in distinguishing the distributions between chemicals that cause cytotoxicity, with or without activating NRF2. This effect was mechanistically rationalized as the preferential interaction between chemicals and either DNA or proteins. Further, we coupled the NRF2-ARE activity and cytotoxicity to arrive at a linked diagram to provide more comprehensive information for chemists in designing safer chemicals. This coupled diagram is the first such toxicity-oriented multi-objective initiative in the molecular design literature. We foresee that linking multiple toxicity endpoints is a general trend in the future of safer chemical design research.

EPA disclaimer
The views expressed in this paper are those of the authors and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency.