Abstract

Assessing the ecotoxicity of pharmaceuticals is of urgent need due to the recognition of their possible adverse effects on nontarget organisms in the aquatic environment. The reality of ecotoxicity data scarcity promotes the development and application of quantitative structure activity relationship (QSAR) models. In the present study, we aimed to clarify whether a QSAR model of ecotoxicity specifically for pharmaceuticals is needed considering that pharmaceuticals are a class of chemicals with complex structures, multiple functional groups, and reactive properties. To this end, we conducted a performance comparison of two previously developed and validated QSAR models specifically for pharmaceuticals with the commonly used narcosis toxicity prediction model, i.e., Ecological Structure Activity Relationship (ECOSAR), using a subset of pharmaceuticals produced in China that had not been included in the training datasets of QSAR models under consideration. A variety of statistical measures demonstrated that the pharmaceutical specific model outperformed ECOSAR, indicating the necessity of developing a specific QSAR model of ecotoxicity for the active pharmaceutical contaminants. ECOSAR, which was generally used to predict the baseline or the minimum toxicity of a compound, generally underestimated the ecotoxicity of the analyzed pharmaceuticals. This could possibly be because some pharmaceuticals can react through specific modes of action. Nonetheless, it should be noted that 95% prediction intervals spread over approximately four orders of magnitude for both tested QSAR models specifically for pharmaceuticals.

1. Introduction

It is evident that humans benefit a lot from pharmaceuticals, a class of active chemicals that can interact with biological receptors in humans. However, also because of their reactive property [1], the dark side of the coin is that they could pose potential ecological threats to nontarget organisms after use and release into the environment [2, 3] because nontarget species like fishes could share similar biological receptors to drugs to humans [48]. Demonstrating and assessing the ecotoxicological relevance of the presence of pharmaceuticals as emerging contaminants in the environment is of increasing concern during recent years [2, 913].

To help understand the potential risk of pharmaceuticals to nontarget organisms, the first step is to determine their occurrence levels in various environmental matrices. Owing to the development of analytical technology for the detection of extremely low concentrations of pharmaceuticals, intensive monitoring has been conducted worldwide especially during recent decades [1419]. While there is still a need for extensive monitoring practices considering the various temporal and spatial scales [14], it is a more urgent demand of a high amount of measured ecotoxicological data that are equivalently salient to assess the environmental risk of the identified and potential pharmaceutical contaminants.

Unfortunately, these kinds of ecotoxicity data for pharmaceuticals are highly scarce at present because the tests for assessing chemical ecotoxicity are usually costly, labor-intensive, and time-consuming. When the measured data are not available, Quantitative Structure Activity Relationship (QSAR) models could offer an effective alternative [2022]. The QSAR model generally describes the relationship between structural features and a toxicity endpoint, which is developed and validated by a rigorous procedure using the existing measured data [23, 24]. With a validated QSAR model, the toxicity of a chemical without measured data can be predicted. QSAR models have been extensively employed to estimate the hazard of pharmaceuticals, especially under the context of ranking and priority screening of highly risky candidates for further study and monitoring [2530].

In the existing applications of QSAR models, the Ecological Structure Activity Relationships (ECOSAR) model was among the most commonly used one for estimating the ecotoxicity of pharmaceuticals [10]. ECOSAR is a well-established modelling tool based on the octanol-water partition coefficient (log Kow) by United States Environmental Protection Agency. This log Kow-based QSAR model worked well for predicting nonspecific narcosis toxicity that is a function of the tendency to dissolve into membranes. However, it could be a problem for pharmaceuticals. On the one hand, pharmaceuticals are complex chemicals with multiple functional groups, which will lead to uncertainty when designing an appropriate chemical class [21]. On the other hand, pharmaceuticals are always thought to cover a wide range of toxic mechanisms. Hence, specific modes of action (MOA) could be expected when they interact with nontarget organisms [10, 31], possibly making the narcotic toxicity prediction model inaccurate. Madden et al. [32] pointed out that toxicity predictions by ECOSAR for pharmaceuticals should be explained with caution, and the applicability domain (AD) should be evaluated carefully, considering the fact that ECOSAR was developed based on a dataset consisting of mainly industrial chemicals of simple structure with only a single functional group.

Recently, some authors argued that it is of necessity to develop the QSAR model of ecotoxicity specifically for pharmaceuticals as a substitute to ECOSAR, considering their special properties. As a result, some of QSAR models have been reported in literature [3336]. In particular, Tugcu et al. [34] developed a QSAR model (hereafter called the “Tugcu’s model”) specifically for pharmaceuticals to estimate the acute toxicity to fish using multiple molecular descriptors. They concluded that their model could provide better predictions of toxicity than ECOSAR. In another study conducted by Sangion and Gramatica [33], QSAR models (hereafter called the “Sangion’s model”) of acute toxicity to three trophic levels of aquatic organisms were developed. The comparison of their model results with predictions obtained by using ECOSAR showed that ECOSAR was also outperformed. Nonetheless, their comparisons were limited to a small subset of pharmaceuticals, and more work should be done to further verify whether the better performance of QSAR models specifically for pharmaceuticals is solid.

The aim of the present study is further to clarify whether a specific QSAR model is indispensable for the accurate prediction of the ecotoxicity of pharmaceuticals. To this end, we predicted the ecotoxicity of a totally new set of pharmaceuticals for which the measured data are available, by using two externally validated QSAR models (i.e., Tugcu’s model and Sangion’s model) and ECOSAR. And then, the predictive ability of each pharmaceutical specific QSAR model was compared with that of ECOSAR, respectively, to demonstrate the goal of the present study. Possible factors affecting the performance of QSAR models were also discussed.

2. Methods

2.1. Description of QSAR Models Specifically for Pharmaceuticals

Two validated QSAR models, specifically for pharmaceuticals (i.e., Sangion’s model and Tugcu’s model), were selected for the comparison of predictive ability with ECOSAR. Sangion’s model was integrated in a well-defined and user-friendly software package, QSARINS (v. 2.2.2) [37, 38]. It was developed for toxicity prediction of pharmaceuticals in three different endpoints, and species were selected: growth rate inhibition at 72 h (EC50) in Pseudokirchneriella subcapitata, immobilization at 48 h (EC50) in Daphnia magna, and mortality at 96 h (LC50) in Pimephales promelas. Tugcu’s model was developed to predict the lowest 96-hour mortality (LC50) in various freshwater fish species. Details of the description for both models can be found in Sangion and Gramatica [33] and Tugcu et al. [34].

2.2. Collection of Measured Ecotoxicity Data

An initial list containing 663 compounds was compiled that consists of pharmaceuticals produced during 2011–2015 in China [3943], by excluding vitamins, electrolyte solution additives, biochemical pharmaceuticals, as pharmaceuticals with a molecular weight over 2000 Da. The organometallics and inorganic compounds were also excluded in the initial list. Basic information was collected, including the Chemical Abstract Service Register Number (CASRN), molecular weight, molecular form, and therapeutic group. The CASRN and name of pharmaceuticals were cross validated by searching in multiple sources, including PubChem (https://pubchem.ncbi.nlm.nih.gov/) and SciFinder (https://scifinder.cas.org/), to assure the accurate citation of the CASRN. We searched measured ecotoxicity data in the ECOTOX database (https://cfpub.epa.gov/ecotox/) for the 663 pharmaceuticals in terms of CASRN, endpoints, and species following the purpose of this study. For the Tugcu’s model, measured data of LC50 for various freshwater fish species were collected. In case of multiple measured toxicity values for the same specific endpoint, the lowest value was selected to maintain consistency with the data selection criteria used in model development.

To make an independent evaluation, pharmaceuticals that have measured toxicity data but had not been included in the training set of the tested QSAR models were selected. For comparison of predictive ability of Sangion’s model with ECOSAR, a total of 41 toxicity data were identified, including 10 EC50s of Pseudokirchneriella subcapitata, 27 EC50s of Daphnia magna, and 4 LC50s of Pimephales promelas. Regarding the case of Tugcu’s model versus ECOSAR, we identified 43 LC50s of various fish species that can be used for comparison. All the data were compiled in Tables S1 and S2. The negative logarithm of the measured LC50 or EC50 values (pT) was calculated on a molar basis and used for comparison with the predicted ones.

2.3. Generation of Predicted Ecotoxicity Data

Regarding Sangion’s model, the structures of pharmaceuticals with measured data were sketched using Spartan 16 (v 2.0.7) and geometrically optimized by the semiempirical AM1 method. The predicted ecotoxicity data were generated by using QSARINS, in which the theoretical molecular descriptors were calculated by the PaDEL-Descriptor software (v 2.21). The leverage approach was used to check whether a given compound is within the scope of AD on the basis of the molecular descriptor space defined according to the training set during model development. The leverage value was automatically calculated by the QSARINS.

In the case of Tugcu’s model, the structures of pharmaceuticals were also drawn using the Spartan 16 software. The conformation was optimized at the lowest energy by employing the semiempirical PM3 method. The molecular descriptors were calculated using the Dragon software (v 7.0.10), the same one as that by Tugcu et al. [34]. The predicted pTs were calculated using the established linear model with the reported molecular descriptors. The leverage approach was used to determine if a tested pharmaceutical is within the scope of AD of the developed model.

ECOSAR incorporated in EPI Suite (v 4.11) was run for the narcotic toxicity prediction by using CASRN or SMILES as inputs. As pharmaceuticals are not a traditional class of structurally or functionally similar compounds, ECOSAR would allocate some pharmaceuticals to more than one class, and consequently, there is more than one result for each of them. We took baseline toxicity predictions that handled pharmaceuticals as neutral organic, considering the goal of the present study. The tested pharmaceuticals were manually checked if they are inside the AD by checking log Kow values. The results showed that all the tested pharmaceuticals were in the AD of ECOSAR.

2.4. Measure of Predictive Ability

To demonstrate whether the predictive ability of the tested QSAR models, a series of summary statistics were calculated for comparison: Pearson correlation coefficient, differences between predicted and measured pT (referred as minimum, mean, and maximum absolute residual, root of mean squared residual (RMSE), and calculated as equations (1)–(4)), and the percentage of pharmaceuticals with differences between predicted and measured pT was less than factors of 2, 5, 10, 100, and 1000. A linear regression analysis of predicted and measured pT was conducted, and the 95% prediction intervals were also calculated. The statistical analysis was accomplished in Microsoft Excel 2016.where pTi and are the measured and predicted values for negative logarithm of the LC50 or EC50.

3. Results

3.1. Sangion versus ECOSAR Model

Statistics characterizing the predictive ability of Sangion’s and ECOSAR models are summarized for comparison in Table 1. It is evident that the pharmaceutical specific model, Sangion’s model, was more able to predict the toxicity trend than the generalized narcotic toxicity prediction model. The predictive ability of Sangion’s model significantly surpassed ECOSAR no matter which parameter was used as a metric. Sangion’s model, which is specific for pharmaceuticals, can give a prediction with uncertainty lower than two orders of magnitude for over 80% of the tested pharmaceuticals, outperforming that of ECOSAR (approximately 60%). If considering only pharmaceuticals in AD, this variable was further improved to around 95%.

For predictions produced by ECOSAR, high RMSE (2.7 log unit, Table 1) was observed in the present study, which is in accordance with [33]. However, our application of Sangion’s model to a new set of pharmaceuticals showed an RMSE of 1.4 log unit (Table 1) that is significantly higher than that (0.6 log unit) applied to the prediction set during model development. This could be attributed to a few influence factors that will be discussed later.

The better performance of Sangion’s model is also clearly demonstrated in Figure 1. Narrower 95% prediction intervals were found for Sangion’s model, while the 95% prediction intervals were approximately eight orders of magnitude wide for ECOSAR. Furthermore, the regression analysis showed that Sangion’s model in general underestimated the toxicity of the most toxic pharmaceuticals (Figure 1(a)), while the baseline toxicity model tended to underestimate the toxicity of most pharmaceuticals (Figure 1(b)).

3.2. Tugcu versus ECOSAR Model

Similarly, the results showed that Tugcu’s model outperformed ECOSAR, although a slightly higher Pearson correlation coefficient was obtained for ECOSAR than Tugcu’s model (Table 2). Nonetheless, rest of the statistics clearly indicated the superiority of Tugcu’s model to ECOSAR. In this application, only three pharmaceuticals that have measured toxicity data were out of the AD of Tugcu’s model. No obvious improvement in model performance was observed by restricting the dataset to include only pharmaceuticals in the AD.

Using the mean absolute residue (MAR) as a metric, Tugcu et al. [34] found in their application that their model and ECOSAR had MAR of 0.348 and 0.872, respectively, indicating that average difference between predicted and measured toxicity was less than one order of magnitude for both models. In contrast, in the present application, the pharmaceutical specific model had a MAR of 1.14, indicating that the mean difference was slightly more than one order of magnitude (Table 2).

The better performance of Tugcu’s model over ECOSAR was also demonstrated by the significantly narrower 95% prediction intervals as shown in Figure 2. Again, the 95% prediction intervals for predictions of ECOSAR were spread nearly eight orders of magnitude. It is also found from Figure 2 that Tugcu’s model underestimated the toxicity for the most toxic compounds but overestimated the toxicity for the least toxic compounds. And, the underestimation still existed for most pharmaceuticals when using ECOSAR to predict fish toxicity.

It should be noted that Tugcu’s model was developed using ecotoxicity data of all fish species, which could influence the accuracy of model predictions. And, this was demonstrated to some extent when we did an intercomparison between the two pharmaceutical specific models. It is obvious that the predictive ability of Sangion’s model was little better than that of Tugcu’s model (Tables 1 and 2) although the intercomparison is not totally reasonable as they are applied to different pharmaceuticals.

4. Discussion

Our results clearly demonstrated that a better performance was observed for either Sangion’s or Tugcu’s model and the QSAR model specifically for pharmaceuticals, compared to the generalized baseline toxicity prediction model of ECOSAR. The present study further adds evidence that we need a QSAR model of ecotoxicity specifically for pharmaceuticals.

Not surprisingly, the model with more descriptors and sophisticated statistical methodology had a better overall performance. ECOSAR have been proved to be a valid tool for the prediction of ecotoxicity of the relatively simpler chemicals, but it could be problematic when being applied to predict the toxicity of chemicals with complex structures and functional groups as those of pharmaceuticals. In other words, the use of too simple equations cannot deal with the complexity of the phenomenon. Assigning a pharmaceutical, which is chemically complex, to one chemical class, means that only one functional subgroup is considered as responsible of the activity, and this can even be false. Modern QSAR models developed and trained on multiclass chemicals are by definition better able to work with more functional groups.

In the present study, the toxicity was generally underestimated by ECOSAR (Figures 1(b) and 2(b)). To further illustrate the predictive ability of ECOSAR, we also considered the lowest predictions in ECOSAR (worst case) in both comparisons. While the average difference between predicted and measured toxicity decreased by approximately 0.3 log unit, the ECOSAR still underestimated the toxicity for nearly half of the tested pharmaceuticals (Figures S1 and S2 in the Supplementary Materials).

The underestimation could be due to that the baseline toxicity, also known as the minimum toxicity, could not reflect the real ecotoxicity, and a specific MOA could exist for those reactive pharmaceuticals [44, 45]. Some authors pointed out that a specific MOA for pharmaceuticals could be indicated by the poor correlation between the acute toxicity and log Kow [10, 31]. However, others have argued that log Kow is a good predictor for most pharmaceuticals, indicating that they are the narcotic [46]. Sanderson and Thomsen [47] further concluded that a nonspecific MOA was suggested by the acute to chronic ratio (ACR) lower than 25 for nearly 70% of analyzed pharmaceuticals. However, it is complicated to recognize the true MOA of pharmaceuticals at present, especially for chronic MOA. The accuracy of using ACR to help determine the MOA of a compound heavily relied on the accurate determination of chronic toxicity. However, lacking high quality chronic data is indisputable, and therefore, high ACRs may also indicate a narcosis and conversely low ACRs may not necessarily indicate a nonspecific MOA [48]. Moreover, the standardized tests are designed for tests of mortality, which may overlook the sublethal effects even if there is any [46]. We have derived the threshold for two pharmaceuticals in different therapeutic groups and found that the traditional endpoints related to population mortality could not be the most sensitive one [49]. Adverse effects may occur via a possibly specific MOA that an organism would die from for a long run. It is found that, over two years, chronic exposure to 17α-ethylestradiol at extremely low concentrations (about 5 ng/L) led to a collapse of fathead minnow populations in a lake [44]. Further complicating this is that even for the same compound, MOAs may vary between organisms, at different life stages and at different exposure concentrations [50], all of which would lead to uncertainties in ACRs’ calculation.

Although the pharmaceutical specific QSAR models outperformed that of ECOASR, it should be noted that the 95% prediction intervals were approximately four orders of magnitude wide for both Sangion’s and Tugcu’s models (Figures 1 and 2).

Moore et al. [21] evaluated the performance of six QSAR packages that predict acute toxicity to fish for industrial chemicals and found even the best package which produced variation of four orders of magnitude, as our assessment for pharmaceuticals’ toxicity prediction by the specific QSAR models. It is unclear if this is a systematic error that is caused by the existing uncertainty in measured toxicity data. We made a perfect assumption for the measured toxicity values that they are accurate, precise, and unbiased reflection of the true toxicity of pharmaceuticals to the aquatic organisms. Nonetheless, this is impossible considering that the measured toxicity data were determined under different testing protocols. Therefore, large residuals can be partly attributed to uncertainty that inherently exists in measured data.

On the other hand, it should also be kept in mind that no perfect model has ever been established. We examined the outlier substances that were defined as those substances that had differences between model predictions and measured values greater than 100 folds. Regarding Sangion’s model, two outlier substances were identified: doramectin and niclosamide. They are both anticestodals that were not included in the training set of Sangion’s model. Seven out forty-one substances were identified as outliers for Tugcu’s model. The outlier substances include the antineoplastic drug, the estrogen, and the anticestodal drug. Again, they were rarely included in the training set during model development. And, estrogen is a well-proven contaminant that has specific MOA [51]. Therefore, limited coverage of pharmaceuticals of different therapeutic groups could be one of potential aspects influencing the performance of the pharmaceutical specific model.

It should be noted that the outlier substances were in the AD of the model, which means that the present molecular descriptors integrated in the model could not totally explain the toxicity variances for pharmaceuticals of different therapeutic uses. The specific QSAR models are “average” statistical results covering limited kinds of pharmaceuticals. They are still MOA-mixed QSAR models considering the fact that the pharmaceuticals of different therapeutic uses could possess different MOAs. It is hard to foresee the performance of a future developed model covering a wider range of pharmaceuticals belonging to various therapeutic groups. Either it will cover more predictors or will it be split to different models based on MOA to improve the model performance. The MOA-based QSAR has been advocated for a long time [52, 53]. However, there are great challenges to thoroughly understand the MOA of pharmaceutical’s toxicity to nontarget organisms. As for pharmaceuticals, it has been argued that the mammalian pharmacology data could be potentially used to help understand their MOA [54]. The MOA-based QSAR model specifically for pharmaceuticals would be one of the approaches for an improved ecotoxicity prediction of these chemicals.

5. Conclusions

QSAR models have been recognized the most common tools for complementing the scarcity of toxicity data in scientific research and management. Our results showed that the generally used baseline toxicity prediction model, ECOSAR, usually underestimated the toxicity of pharmaceuticals in the present practice. The better performance of both Sangion’s and Tugcu’s models clearly demonstrated the necessity for QSAR models specifically for pharmaceuticals with complex structures and multiple functional groups. Modern QSAR models that are able to work with more functional groups are urgently needed. However, there is also a need for improving the performance of the present pharmaceutical-specific QSAR models either by covering a large pool of pharmaceuticals in different therapeutic groups or finally by developing an MOA-based model with the aid of existing pharmacology data.

Data Availability

The numerical data used to support the findings of this study are included within the article and the associated supporting information.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was cosupported by the Beijing Natural Science Foundation (8162037) and the Fundamental Research Funds for the Central Universities in China (2021YJSHH12 and C202003288). Qingwei Bu was funded by Yue Qi Young Scholar Project, China University of Mining & Technology, Beijing (2017QN15). The authors also appreciate Prof. Paola Gramatica, University of Insubria, for providing the free license of the QSARINS software package.

Supplementary Materials

Table S1: the measured (the lowest observed value) and predicted toxicity data for pharmaceuticals used for performance comparison of Sangion’s model versus ECOSAR. Table S2: the measured (the lowest observed value) and predicted toxicity data for pharmaceuticals used for performance comparison of Tugcu’s model versus ECOSAR. Figure S1: measured toxicity versus the lowest model prediction of ECOSAR for pharmaceuticals used for performance comparison of Sangion’s model versus ECOSAR. The best-fit linear regression and 95% prediction intervals are shown. Figure S2: measured toxicity versus the lowest model prediction of ECOSAR for pharmaceuticals used for performance comparison of Tugcu’s model versus ECOSAR. The best-fit linear regression and 95% prediction intervals are shown. (Supplementary Materials)