Article pubs.acs.org/est
Iterative Fragment Selection: A Group Contribution Approach to Predicting Fish Biotransformation Half-Lives Trevor N. Brown,†,‡,§ Jon A. Arnot,‡ and Frank Wania*,†,‡ †
Department of Chemistry, University of Toronto Scarborough, 1265 Military Trail, Toronto, Ontario, Canada M1C 1A4 Department of Physical and Environmental Sciences, University of Toronto Scarborough, 1265 Military Trail, Toronto, Ontario, Canada M1C 1A4
‡
S Supporting Information *
ABSTRACT: There are regulatory needs to evaluate thousands of chemicals for potential hazard and risk with limited available information. An automated method is presented for developing and evaluating Quantitative Structure−Activity Relationships (QSARs) for a range of chemical properties that can be applied for screening level chemical assessments. The method is an integrated algorithm for descriptor generation, data set splitting, cross validation, and model selection. Resulting QSARs are two-dimensional (2D) fragmentbased group contribution models. The QSAR development and evaluation method does not require previous expert knowledge for selecting 2D fragments associated with the chemical property of interest. The method includes information on the domain of applicability (structural similarity to the training set) and estimates of the uncertainty in the QSAR predictions. As a demonstration, the method is applied to generate novel QSARs for fish primary biotransformation half-lives (HLN). Results from the new HLN QSARs are compared to another 2D fragment-based HLN QSAR developed with expert judgment, and the predictive powers of the models are found to be similar. The relative merits and limitations of each method are investigated and the new QSAR is found to make comparable predictions with significantly fewer fragments. A coefficient of determination (R2) of 0.789 and a root mean squared error (RMSE) of 0.526 were obtained for the training data set and an R2 of 0.748 and an RMSE of 0.584 were obtained for the validation data set, along with a concordance correlation coefficient (CCC) of 0.857 showing good predictive power.
■
INTRODUCTION National and international regulatory programs and a global treaty mandate chemical hazard and risk assessments to protect human and environmental health.1−3 Chemical assessments require monitoring, exposure, toxicity, and chemical property information such as the octanol−water partition coefficient (KOW) and reaction half-lives (HLs). Compared to the number of chemicals legislated for assessment there is a paucity of available measured data, making it necessary to use theoretical approaches to estimate chemical properties, exposure, and toxicity.4−6 Chemical distribution properties and HLs can be predicted from two-dimensional (2D) structure using quantitative structure−property relationships (QSPRs) or quantitative structure−activity relationships (QSARs) such as those implemented in the U.S. Environmental Protection Agency’s Estimation Program Interface (EPI) Suite7 and elsewhere.8,9 More sophisticated methods for predicting chemical properties have also been developed but are often proprietary or difficult to apply to very large data sets.10−12 For some physical chemical properties, such as KOW, multiple prediction tools are available; however, the uncertainty and variability between predictions can be substantial.13,14 For other chemical parameters, such as HLs, the availability of high throughput screening-level tools is very limited.7 There is thus a need to develop, evaluate, and compare more QSARs to address data © 2012 American Chemical Society
gaps, identify possible sources of uncertainty in training and testing data sets and in QSAR methodologies, and quantify and reduce uncertainty in applying these predictive methods for chemical assessments. Biotransformation HLs are important chemical parameters for calculating ecological and human exposure to chemicals.15−18 In general, the influence of biotransformation on human exposures is much greater than the influence of chemical partitioning properties.19 Currently, there is only one publicly available QSAR for whole body biotransformation HLs in vertebrates (fish).7,20 Experimental data and expert knowledge for whole body biotransformation HLs are limited and, due to the complex and variable nature of this process, also highly uncertain.21 Expert knowledge and “rules of thumb” of environmental (microbial) degradation mechanisms and recalcitrant chemical groups have been built into chemical class specific models (see review by Pavan and Worth22) and more general models such as the BIOWIN set of models in EPI Suite.23 Mechanistic knowledge of environmental microbial degradation was initially applied in the development and Received: Revised: Accepted: Published: 8253
March 28, 2012 June 5, 2012 July 10, 2012 July 10, 2012 dx.doi.org/10.1021/es301182a | Environ. Sci. Technol. 2012, 46, 8253−8260
Environmental Science & Technology
Article
evaluation of a fish biotransformation HL-QSAR and the results suggest some overlap in the importance of selected chemical descriptors (fragments) and biotransformation mechanisms in microbes and fish.7,20 The success in extrapolating this knowledge is encouraging; however, it is possible that other mechanisms and molecular fragments contribute to the rate of chemical biotransformation, for which there is no current “expert recognition”. A tool that automatically generates hypotheses for mechanisms and molecular fragments related to a chemical property of interest could provide new and deeper insights into QSARs and expand expert knowledge. The primary objective of this research is to develop a new automated tool for the prediction of chemical properties or activities, and to apply it to the development of a fish biotransformation HL-QSAR. The overall method of fragment generation, data set splitting, and model selection is referred to as Iterative Fragment Selection (IFS). This new QSAR method does not require pre-existing expert knowledge (or “rules-ofthumb”) on structural information related to a specific chemical property. 2D fragments are generated automatically from the data set and the method iteratively selects fragments and QSARs through an automated process of fragment selection and model testing. To our knowledge only two other methods, both either commercial or patented products, generate fragments from the chemical structures of the data set without applying any information about the property to be predicted before or during the fragment generation: MULTICASE12 and Discrete Substructural Analysis.24 Predictive power of the QSARs generated is of course limited by the quality of the training data sets available, the selection and refinement of which still primarily relies on expert judgment. A new biotransformation HL-QSAR (referred to as IFS-HLN) is developed and evaluated and then compared with an existing fragment-based HL-QSAR.7,20 Similarities and differences between the QSARs are discussed with a focus on model performance and fragment selection.
After the initial fragmentation each chemical has been broken into small fragments referred to as single unit fragments. Complex fragments are subsequently created by recursively recombining the single unit fragments with their neighbors until all substructures of every chemical are in theory generated. However, in two special cases the recursive recombination is truncated: first in the case of a fragment that is unique in the data set, and second in the case of multiple fragments with perfectly collinear fragment counts in all of the chemicals in which they are present. Each represents a case where little further useful information is gained by continuing the recursive recombination. For every chemical the unique fragment(s) and the set(s) of collinear fragments are each consolidated into larger structures which contain all of the individual fragments. Further details are provided in the Supporting Information (SI), section SI1. Fragments are divided into three classes of increasing selection priority: rare fragments, normal fragments, and basic fragments. Rare fragments are defined as those fragments that occur in only one or two chemicals. In the case of crossvalidation a fragment is considered rare if it is rare in any one of the cross-validation subsets. In structural terms, basic fragments are essentially those fragments which are never contained as a subfragment in any other fragment. However the exact definition is based on fragment counts rather than structure, in order to make the definition of fragment priorities independent from the method used to generate the fragments. Finally, normal fragments are those which do not fall under the other two classifications. Further details are provided in the SI, section SI1. Model Selection. Model selection proceeds according to the IFS algorithm. Fragments are added iteratively to the multiple linear regression (MLR) starting from a model containing only the intercept (the average value of the observations). Forward selection and backward removal of fragments are based on improvements to the goodness of fit (GoF), with the fragments that produce the largest improvement in GoF given priority. Fragments are added one at a time and after each addition none or many fragments may be removed if this improves the GoF. A selected fragment may also replace one or more fragments already included in the model if selection rules prohibit the fragments from being included together. There are three modes in which fragments are added or removed according to their priority. In basic fragment mode only basic fragments are added, but all fragments may be removed. In normal mode only normal fragments may be added and only normal and rare fragments may be removed. In rare mode only rare fragments may be added or removed. Once no further improvements in GoF are obtained by adding or removing fragments in one mode the algorithm proceeds to the next mode in the order basic, then normal, then rare, and then back to basic. All fragment counts are regressed against the property to be predicted in a single multiple linear regression; no fitting is done on residuals. Goodness of Fit Metric. The Akaike Information Criterion corrected for data set size (AICC) is used to evaluate model performance and GoF.30 AICC penalizes model complexity and can be used to limit the number of fragments selected to a reasonable fraction of the number of chemicals in the training data set. The assumptions on which the AICC is based are well suited to fitting measured data with fragments derived from 2D structures; namely that a model which perfectly fits the data does not exist and instead the goal is to find the simplest model
■
METHODS Overview of the IFS Method. The methods adhere to good QSAR development practice,25 while developing an automated method of descriptor generation, data set splitting, cross validation, and model selection. In the first step of the method fragments are generated using a modified version of the OpenBabel chemistry library version 2.3.0.26,27 All potentially useful fragments contained in the chemicals of a data set are generated automatically and used as the pool of descriptors during model selection. The model selection algorithm is composed of nested iterations of fragment selection and removal. A simple set of selection rules based on fragment counts determines the priority of adding and removing fragments during model selection. Multiple linear regression and other matrix algebra operations are performed with the Armadillo library version 2.0.1,28 and other supporting code is written in C++. Generation of Descriptors. The pool of descriptors drawn for model selection is composed entirely of fragments created from the data set. All of the chemicals in the data set are fragmented by breaking single and aromatic bonds, with the exception of those bonds involving hydrogen. Double and triple bonds are not broken. Fragments are converted to canonical Daylight SMARTS notation29 to allow for rapid substructure searching in OpenBabel and other software that supports this feature. 8254
dx.doi.org/10.1021/es301182a | Environ. Sci. Technol. 2012, 46, 8253−8260
Environmental Science & Technology
Article
which neglects the least information.30 Many fragments can be generated, many of which likely contain little useful information, thus there is a risk of overfitting; however, a further benefit to using the AICC instead of another GoF metric is that it reduces the instances of overfitting by limiting the number of fragments included in the regression. Fragment Selection Rules. A fragment selection rule has been defined to prevent combinations of fragments from being included together during model selection that have the potential to lead to instability and collinearity. In general the fragment selection rule excludes instances of fragment coincidence; meaning fragments which are not perfectly collinear but which are present in many or all of the same chemicals in the training data set. Coincident fragments create the possibility of additively combining fragment counts to create perfect collinearity among three or more fragments. A pair of fragments is defined as coincident if either of them occurs independently in fewer than three chemicals. Many fragments are coincident with basic fragments and therefore applying selection rules to basic fragments would prevent most other fragments from being included in the model selection, so the fragment selection rule is not applied when one or both of a pair of fragments compared is a basic fragment. An additional selection rule is applied when considering rare fragments. The standard error of the MLR coefficient for a rare fragment must be less than half of the absolute magnitude of the coefficient, meaning that p < 0.05. This must be true when the fragment is first added and must still be true each time fragments are tested for removal from the MLR. If the standard error becomes larger than half the absolute value of the coefficient, indicating p > 0.05, the rare fragment will be removed regardless of how this affects the GoF. Data Splitting, Model Validation, and Cross-Validation. After the initial fragmentation of the full data set and identification of basic fragments, the data set is split into training and validation data sets. The training data set is seeded with the two chemicals with the highest and lowest values for the property to be predicted, and the validation data set is seeded with the chemical that has the highest property and structural similarity to these two chemicals. Allocation of the remaining chemicals into the training and validation data sets proceeds by adding chemicals alternatively to each data set. In each case the chemical with the lowest property and structural similarity to the chemicals already in the data set is selected. More information on the definition of property and structural similarity is provided in the SI section SI2. During model selection k-fold cross validation is applied. The training data set is split into multiple internal training and validation data sets, referred to hereafter as fitting and testing data sets (internal validation) in order to differentiate them from the overall training and validation data sets (external validation). The result is k fitting data sets each composed of k − 1 folds with k corresponding testing data sets each composed of the fold missing from its fitting data set. A separate model is created by MLR on each fitting data set using the same set of fragments, and the final model presented is the average of the models developed on the individual fitting data sets. Multiple models are created simultaneously; therefore, an overall AICC value is calculated to simultaneously test the GoF for the full training data set. The standard AICC formula30 is applied but instead of using the residual sum of squares of each fitting data set the predictive sum of squares of all of the testing data sets is used, so that square of the residual of every chemical
in the training data set is counted once as it is predicted by the corresponding fitting data set. This is analogous to using the leave many out Q2 cross validation GoF metric but has the additional advantages of using the AICC. Division of the full training data set into the k folds was done using the same algorithm used to divide the data set into training and validation sets, but with a slight modification that causes each individual k fold to be composed of a cluster of very similar chemicals. More details are provided in the SI section SI2. Domain of Applicability. A domain of applicability is defined for developed QSARs based on structural similarity to chemicals in the training data set and the residuals of the fitted values for the training data set. Structural similarity is defined as the degree of fragment overlap between two chemicals, considering only those fragments that are in the developed QSAR. For any external chemical the top five most structurally similar training data set chemicals are taken along with the residuals of their fitted values and combined into a chemical similarity score (CSS) with a value between 0 and 1. The CSS is calculated for all chemicals in the validation and training data sets, and a cutoff value is defined such that 95% of the training chemicals fall within the domain of applicability of the developed QSAR. Further details regarding the definition of the domain of applicability are provided in the SI section SI3. Application of the IFS Method to Experimental HL Data. The IFS model selection algorithm was applied to a database of fish biotransformation half-lives (log10 HL N according to previous nomenclature and methods) for 632 organic chemicals7,20,21 to develop and test two new HLQSARs. The first new HL-QSAR (referred to as IFS-EPI) uses the training and validation data sets from the original publication.7,20 The second new HL-QSAR (referred to as IFS-HLN) uses training and validation data sets created by the method based on property and structural similarity described above and in the SI. Two-dimensional chemical structures (SMILES notation31) were retrieved from the EPI Suite database7 and manually checked for correctness. A number of chemicals in the HLN database differ only by their stereochemistry, the effects of which are not captured by the 2D fragments generated in this work. Therefore, for each unique 2D structure the log HLN values were averaged and treated as a single data point. This reduced the total number of chemical structures from 632 to 619. The original training data set used in the BCFBAF (HL-QSAR) model in EPI Suite,7,20 hereafter referred to as the EPI-HLN model, contained 421 chemicals but this has been reduced to 417 for the IFS-EPI data set by averaging the expected values for chemicals with identical 2D structures. The number of chemicals in the testing data set likewise was reduced from 211 to 209 chemicals. The training and testing data sets for the IFS-HLN QSAR contain fewer total chemicals (n = 412 and n = 207, respectively) because averaging of identical 2D structures was done before splitting the data set into training and testing data sets instead of afterward. The distribution of chemicals into the various training and validation data sets is summarized in SI Table SI1.
■
RESULTS Fitting and Validation Statistics. Cross validation was performed with values of k ranging from 3 to 10, and the resulting QSAR with the best statistics for the validation data set was selected and presented here. Selection was primarily based on the concordance correlation coefficient (CCC)32 of 8255
dx.doi.org/10.1021/es301182a | Environ. Sci. Technol. 2012, 46, 8253−8260
Environmental Science & Technology
Article
the prediction for the validation data set, but in each case the selected QSAR also had among the highest R2, and lowest root mean squared error (RMSE). For IFS-HLN the QSAR with k = 8 was selected, and for IFS-EPI the QSAR with k = 6 was selected. Statistics measuring the GoF for the validation data sets of the IFS-HLN and IFS-EPI data sets showed no obvious trend with increasing values of k. Table 1 outlines the statistics of the training and validation data sets for the two selected QSARs (IFS-HLN, IFS-EPI) and Table 1. Regression and Prediction Statistics of the Three Selected QSARsa IFS-HLN k fragments
8 36
n AICC R2 RMSE slope
412 −146.0 0.789 0.526 0.791
n CCC R2 RMSE slope
207 0.857 0.748 0.584 0.758
IFS-EPI 6 38 training data set 417 −150.6 0.804 0.517 0.797 validation data set 209 0.846 0.717 0.622 0.817
EPI-HLN 59 421 −115.2 0.821 0.494 0.821 211 0.856 0.734 0.600 0.819
a
k = number of cross validation folds; AICC = Akaike information criterion corrected for data set size; CCC = concordance correlation coefficient; R2 = coefficient of determination; RMSE = root mean squared error.
those of EPI-HLN for comparison, and Figure 1 shows the fits for the training and validation data sets of IFS-HLN. IFS-HLN and IFS-EPI have very similar fitting statistics for their respective training data sets, but for the validation data sets IFS-HLN has superior validation statistics (CCC, R2, and RMSE, see Table 1). EPI-HLN has a slightly better fit to its training data set than IFS-HLN or IFS-EPI but the predictive power for their respective validation data sets is slightly better (IFS-HLN) and slightly worse (IFS-EPI). The number of chemicals with absolute residual values greater than one log unit is 28 in the training data set and 20 in the validation data set for the IFS-HLN model, 30 and 21 for the IFS-EPI model, and 21 and 22 for the EPI-HLN model. A y-scramble was performed on the IFS-HLN QSAR to test the possibility of chance correlations. The experimental values of the training data set were randomized and the QSAR was refitted and cross-validated leaving all other parameters unaltered. Randomization and refitting were performed 50 times and the fitting statistics were compared to the IFS-HLN QSAR. The average RMSE of the y-scramble for the training and validation data sets were 1.104 (minimum 1.085) and 1.204 (minimum 1.092), respectively. The average AICC of the training data set was 129.1 (minimum 123.0) and the average CCC of the validation data set was 0.013 (maximum 0.190). Fitting statistics for the y-scramble are very poor, indicating that very little of the statistical goodness of fit for the IFS-HLN QSAR is due to chance correlations. Descriptors Selected. Both IFS-HLN and IFS-EPI have a greater than 1:10 descriptor to training data-point ratio, and all of the other k-fold cross-validations created for this work are in
Figure 1. IFS-HLN fits for training data set (A) and validation data set (B). Purple, red, yellow, light green, and dark green circles represent chemicals with CSS scores better than
