Predictive Chemometric Modeling and Three-Dimensional

Article pubs.acs.org/IECR

Predictive Chemometric Modeling and Three-Dimensional Toxicophore Mapping of Diverse Organic Chemicals Causing Bioluminescent Repression of the Bacterium Genus Pseudomonas Supratik Kar and Kunal Roy*,† Drug Theoretics and Cheminformatics Laboratory, Department of Pharmaceutical Technology, Jadavpur University, Kolkata 700032, India S Supporting Information *

ABSTRACT: Classification and regression-based quantitative structure−activity relationship (QSAR) as well as threedimensional (3D) toxicophore models were developed for toxicity prediction of 104 organic chemicals causing bioluminescent repression of the bacterium genus Pseudomonas isolated from industrial wastewater. Statistically significant and interpretable in silico models were obtained using linear discriminant analysis (classification), genetic partial least-squares (regression), and 3D toxicophore models. The QSAR and toxicophore models were scrupulously validated internally as well as externally along with the randomization test to avoid the possibilities of chance correlation. Features such as octanol−water partition coefficient, thirdorder branching, CH2 fragment or unsaturation, and the presence of a higher number of electronegative atoms (specifically halogen atoms) and their contribution toward hydrophobicity have been identified as major responsible structural attributes for higher toxicity from the developed in silico models. The present approaches can provide rich information in the context of virtual screening of relevant chemical libraries for aquatic toxicity prediction.

■

INTRODUCTION Industrial and domestic wastewaters are majorly treated using the activated sludge process (ASP).1 High concentration of organic matter and salts and the presence of poorly biodegradable organic compounds and their large variation in volume as well as compositions make the process a difficult one.2 This leads to accumulation of large organic toxicants in municipal wastewater3 and a potential threat to the potable drinking water.4 Operation of the ASP is strongly inhibited due to the presence of large number of toxicants in the ASP influent. Problems such as a decrease in waste organics removal, reduction of solids separation efficiency, and modification of sludge compacting properties are also added ones.5 These problems can be avoided if ASP influent wastewater is tested for toxicity and protecting actions are taken as early as possible. Testing of wastewater toxicity can involve a number of chemical and biological procedures. Bacteria-based assays are generally used for influent wastewater toxicity monitoring, and these assays are majorly considered simple and time- and costeffective.6 Moreover, biological treatment can influence the overall micropollutant removal.7 Various microbial toxicity tests such as growth inhibition, enzymatic activity (dehydrogenase), ATP contents, bacterial luminescence (e.g., Microtox), metabolic heat production (e.g., microcalorimetric techniques), and respiration rate have been developed to determine toxicity toward activated sludge.8−10 Among all these methods, the bioluminescent bacteria-based microtox assay is one of the important and mostly studied assays, but the response to toxicants of Photobacterium phosphoreum is different from that of the activated sludge microorganisms. Ren and Frymier11−13 developed the Shk1 assay which is an improved influent wastewater monitoring system based on the bioluminescent bacterium Shk1, a genetically modified Pseudomonas isolated © 2013 American Chemical Society

from the activated sludge of industrial wastewater treatment plants. As the Shk1 assay is less sensitive than the microtox assay, it could therefore be more suitable for influent wastewater toxicity monitoring. With the rapid encroachment of the world’s population accompanied by industrial progress, there is a significant impact of chemicals including pharmaceuticals on the world’s ecosystems.14 Since human society depends on so many types and classes of chemicals, it is very important to know the way in which the chemicals interact with the environment as well as mankind. Every year a huge number of chemicals are submitted to Environmental Protection Agency (EPA) for toxicity studies. In these particular consequences, applications of bacterial-based assays are not only time-consuming but also costly. Again, increasing pressure from social and economic backgrounds to eliminate the use of animal testing is an important reason to develop alternative methods,15 such as in silico QSAR models which also support the 3Rs (replacement, refinement and reduction in animals in research)16 and Registration, Evaluation, Authorisation and Restriction of Chemicals (REACH) policies.17 Globally, it is difficult for the chemical industry and regulatory agencies to make decisions regarding exposure guidelines for environmental contaminants when such experimental data are not available regarding systemic toxicity.18 In such circumstances, in silico approaches, specifically quantitative structure−activity relationship (QSARs) and three-dimensional (3D) toxicophore models have the ability to predict such potential health hazards from Received: Revised: Accepted: Published: 17648

August 26, 2013 November 14, 2013 November 18, 2013 November 18, 2013 dx.doi.org/10.1021/ie402803h | Ind. Eng. Chem. Res. 2013, 52, 17648−17657

Industrial & Engineering Chemistry Research

Article

toxicity data is presented in Figure S1 in Supporting Information. Descriptor Calculation. A pool of 218 descriptors was calculated using Cerius 2 version 4.10,27 Dragon 6,28 and PaDEL-Descriptor version 2.1129 software. Using Descriptor+ module of the Cerius 2 software, 232 descriptors belonging to various categories were thus calculated for the present work which included: (i) spatial, (ii) topological, (iii) thermodynamic, (iv) electronic and (v) structural and (vi) E-state parameters. After those descriptors having variance less than 0.0001 were excluded, a total pool of 147 descriptors was chosen from Cerius 2. A set of 49 descriptors (constitutional and functional group counts) was selected after considering 99% intercorrelation cutoff among 197 descriptors obtained from Dragon software. Additionally, 22 extended topochemical atom indices (ETA) (both the first and second generations) were calculated using the PaDEL-Descriptor software. Along with the theoretically calculated 218 descriptors, experimental octanol−water partition coefficient (log Kow) values reported by Ren and Frymier11−13 have also been utilized as one predictive descriptor. Data Set Splitting. Selection of the training and test set plays an important role in the development of a statistically significant QSAR model. The selection should be such that the test set molecules should lie within the chemical space occupied by the training set molecules. In this study, the entire data set was randomly divided between training and test sets keeping the compounds from both classes (higher and lower toxic groups). In the case of the descriptor-based QSAR, 30% of the compounds were selected as the test set (ntest = 31), while the remaining 70% (ntraining = 73) were used as the training set. On the contrary, in the case of 3D toxicophore models, the test set of QSAR models was used as the training set, and the training set of the QSAR model was used as the validation or test set. The major reasons behind this step are (i) utilizing both sets for model development to extract all possible features from the total data set and (ii) generally small number of compounds are used for model development in 3D toxicophore generation. The training-set molecules were used to develop the in silico models, while the predictive abilities of the models were assessed using the test set. Model Development. In this present work, we have attempted to build three different types of models: (a) a classification-based QSAR model to identify the major discriminatory features between higher and lower toxicants, (b) regression-based QSAR models to quantify the contributions of the structural attributes to toxicity, and (c) 3D toxicophore models to identify the features responsible for the toxicity. The descriptor-based QSAR models were built using four different chemometric tools: (a) stepwise multiple linear regression (stepwise MLR),30 (b) genetic function approximation (GFA),31 (c) genetic partial least-squares (G/PLS),32 and (d) using linear discriminant analysis (LDA)33 techniques. The 3D toxicophore model was built using the conformers obtained from the BEST method of conformer generation where conformational analysis of the molecules was done using the poling algorithm.34 The 3D toxicophore was generated using the HypoGen module implemented in Discovery Studio 2.1 software.35 A minimum of 0 to a maximum of 5 features including hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), hydrophobic aliphatic (HYD-Aliphatic), hydrophobic aromatic (HYD-Aromatic), and ring aromatic (RA) features were selected for the 3D toxicophore generation.

chemical exposure. The QSARs not only save time but also provide valuable resources which could be endowed more sensibly.19,20 In recent times, QSAR has been applied to predict the toxicity of different chemicals, specifically toward the environment as well as humans. Xu and Nirmalakhandan21 have applied QSAR models to predict joint toxicity of mixtures of organic chemicals to microorganisms. Schultz et al. developed QSAR models to predict aquatic toxicity of a set of halosubstituted carbonyl compounds.22 Katritzky et al. developed a four parameter QSAR model to predict organic compounds causing bioluminescent repression of the bacterium genus Pseudomonas isolated from industrial wastewater.23 Pharmaceutical toxicants are also modeled by different groups of researchers. Christen et al. (2010) evaluated the mode-of-action concept, the fish plasma model, and a QSAR model to predict the effects of pharmaceuticals in the aquatic system.24 Kar and Roy have developed robust quantitative interspecies toxicity correlation models for Daphnia magna and fish assessing the ecotoxicological hazard potential of structurally diverse pharmaceuticals.25 With this perspective, we have attempted in the present report to develop classification and regression-based QSAR models along with 3D toxicophore models for 104 organic chemicals causing bioluminescent repression of the bacterium genus Pseudomonas isolated from industrial wastewater. All the developed models have been constructed in consonance with the spirit of the Organization for Economic Cooperation and Development (OECD) guidelines.26 The constructed in silico models offer an important tool toward predicting the aquatic toxicity of large number of chemicals. The aim of this work has been to develop integrated in silico models for preliminary identification of higher toxic and lower toxic classes of compounds with the classification-based model and then to quantify the prime structural attributes which are majorly responsible for aquatic toxicity of the studied toxicants by regression-based QSAR and 3D toxicophore models. The developed models can provide rich information in the context of virtual screening of relevant chemical libraries regarding aquatic toxicity.

■

MATERIALS AND METHODS Data Set. The experimental toxicity values of 104 organic chemicals causing bioluminescent repression of the bacterium genus Pseudomonas isolated from an industrial wastewater were taken from literature sources11−13 to build up classification and regression-based QSAR and 3D toxicophore models. Potency Threshold. Except for the 3D toxicophore, half maximal effective concentration (EC50) of the molecules was converted to the negative logarithmic scale (pEC50) in nanomolar (nM) units for the development of the QSAR models. On the contrary, EC50 values of these molecules were expressed in nM to develop the 3D toxicophore models. The mean of the pEC50 toxicity values (obtained from the normality distribution plot) was 2.8. Compounds having a pEC50 value of 2.8 or more were classified as higher toxic agents or positive (H), and those with pEC50 values less than 2.8 were termed as lower toxic agents or negative (L) substances. For the toxicophore model, compounds with an EC50 value of approximately 1500 nM (≈2.8 in logarithmic scale) or above were categorized into the lower toxicity group, and compounds with EC50 values between 0 to 1500 nM were termed as the higher toxicity group. The normality distribution plot for the 17649

dx.doi.org/10.1021/ie402803h | Ind. Eng. Chem. Res. 2013, 52, 17648−17657


Article

A value of 2 was ascertained to the uncertainty parameter.36 The 10 hypotheses thus generated were analyzed in terms of their correlation coefficients and the cost function values. Software. Software tools such as Discovery Studio 2.1,35 STATISTICA 7.0,37 SPSS 9.0,38 MINITAB 14,39 and SIMCAP 10.040 were used in the present study for developing in silico models. Validation Metrics. Statistical metrics were employed to check the fitness of the in silico models, and different internal, external, and overall validation strategies were subsequently employed for model validation. Metrics for Classification-Based QSAR Models. The common classification evaluation metrics, namely, accuracy, sensitivity, specificity, precision, and F-measure were calculated along with the more advanced methods of geometric means (G-means) and Cohen’s kappa that are not affected by imbalanced data sets.41,42 Classification capability of the developed model was judged by a number of statistical tests. Wilks’ λ statistics,43 Canonical index (Rc),44 Matthews correlation coefficient (MCC), squared Mahalanobis distance43 were computed, and a receiver operating characteristic (ROC) curve45 was plotted to check the model quality. The developed classification model was also judged for its quality using the chisquare (χ2) statistic that detects the independence between two groups or classes signifying that higher value of this parameter will indicate greater separability between groups. The ROC curve obtained by plotting the sensitivity and (1-specificity) indices along the Y and X axes respectively identifies the discrimination ability of the classification system. The performance of a diagnostic variable can be quantified by calculating the area under the ROC curve (AUROC). The ideal test would have an AUROC of 1, whereas a random guess would have an AUROC of 0.5. Additionally, we have calculated two parameters,46 namely, the ROC graph Euclidean distance (ROCED) and the ROC graph Euclidean distance corrected with Fitness Function (FIT(λ)) (ROCFIT) to have better explainable results. In the present study, we have used another validation diagram named pharmacological distribution diagram (PDD)47 to verify the degree and extent of discrimination achieved in the training as well as test set observations. Metrics for Regression-Based QSAR Models. The goodness-of-fit of the equations was judged by the quality metric R2, as well as the internal validation metric leave-one-out cross-validation parameter QLOO2 and external validation metric Rpred2 or Qext(F1)2. The rm2 metrics,48 namely, rm 2 and Δrm2 developed by the present authors’ group for internal, external, and overall validation of models are also employed for the present work. The calculation of the rm2 metrics for the test set data ( rm 2(test), Δrm2(test)) additionally estimated the closeness

Q2ext(F2)54 and Golbraikh and Tropsha’s55 criteria to check the model reliability. Metrics for 3D Toxicophore Model. The toxicophore model was validated using two methods: Fischer’s validation for the training set, and test set prediction using the external validation method. In order to judge the external predictive quality of the toxicophore model, classification-based statistical parameters (sensitivity, specificity, accuracy, precision and Fmeasure) were computed. Y-Randomization Test. The final regression-based QSAR model was subjected to a randomization test. Additionally, in the case of the toxicophore model, the Fischer randomization technique at 95% confidence level was employed to judge the possibility of chance correlation and whether the toxicophore model was really a significant one or a mere outcome of chance only. Accordingly, we have calculated the metric cRp2 using the following formula56 for both QSAR and toxicophore models: c 2 Rp

(1)

■

RESULTS Results Obtained from Classification-Based QSAR Analysis. As the initial number of descriptors was too high, we have applied a descriptor-thinning strategy prior to the classification approach. The linear discriminant analysis (LDA) was applied for each class of descriptors to find out the best discriminating features for individual classes. Including log KOW, a standardized pool of 10 descriptors was selected for the final discriminant analysis. As the values of the descriptors varied widely in magnitude, to determine the relative importance of a given descriptor easily from the magnitude of its coefficient, each descriptor set (discrimination as well as the test set) was scaled within the range (0 to 1) according to the following equation: D′ =

D − Dmin Dmax − Dmin

(2)

In eq 2, D′ is the scaled descriptor, D is the unscaled value of the descriptor, and Dmax and Dmin are the maximum and minimum values of the particular descriptor considering all training set compounds. LDA was performed using stepwise method of variable selection with objective function F = 4 for inclusion, F = 3.9 for exclusion with a tolerance value of 0.01, followed by the discriminant analysis. The calculated and predicted classification categories of the discrimination and test sets respectively were determined at the 50% probability level. The discriminant function ΔP is represented with the following equation:

rm 2(LOO) and Δrm2(LOO) parameters

were used for the training set and,

R2 − R r 2

For an acceptable model, the value of cRp2 should be more than 0.5. Applicability Domain (AD) Test. According to the OECD principle 3, a QSAR model should have a defined domain of applicability. Technically, AD represents the chemical space defined by the structural information of the chemicals used in model development, i.e., the training set compounds in a QSAR analysis. Here, we have tried three different approaches to assess AD. The applicability domain of the models were checked using (a) the leverage approach,57 (b) the DModX (distance to the model in X-space) approach,58 and (c) the Euclidean distance approach.59

between the values of the predicted and the corresponding observed activity data of the test set. It has been shown that for an acceptable model, the value of Δrm2(test) should preferably be lower than 0.2 provided that the value of rm 2(test) is more than 0.5.49 Similarly,

=R×

rm 2(overall) and Δrm2(overall)

were used for the overall set.49 The rm2 metrics have been used widely by our group50,51 as well as different other groups52,53 to check predictive ability of QSAR models. The models were also subjected to additional external validation parameters such as 17650



Article

discrimination and test sets are represented in Figure S2a,b, Supporting Information. The discrimination of toxicity of wastewater toxicants was carried out to show that the obtained discriminant function values in terms of % probability activity (PA) of the LDA model for both the classes make it possible to identify the two populations. In order to design the PDD, we observed that the maximum of the Ei (expectancy to get lower toxicants) and Ea (expectancy to get higher toxicants) values are distributed on different sides of PA = 50%. We obtained higher toxicants (with a PA value 50% or more) and lower toxicants (with a PA value less than 50%), in both the discrimination and the test groups. On analysis of Figure 1, although overlapping of Ei can be seen in the Ea region as well as vice versa for the discrimination set of compounds, the overlapping block is extremely lower in heights. In the case of the test set, not a single Ea block is found in the Ei region though two lower heights blocks are found in the active zone. The lower overlapping between Ei and Ea signifies meaningful PDD for a robust classification model. The studied PDD for both the sets are significant and distinctly dependable. Using the PDD, it is possible to discriminate between two groups within a structurally heterogeneous set of compounds, and it constitutes a valuable tool in the validation of discrimination analysis for our study. A contribution plot (Figure 2) was developed for the identification of best discriminating descriptors (eq 3) by taking

ΔP = (25.799 × log K OW ) + (14.638 × 3χpv ) + (20.617 × S_dCH 2) + (10.785 × 3χcv ) − 14.193 NTr = 73, λ = 0.395, F(df = 4.68) = 25.960; (p < 0.0000), R c = 0.777, Squared_Mahalanobis_distance = 5.942, MCCTr = 0.836, AUROCTraining = 0.974, χ 2 (df = 4) = 63.967; (p < 0.0000); ρ = 18.25; G‐means Training = 0.918, Cohen’s_κTraining = 0.840, NTest = 31, MCCTest = 0.808, AUROCTest = 0.962, G‐means Test = 0.901, Cohen’s_κTest = 0.810, ROCED = 0.312, (3)

ROCFIT = 0.788

KOW,3χvp,

The LDA analysis yielded a four-parameter (log S_dCH2, and 3χvc) discriminant equation with encouraging values of statistical metrics. The statistics for the different validation metrics strongly accounted for the significance of the derived model. All the metrics were within the acceptable limit and thus reflected the reliability and acceptability of the LDA model. The highly acceptable values (Table 1) were obtained Table 1. Qualitative Prediction of Classification-Based QSAR and 3D-Toxicophore Model In silico tools

metrics

training set (%)

test set (%)

classification QSAR (LDA)

sensitivity specificity precision accuracy F-measure sensitivity specificity precision accuracy F-measure

91.89 91.67 91.89 91.78 91.89 68.75 93.33 91.67 80.65 78.57

93.75 86.67 88.24 90.32 90.91 56.76 86.11 80.77 71.23 66.67

3D toxicophore

for both the training and test sets which allowed evaluation of the model’s ability to classify correctly higher and lower toxicants and to provide a complete view of the classification models’ predictive abilities. The ROC curves for the

Figure 2. Average contributions of indices to the discriminant functions for higher and lower toxic groups.

Figure 1. PDD for the discrimination set and the test set. 17651



Article

Figure 3. Test of AD by Euclidean distance method (mean normalized distance from the training set compounds is shown) for (a) classification based QSAR approach and (b) regression based QSAR approach.

Table 2. Statistical Qualities of Different Regression-Based QSAR Modelsa model stepwise regression GFA-linear GFA-spline G/PLS-linear G/PLS-spline a

number of descriptors

R2

Q2

0.713

0.666

3

0.714 0.723 0.676 0.736

0.654 0.673 0.641 0.685

4 4 3 (LV = 1) 4 (LV = 2)

2 rm(LOO)

Δr2m(LOO)

r2

R2pred

0.543

0.204

0.736

0.719

0.530 0.551 0.508 0.565

0.189 0.200 0.236 0.194

0.701 0.705 0.687 0.712

0.684 0.687 0.682 0.693

2 rm(test)

Δr2m(test)

QF22

GTCb

2 rm(overall)

0.623

0.216

0.716

passed

0.559

0.211

0.588 0.589 0.568 0.602

0.173 0.218 0.197 0.168

0.678 0.683 0.677 0.689

passed passed passed passed

0.539 0.589 0.521 0.569

0.187 0.218 0.230 0.189

Δr2m(overall)

The best model is bold faced. bGolbraikh and Tropsha’s criteria.

using different chemometric tools. A detailed report of the statistical quality of various models is elaborated in Table 2. Among the various developed models, the best results in terms of internal and external validation were obtained from the G/ PLS spline method. The developed equation is as follows:

the product of their average values and their corresponding coefficients. Analyzing the contribution plot, we can say that the index octanol−water partition coefficient (log KOW) makes a significant positive contribution toward the toxicity. The second most important discriminating descriptor is 3χvp which also positively contributes toward the toxicity. The parameter 3χvp is a topological descriptor belonging to the category of molecular connectivity indices and defined as the third-order valencemodified path connectivity index which indicates the impact of branching and shape on the toxicity of the studied molecules. Indices such as 3χvc and CH2 fragments also have positive contributions toward the toxicity, but the contributions are significantly lower than those of the previous two descriptors. The parameter 3χvc is a topological descriptor belonging to the category of molecular connectivity indices and defined as the third-order cluster index based on valence count. It encodes the number of branch points in a molecule and indicates that a decrease in branching of the molecule leads to decreases in its toxicity value. The S_dCH2 descriptor refers to the summation of E-state values for the CH2 fragments present within the molecular structure. Thus, the presence of this descriptor signifies the influence of this structural fragment for higher toxicity profile. Therefore, CH2 fragments or unsaturation should be avoided for lower toxicity profile. The AD of the discriminant model has been studied by using a diversity validation measure performed by a program EUCLIDEAN59 developed in our laboratory. A scatter plot of the mean normalized distance of the compounds (Figure 3a) showed that not a single toxicant falls outside of the AD of the model. Thus, predictions for 100% of the test set compounds are highly reliable. Results Obtained from the Regression-Based QSAR Analysis. Regression based QSAR models were developed

pEC50 = 3.039 − 0.411 × < 2.93 − log K OW > −0.523 × < 1.512 − 3χpv > + 3.435 × + 0.621 × Atype − H − 49 2 n Training = 73, LV = 2, R2 = 0.736, Q LOO = 0.685, 2 2 rm(LOO)Scaled = 0.565, Δrm(LOO)Scaled = 0.194, n Test 2 2 2 = 31, Q F1 = R pred = 0.693, Q F2 = 0.689, 2 2 rm(test)Scaled = 0.602, Δrm(test)Scaled = 0.168, 2 2 rm(overall)Scaled = 0.569, Δrm(overall)Scaled = 0.189

(4)

Equation 4 involving four descriptors and only two latent variables (LVs) could explain 73.6% of the variance. Though the total data set is extremely diverse both in terms of structure and mechanism of actions of the molecules, satisfactory measures of goodness-of-fit, robustness, and predictivity for eq 4 were obtained in compliance with OECD principles. Moreover, least possible deviation of the predicted activity data from the corresponding observed ones is further implied from the satisfactory values of all the rm2 metrics. Quite identical values for the Qext(F1)2 (0.685) and Qext(F2)2 (0.689) metrics indicate that the test set selected for the QSAR model development has similar distribution of response as the training 17652



Article

Figure 4. Williams plot for the G/PLS spline model.

Like the classification-based approach, 3χvp is evolved as the second important descriptor. As previously mentioned, the parameter 3χvp is a topological descriptor belonging to the category of molecular connectivity indices and defined as the third-order valence-modified path connectivity index which indicates the impact of branching and shape on the toxicity of the molecules. According to the spline term and its coefficient, if the 3χvp index has a value of 1.512 or more, the corresponding compounds will be classified into higher toxicity ones, and it is true for all 17 (for example, 8-quinolinol, 2,3,6-trichlorophenoand, and 2,4,5-trichlorophenol) compounds which have 3χvp values of more than 1.512 in our study. In order to lower the toxicity of a particular compound, the 3χvp value should be less than 1.512. A lower value will lead to lower toxicity. This suggests a lower degree of branching, smaller size, and presence of a smaller number of heteroatoms as the essential structural attributes for lower toxicity. The Σβs′ (= [Σβs]/NV) is defined as a relative measure of number of electronegative atoms in the substructure.61 Σβs is calculated considering sigma bonds in the substructure and vertex count excluding hydrogens (NV). The index bearing a value of less than 0.438 lowers the toxicity. The index Atype− H−49 signifies hydrophobicity of H attached to C3sp3, C2−3sp2, C1−3sp where the subscripts represents hybridization and the superscript represents formal oxidation number. The forms of hydrogen may be represented as follows:

set. External predictivity was further judged by Golbraikh and Tropsha’s criteria (GTC). For the G/PLS spline model, these statistical parameters yielded the following results which are highly acceptable: 2 R CV,ext = 0.685, r 2 = 0.712, |r02 − r0′ 2| = 0.082, 2

= 0.001, or

r − r2

r′20

r 2 − r02 r2

= 0.116, 0.85 ≤ k = 1.04 ≤ 1.15

or 0.85 ≤ k′ = 0.94 ≤ 1.15

The variable importance projections (VIPs) and the coefficient histogram of the original descriptors for G/PLS spline model are presented in Figures S3 and S4 in Supporting Information. In order to comply with the OECD Principle 5, we explain the mechanistic interpretation and importance of each descriptor appearing in the best regression equation with suitable examples. The descriptors are ranked according to the descending VIP values: , , and Atype−H −49. Decimal logarithm of octanol−water partition coefficient (log KOW) has a direct nonlinear relationship with the toxicity and evolved as one of the major predictive descriptors for the regression-based QSAR model. A strong relationship between toxicity and log KOW has been observed in lower class of animals in previous reports also.50,60 The spline term of partition coefficient has a negative contribution toward the aquatic toxicity. A partition coefficient value of 2.93 or more will make the spline term zero, and the contribution of this descriptor will be null. Then the toxicity value will be maximum for these compounds. As the value decays from 2.93, the toxicity value will be lower. So, one can conclude that for lower toxicity, partition coefficient should be less than 2.93. Compounds such as triethanolamine (49), ethanol (65), methanol (67), and formamide (79) have partition coefficient values in negative range, which makes them lower toxicants within the studied molecules.

where X is any electronegative atoms and R3 is C or H or X. The presence of any of these fragments and their contributions toward hydrophobicity make a compound more toxic. Higher toxicity of compounds such as 4-chlorobenzaldehyde, 2chlorobenzaldehyde, and bromoform can be demonstrated by the presence of the mentioned fragments. Randomization test was also done by randomly reordering (100 permutations) response data using SIMCA-P 10.0. The 17653



Article

best ranking three feature (HYD Aliphatic, HYD Aliphatic, and HYD Aromatic) toxicophore based on the values of the correlation coefficient (0.724) and the cost functions (total cost = 156.254). The values of fixed cost and null cost, expressed in bits, differed significantly from each other by 105.610 bits, and such a difference implied the existence of more than 95% chance of true correlation for the developed 3D toxicophore. Again, the difference between total cost and null cost was 45.563, which is quite acceptable for a robust toxicophore model. The configuration cost parameter yielded a value of 5.087 bits, which was much lower than the threshold value of 17, and the rms deviation for the hypothesis was 1.853. The robustness of the developed toxicophore model was also evaluated based on qualitative validation parameters (Table 1). The Fischer validation test checked at the 95% confidence level also showed that the value of the average correlation coefficient of the randomized models (Rr) was much lower than the corresponding correlation coefficient (R) of the nonrandomized model (Rr = 0.204, R = 0.932), implying the robustness of the developed model and denoting the existence of a true correlation with a cRp2 value of 0.503. Again, total cost of the randomized model (195.165) was close to the null cost (201.817) of the original model, which proved that the developed model was not obtained by chance. External validation of the developed model was performed by mapping the test set molecules, and the predicted activity data was matched with those of the observed ones based on the classification technique. The ability of the model to ideally distinguish between the two classes of the test set was determined based on different qualitative validation parameters (Table 1). Acceptable values for all these validation parameters ensured the predictive potential of the developed 3D toxicophore model and reflected the statistical significance of the developed model. Hypothesis 1 is a three-feature toxicophore displaying the significance of the following features arranged at specific distances: HYD Aliphatic, HYD Aliphatic, and HYD Aromatic. The hydrophobic aliphatic and the hydrophobic aromatic features denote the regions favorable for substitution by aliphatic hydrophobic groups and aromatic hydrophobic groups, respectively. A triangle having a definite shape is obtained by placing the three features at the vertices. The distances between the two consecutive HYD Aliphatic features is 6.515 Å, while the HYD Aromatic feature is placed at a distance of 4.104 Å and 3.595 Å from the two hydrophobic aliphatic features. Two HYD Aliphatic features and HYD Aromatic feature are placed at an angle of 113.984°. Molecules bearing substitutions that match with the vertices of the triangle thus formed and thereby capture all the toxicophoric features may exhibit maximum toxicity. Figure 6a shows the mapping of one of the most toxic compounds 2,3,4,5-tetrachlorophenol (46) to the best hypothesis. The phenyl ring and two chlorine atoms at the 3 and 5 positions of the 2,3,4,5-tetrachlorophenol molecule capture the HYD Aromatic and two HYD Aliphatic features respectively. Thus, the molecule capturing all the features exerts maximum toxicity. On the contrary, compound nos. 59, 68, and 76 lacking the detrimental substituents fail to map with all the toxicophoric features which in turn account for their least toxic character. Figure 6b shows the mapping of one of the least toxic compounds (76) to the developed toxicophore. Comparison and Consensus Results of All Three Models. The molecular attributes required for the toxicity

randomization parameter of the best model is well above the permissible limit, indicating that the model is not obtained by chance. The toxicity intercept values are R2 = (0.0, −0.007), Q2 = (0.0, −0.172). The randomization plot is presented in Figure S5 in Supporting Information. The lack of chance correlation in the G/PLS spline model is also well reflected from the value of c 2 Rp (0.729), which is higher than the acceptable threshold value of 0.5. Both results suggest that the obtained model is not derived by chance. The applicability domain for the G/PLS spline model was checked using the leverage, DModX, and diversity validation approaches. The leverage approach showed a critical HAT diagonal value (h*) of 0.205, and we have used a ± 3σ standard deviation unit of the predicted residual values for defining the domain of applicability of the model. Apart from acetic acid (66), the remaining of all the training set compounds bear values of standardized residual within the limit of ±3σ indicating that only one compound was a prediction outlier. However, leverage values of formaldehyde (47) of the training set being greater than the critical value of 0.205 (h > h*), this compound behaves as an influential observation (X outlier) although it is not a response outlier (not Y outlier). Out of the 31 test set compounds, 27 test compounds were found to be within the domain of applicability. Chemicals chloroform (3), glutaraldehyde (25), carbon tetrachloride (53), and 1,1,1trichloroethane (73) were (< ± 3σ) completely outside the AD of the model, as defined by the Hat vertical line (high h leverage value). Further, at the 95% confidence level, DModX values of 26 test compounds were within the critical value of 1.986 (chemical 10 [1,2-dibromoethane] along with compound number 3, 25, 53, and 73 are outside the critical value), though in the case of the diversity validation approach employing Euclidean distance measures, all the compounds were found to lie within the domain defined by the training set. Considering all processes of AD tests, we conclude that 26 test compounds were inside of the AD, and their predictions are highly reliable. Therefore, we can confidently predict 83.9% of test compounds based on the developed model after rigorous tests for validation and applicability domain check. Euclidean, Williams, and DModX plots are presented in Figures 3b, 4, and 5 respectively. Analysis of the 3D Toxicophore Model. Ten toxicophore hypotheses were based on the conformers obtained from the BEST method of conformer generation using the training set compounds. Subsequently, hypothesis 1 was selected as the

Figure 5. DModX plot of test set compounds at 95% confidence level of the best G/PLS-spline model. 17654



Article

atoms at the 3 and 5 position of the 2,3,4,5-tetrachlorophenol molecule which is one of the most toxic compounds capture the HYD Aromatic and two HYD Aliphatic features, respectively. Here also, highly electronegative chlorine atoms and their contribution toward hydrophobicity are mainly responsible for higher toxicity profile of this molecule. Fragments such as CH2 and unsaturation are depicted by the aromatic feature of the molecule. Compounds such as ether (59), 2-propanol (68), and bromoethane (76) lack unsaturation, higher degree of branching, and higher number of electronegative atoms to map the requisite features of the toxicophore. As a result, they show a lower toxicity profile.

■

CONCLUSIONS We have presented QSAR and toxicophore models for the prediction of toxicity and identification of required structural features and fragments of a diverse class of water toxicants. Like earlier studies, lipophilicity is proved to be one of the major descriptors to define the toxicity of these toxicants. Again, it is also true that lipophilicity is not a sole descriptor for toxicity prediction. The third-order branching, flexibility, shape of the molecule, CH2 fragments, and lipophilic contribution for particular fragments present in toxicants are other important predictors for these classes of chemicals which also support the previously studied results. These models therefore provide an understanding of important structural attributes of molecules and may provide important information for the design of chemicals for synthesis of future molecules with diminished systemic toxicity profile. Moreover, the 3D toxicophore model developed in the present work may serve as an efficient 3D query tool for screening of large databases and identifying the molecules bearing the detrimental structural features associated with their toxicity potential. It is conceivable that our in silico approach offers an efficient screening method for toxicity and can be used in integration with existing in vitro methods to increase overall predictivity.

Figure 6. (a) Mapping of one of the most toxic compounds (compound no. 46) to the developed toxicophore and (b) mapping of one of the least toxic compounds (compound no. 76) to the developed toxicophore. Hydrophobic aliphatic group and hydrophobic aromatic group features are represented with blue and red colors respectively.

■

ASSOCIATED CONTENT

S Supporting Information *

and their degree of contribution to aquatic toxicity of the molecules were determined from QSAR models. Additionally, the crucial features constituting the toxicophore of the molecules were analyzed based on the developed 3D toxicophore model. Satisfyingly, the results obtained from the regression and classification-based QSAR and toxicophore models were in compliance with each other with respect to the obtained requisite features and structural fragments that enable one to discriminate between higher and lower aquatic toxic potency of the studied chemicals. Summarizing both the classification and regression-based QSAR models, features such as octanol−water partition coefficient, third-order branching, CH2 fragments or unsaturation and the presence of higher number of electronegative atoms (specifically halogen atoms) and their contribution toward hydrophobicity are identified majorly responsible structural attributes for higher toxicants. For lower toxicity, the log KOW and 3χvp value should be lower than 2.93 and 1.512, respectively which can be well inferred from the classification-based QSAR model. Additionally, one should avoid the following fragments:

Figure S1. Normality distribution diagram. Figure S2. Receiver operating characteristics (ROC) curves for the discrimination set and the test set. Figure S3. Histogram of VIPs of the descriptors used in the best G/PLS spline model. Figure S4. Histogram of coefficients of the original descriptors used in the best G/PLS spline model. Figure S5. Validation of the best G/ PLS spline model. Table S1. Results obtained from classification and regression based QSARs and toxicophore model for aquatic toxicity. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected], [email protected], or [email protected]. URL: http://sites.google.com/ site/kunalroyindia/. Phone: +91-98315 94140. Fax: +91-332837 1078. Present Address †

K.R. is currently at Manchester Institute of Biotechnology, Manchester M1 7DN, United Kindom Notes

and lower the number of electronegative substitutions for lower toxicity profile of chemicals. The phenyl ring and two chlorine

The authors declare no competing financial interest. 17655



■

Article

Washington, DC, 2009; (accessed 21 August 2012). (19) Kar, S.; Roy, K. Predictive toxicology using QSAR: a perspective. J. Indian Chem. Soc. 2010, 87, 1455−1515. (20) Kar, S.; Roy, K. Risk assessment for ecotoxicity of pharmaceuticals-an emerging issue. Expert. Opin. Drug Saf. 2012, 11, 235−274. (21) Xu, S.; Nirmalakhandan, N. Use of QSAR models in predicting joint effects in multi- component mixtures of organic chemicals. Water Res. 1998, 32, 2391−2399. (22) Schultz, T. W.; Ralston, K. E.; Roberts, D. W.; Veith, G. D.; Aptula, A. O. Structure- activity relationships for abiotic thiol reactivity and aquatic toxicity of halosubstituted carbonyl compounds. SAR QSAR Environ. Res. 2007, 18, 21−29. (23) Katritzky, A. R.; Kasemets, K.; Slavov, S.; Radzvilovits, M.; Tämm, K.; Karelson, M. Estimating the toxicities of organic chemicals in activated sludge process. Water Res. 2010, 44, 2451−2460. (24) Christen, V.; Hickmann, S.; Rechenberg, B.; Fent, K. Highly active human pharmaceuticals in aquatic systems: a concept for their identification based on their mode of action. Aquat. Toxicol. 2010, 96, 167−81. (25) Kar, S.; Roy, K. First report on interspecies quantitative correlation of ecotoxicity of pharmaceuticals. Chemosphere 2010, 81, 738−47. (26) Guidance Document on the Validation of (Quantitative) 1226 Structure-Activity Relationships (Q)SARs] Models; ENV/JM/ MONO(2007)2; OECD Environment Health and Safety Publications Series on Testing and Assessment, No. 69; Organisation for Economic Co-operation and Development: Paris, 2007. (27) Cerius 2, Version 4.10 Software; Accelrys Inc.: San Diego, CA, USA; http://www.accelrys.com/cerius2. (28) DRAGON, ver. 6 software; TALETE srl: Italy, http://www. talete.mi.it/products/dragon_molecular_descriptors.htm. (29) Yap, C. W. PaDEL-Descriptor: An open source software to calculate molecular descriptors and fingerprints. J. Comput. Chem. 2011, 32, 1466−1474. (30) Darlington, R. B. Regression and Linear Models; McGraw-Hill: New York, 1990. (31) Fan, Y.; Shi, L. M.; Kohn, K. W.; Pommier, Y.; Weinstein, J. N. Quantitative structure−antitumor activity relationships of camptothecin analogues: cluster analysis and genetic algorithm-based studies. J. Med. Chem. 2001, 44, 3254−3263. (32) Rogers, D.; Hopfinger, A. J. Application of Genetic Function Approximation to Quantitative Structure Activity Relationships and Quantitative Structure Property Relationships. J. Chem. Inf. Comput. Sci. 1994, 34, 854−866. (33) Mitteroecker, P.; Bookstein, F. Linear Discrimination, Ordination, and the isualization of Selection Gradients in Modern Morphometrics. Evol. Biol. 2011, 38, 100−114. (34) Smellie, A.; Teig, S. L.; Towbin, P. Poling: promoting conformational variation. J. Comput. Chem. 1995, 16, 171−187. (35) Discovery Studio 2.1; Accelrys Inc: SanDiego, CA, 2010. (36) Poptodorov, K.; Luu, T.; Hoffmann, R. D. In Methods and Principles in Medicinal Chemistry, Pharmacophores and Pharmacophores Searches; Langer, T.; Hoffmann, R. D., Eds.; Wiley-VCH: Weinheim: Germany, 2006; Vol. 2, pp 17−47. (37) STATISTICA Statistical Software; StatSoft Inc.: Tulsa, OK; http://www.statsoft.com/. (38) SPSS Statistical Software; SPSS Inc.: Chicago, IL; http://www. spss.com. (39) MINITAB Statistical Software; Minitab Inc.: State College, PA; http://www.minitab.com. (40) SIMCA-P 10.0; UMETRICS: Umea, Sweden, 2002; www. umetrics.com. (41) Roy, K.; Mitra, I. On Various Metrics Used for Validation of Predictive QSAR Models with Applications in Virtual Screening and Focused Library Design. Comb. Chem. High Throughput Screen 2001, 14, 450−474.

ACKNOWLEDGMENTS S.K. thanks the Department of Science and Technology, Government of India, for awarding him a Research fellowship under the INSPIRE scheme. K.R. thanks the Council of Scientific and Industrial Research (CSIR), New Delhi, for awarding a major research project.

■

REFERENCES

(1) Chelliapan, S.; Wilby, T.; Yuzir, A.; Sallis, P. J. Influence of organic loading on the performance and microbial community structure of an anaerobic stage reactor treating pharmaceutical wastewater. Desalination 2011, 271, 257−264. (2) Chang, Ch. Y.; Chang, J. S.; Vigneswaran, S.; Kandasamy, J. Pharmaceutical wastewater treatment by membrane bioreactor process − a case study in southern Taiwan. Desalination 2008, 234, 393−401. (3) Vieno, N.; Tuhkanen, T.; Kronberg, L. Elimination of pharmaceuticals in sewage treatment plants in Finland. Water Res. 2007, 41, 1001−1012. (4) Zeilinger, J.; Steger-Hartmann, T.; Maser, E.; Goller, S.; Vonk, R.; Länge, R. Effects of synthetic gestagens on fish reproduction. Environ. Toxicol. Chem. 2009, 28, 2663−2670. (5) Noble, J. GE ZeeWeed MBR technology for pharmaceutical wastewater treatment. Membr. Technol. 2006, 2006, 7−9. (6) Dutka, B. J.; Kwan, K. K. Application of four bacterial screening procedures to assess changes in the toxicity of chemical in mixtures. Environ. Pollut. 1982, 29, 125−134. (7) Stumpf, M.; Ternes, T. A.; Wilken, R. D.; Rodrigues, S. V.; Baumann, W. Polar drug residues in sewage and natural waters in the state of Rio de Janeiro, Brazil. Sci. Total Environ. 1999, 225, 135−141. (8) Wong, K. Y.; Zhang, M. Q.; Li, X. M.; Lo, W. A luminescencebased scanning respirometer for heavy metal toxicity monitoring. Biosens. Bioelectron. 1997, 12, 125−133. (9) Carballa, M.; Omil, F.; Lema, J. M.; Llompart, M.; García-Jares, C.; Rodriguez, I.; Gómez, M.; Ternes, T. Behavior of pharmaceuticals, cosmetics and hormones in a sewage treatment plant. Water Res. 2004, 38, 2918−2926. (10) Zorita, S.; Martensson, L.; Mathiasson, L. Occurrence and removal of pharmaceuticals in a municipal sewage treatment system in the south of Sweden. Sci. Total Environ. 2009, 407, 2760−2770. (11) Ren, S.; Frymier, P. D. Estimating the toxicities of organic chemicals to bioluminescent bacteria and activated sludge. Water Res. 2002, 36, 4406−4414. (12) Ren, S.; Frymier, P. D. Comparative study of two bioassays for applications in influent wastewater toxicity monitoring. J. Environ. Eng. 2003, 129, 216−221. (13) Ren, S.; Frymier, P. D. Toxicity estimation of phenolic compounds by bioluminescent bacterium. J. Environ. Eng. 2003, 129, 328−335. (14) Dearden, J. C. Prediction of Environmental Toxicity and Fate Using Quantitative Structure- Activity Relationships (QSARs). J. Braz. Chem. Soc. 2002, 13, 754−762. (15) Directive 2006/121/EC of the European Parliament and of the Council of 18 December 2006 amending Council Directive 67/548/EEC on the approximation of laws, regulations and administrative provisions relating to the classification, packaging and labelling of dangerous substances in order to adapt it to Regulation (EC) No. 1907/2006 concerning the REACH and establishing a European Chemicals Agency; European Commission, Official Publications of the European Communities (OPOCE): Luxembourg, 2006. (16) Benigni, R.; Giuliani, A. Putting the Predictive Toxicology Challenge into Perspective: Reflections on the Results. Bioinformatics 2003, 19, 1194−1200. (17) Williams, E. S.; Panko, J.; Paustenbach, D. J. The European Union’s REACH regulation: a review of its history and requirements. Crit. Rev. Toxicol. 2009, 39, 553−675. (18) Integrated Risk Information System; U.S. Environmental Protection Agency, National Center for Environmental Assessment: 17656



Article

(42) Cohen, J. A. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 1960, 20, 37−46. (43) Gálvez-Llompart, M.; Recio, M. C.; García-Domenech, R. Topological virtual screening: a way to find new compounds active in ulcerative colitis by inhibiting NF-kB. Mol. Divers. 2011, 15, 917−926. (44) Prado-Prado, F. J.; Uriarte, E.; Borges, F.; González-Díaz, H. Multi-target spectral moments for QSAR and complex networks study of antibacterial drugs. Eur. J. Med. Chem. 2009, 44, 4516−4521. (45) Fawcett, T. An introduction to ROC analysis. Pattern Recogn. Lett. 2006, 27, 861−874. (46) Perez-Garrido, A.; Helguera, A. M.; Borges, F.; Cordeiro, M. N. D. S.; Rivero, V.; Escudero, A. G. Two New Parameters Based on Distances in a Receiver Operating Characteristic Chart for the Selection of Classification Models. J. Chem. Inf. Model 2011, 51, 2746− 2759. (47) Murcia-Soler, M.; Pérez-Giménez, F.; Garćıa-March, F. J.; Salabert-Salvador, M. T.; D́ ıaz-Villanuev, W.; Medina-Casamayor, P. Discrimination and selection of new potential antibacterial compounds using simple topological descriptors. J. Mol. Graphics Model 2003, 21, 375−390. (48) Ojha, P. K.; Mitra, I.; Das, R. N.; Roy, K. Further exploring rm2 metrics for validation of QSPR models. Chemom. Intell. Lab. Syst. 2011, 107, 194−205. (49) Roy, K.; Mitra, I.; Kar, S.; Ojha, P.; Das, R. N.; Kabir, H. Comparative studies on some metrics for external validation of QSPR models. J. Chem. Inf. Model. 2012, 52, 396−408. (50) Kar, S.; Harding, A. P.; Roy, K.; Popelier, P. L. A. QSAR with quantum topological molecular similarity indices: toxicity of aromatic aldehydes to Tetrahymena pyriformis. SAR QSAR Environ. Res. 2010, 21, 149−168. (51) Kar, S.; Roy, K. Development and validation of a robust QSAR model for prediction of carcinogenicity of drugs. Indian J. Biochem. Biophys. 2011, 48, 111−122. (52) Shahlaei, M.; Sabet, R.; Ziari, M. B.; Moeinifard, B.; Fassihi, A.; Karbakhsh, R. QSAR study of anthranilic acid sulfonamides as inhibitors of methionine aminopeptidase-2 using LS-SVM and GRNN based on principal components. Eur. J. Med. Chem. 2010, 45, 4499− 4508. (53) Toropova, A. P.; Toropov, A. A.; Lombardo, A.; Roncaglioni, A.; Benfenati, E.; Gini, G. A new bioconcentration factor model based on SMILES and indices of presence of atoms. Eur. J. Med. Chem. 2010, 5, 4399−4402. (54) Schüürmann, G.; Ebert, R.-U.; Chen, J.; Wang, B.; Kühne, R. External validation and prediction employing the predictive squared correlation coefficient − test set activity mean vs. training set activity mean. J. Chem. Inf. Model 2008, 48, 2140−2145. (55) Golbraikh, A.; Tropsha, A. Beware of q2! J. Mol. Graphics Model 2002, 20, 269−276. (56) Mitra, I.; Saha, A.; Roy, K. Exploring quantitative structureactivity relationship (QSAR) studies of antioxidant phenolic compounds obtained from traditional Chinese medicinal plants. Mol. Simul. 2010, 36, 1067−1079. (57) Gramatica, P. Principles of QSAR Models Validation: Internal and External. QSAR Comb. Sci. 2007, 26, 694−701. (58) Wold, S.; Sjostrom, M.; Eriksson, L. PLS-regression: a basic tool of chemometrics. Chemom. Intell. Lab. Syst. 2001, 58, 109−130. (59) EUCLIDEAN (a program written in Java) is developed and validated on known data sets by Pravin Ambure (Email: ambure. [email protected]) of Drug Theoretics and Cheminformatics Laboratory, Jadavpur University, 2013. (60) Schultz, T. W.; Sinks, G. D.; Bearden, A. P. QSAR in aquatic toxicology: a mechanism of action approach comparing toxic potency to Pimephales promelas, Tetrahymena pyriformis, and Vibrio fischeri. In Comparative QSAR; Devillers, J., Ed.; Taylor & Francis: London, 1998; pp 51−109. (61) Roy, K.; Das, R. N. QSTR with extended topochemical atom (ETA) indices. 15. Development of predictive models for toxicity of organic chemicals against fathead minnow using second generation ETA indices. SAR QSAR Environ. Res. 2012, 23 (606), 125−140. 17657


Predictive Chemometric Modeling and Three-Dimensional

Recommend Documents