Integral estimation of xenobiotics’ toxicity with regard to their metabolism in human organism

Drug metabolism prediction

To predict the particular metabolites and generate the complete metabolic tree, we used the computational method developed earlier, which is based on the prediction of possible classes of nine biotransformation reactions and the subsequent estimation of metabolic sites (SOM) for each reaction. We considered nine classes of biotransformation reactions: aliphatic and aromatic hydroxylation, N- and O-glucuronidation, N-, S- and C-oxidation, and N- and O-dealkylation. These reactions are catalyzed by five major isoforms of human cytochromes P450s (1A2, 2C19, 2C9, 2D6, 3A4) and by human UDP-glucuronosyltransferase. The training set was created using information about drug-like compounds’ biotransformations from the Biovia Metabolite database [22]. The procedure of metabolite generation includes three steps: (1) prediction of the class of reaction for the entire molecule; (2) prediction of SOM for each reaction class; (3) calculation of the probability of metabolite formation. The probability of metabolites formation is calculated using the Bayesian approach based on the analyses of “structure-biotransformation reactions” and “structure-modified atoms” relationships. By a stepwise application of this procedure, we obtained a set of metabolites generated at each particular stage of biotransformation with the estimated probabilities of their formation. This method is described in detail in our publications [15], [16], [17], [20].

Rat acute toxicity prediction

We have earlier proposed the method for rodents’ acute toxicity prediction and considered its applicability for prediction of rat acute toxicity under four routes of administration: oral, intravenous, intraperitoneal and subcutaneous [19]. The method involves the electrotopological quantitative neighbourhoods of atoms (QNA) [18], atom centric substructural multilevel neighbourhoods of atoms (MNA) descriptors, topological length and volume of molecules for representation of structures and self-consistent regression for establishing of QSAR [23]. As shown earlier, the method proposed provides the reasonable accuracy of prediction and better performance comparing to the US EPA T.E.S.T. program [19]. This method was implemented in the GUSAR program [23] used in this study. For a detailed description of this method, see [19], [23]. The on-line version of web-service for prediction of LD₅₀ values of rat acute toxicity with four routes of administration is freely available via the web interface [24].

In this work, we consider only one intravenous route of administration because in this case the immediate toxic reaction induced by the interaction of a drug-like compound with the organism is observed. Furthermore, the proper statistical parameters regarding both the R2 and coverage values had been previously found for the intravenous route of administration [19].

To extend the applicability domain of QSAR model for prediction of rat acute intravenous toxicity, we significantly increased the training set using data from the BIOVIA Toxicity Database [25]. Currently, it includes the information about 3632 chemical structures with the data on 3965 LD₅₀ values of rat acute intravenous toxicity. We systematize the data in accordance with the current “good QSAR modeling practice” [26], [27] by removing salts, mixtures, inorganic compounds, and polymers. The median LD₅₀ value was calculated for those structures, which had more than one experimental LD₅₀ value. To provide the comparability of toxicity data, we transformed the toxicity end-points values from the commonly applied in toxicology mg/kg units to the log10(1/LD₅₀) (mol/g) representation of LD₅₀ values in mol/g. The final training set included 2978 structures.

GUSAR created 320 QSAR models varying parameters of using QNA and MNA descriptors. The complete description of models is presented in supplementary materials. All these models were used as a single averaged consensus model. Table 1 shows the statistical parameters of the model.

Both Table 1 and the plot in Fig. 1 with the observed vs. predicted LD₅₀ values for the training set (Fig. 1) demonstrate a reasonable range of the toxicity end-point values and reasonable “goodness-of-fit” characteristics.

Fig. 1:

The observed vs. predicted log10(1/LD₅₀) (mol/g) values for the training set of rat acute intravenous toxicity.

Test set

To validate the quality of the QSAR model developed for the rat acute intravenous toxicity and to compare the performance of four different numerical parameters of the integral toxicity estimation, we created the external test set included 37 drug structures from the 22 ChEMBL database release [28]. The experimental data on metabolic reactions of these drugs, structures of about 200 of their metabolites and the LD₅₀ values for the parent compounds (rats, intravenous route of administration), were extracted from the BIOVIA Toxicity Database [25]. The list of these drugs and their LD₅₀ values calculated as a median lethal dose in the N independent studies in rats when administered intravenously are given in Table 2.

Methods proposed for evaluation of integral toxicity

The definitive goal of our studies is to improve the general quality of the toxicity risk assessment of drug-like compounds taking into account the complexity of the effects produced in an organism by a particular drug and all its metabolites.

We assume that the doses of metabolites (D_i) may be represented from the dose of the initial compound (D) and probabilities (P_i) of its metabolite formation as follows:

D i ≈ D ∗ P i ,  thus:  ∑ i D i ≈ D

Based on the probability of metabolite formation and their acute toxicity values calculated by our method [20], we proposed four different methods for calculating the integral toxicity.

Method 1. The integral toxicity is calculated as the effect of toxicity of all metabolites without taking into account the toxicity of the parent compound. It is assumed that each of the parent compounds is metabolized and all its metabolites are formed with an equal probability P_i=1/N (N is an amount of metabolites). We estimate an integral LD₅₀ value LD_M1 using the LD₅₀ values predicted for metabolites (LDi) as follows:

L D M 1 = N ∑ i 1 L D i

Method 2. The integral toxicity is calculated as the toxic effects of the parent compounds and all its metabolites. We assumed that the parent compound is not completely metabolized and may be considered together with its metabolites. We estimate an integral LD₅₀ value LD_M2 using the LD₅₀ values predicted for metabolites (LD_i) and for the parent compound (LD₀) as follows:

L D M 2 = N + 1 ∑ i 1 L D i + 1 L D 0

Method 3. The integral toxicity is calculated as the toxic effects of the most toxic metabolite. We assumed that the effect of the most toxic metabolite is crucial, even if it is formed in a small quantity. We estimate an integral LD₅₀ value LD_M3 as follows:

L D M 3 = L D most toxic

Method 4. The integral toxicity is calculated as the toxic effects of the most toxic metabolite and the parent compound. We assumed that both effects of the most toxic metabolite and the parent compound are important. We estimate an integral LD₅₀ value LD_M4 as follows:

L D M 4 = 2 1 L D most toxic + 1 L D 0

Thus, for each of the parent compounds from the test set, we calculated all the four values mentioned above and compared them with the experimental LD₅₀ values (rats, i.v.).

Results of integral toxicity estimation for the test set

Predicted values of rat acute intravenous toxicity for the test set estimated using four different methods are presented in Table 2.

To compare the performance of the integral toxicity estimation, we calculated the R2 and RMSE values characterizing the conformity between the experimental and predicted data presented in log(mol/g) (Table 3). From the data presented in Table 3, we can conclude that the values of the R2 test are equal or exceed 0.759 and the RMSE test values are <0.487 for all four methods.

It is clear from the data presented in Table 3 that the relative performance of different methods arranged in the ascending order, the R2 values are as follows: M1<M2<Parent compound<M3<M4. In terms of RMSE the order is slightly different: M2<M4<M1<M3<Parent compound. Thus, we can conclude that the utilization of the M4 method (the effect of the most toxic metabolite and the parent compound on integral toxicity) provides the best performance. The accuracy obtained is close to the accuracy of estimation that takes into account the toxicity of the parent compound. It may be explained by the fact that the measured experimental values of LD₅₀ (i.v. rats) of the training set already include the hidden values of metabolite toxicity and the most toxic metabolites play a major role.

The results show that the consideration of the LD₅₀ values predicted for both the parent compound and the most toxic metabolite (method 4) provides the better estimates of integral toxicity. The prediction results obtained by method 1 (consideration of the predicted acute toxicity for all metabolites without taking into account those of the parent compound) are less accurate. It allows concluding that the role of the parent compound in integral toxicity is significant in most cases. The values of toxicity obtained by methods 1 and 2 are lower in comparison with those obtained by the other three methods for the test set. It may be explained by the commonly known fact that metabolism usually leads to a decrease in xenobiotics toxicity. An opposite situation is observed in the case of method 3, when the toxicity of the most toxic metabolite is used as the estimate of integral toxicity.

The example of the paroxetine metabolic pathway (see Fig. 2) and the results obtained (see Table 4) for antidepressant drug paroxetine are presented below.

Fig. 2:

Paroxetine metabolism pathway [29].

As can be seen from the data presented in Table 4, the predicted acute toxicity for paroxetine significantly exceeds the experimental LD₅₀ value. Consideration of the paroxetine biotransformation leads to a decrease in the predicted values in all cases; the best correspondence between the predicted and experimental values is observed for method 3 that estimates the integral toxicity on the basis of the values predicted for the most toxic metabolites. The integral toxicity calculated by method 4 that takes into consideration not only the toxicity of the most toxic metabolite but also the predicted toxicity of the parent compound does not show the best performance due to the low toxicity predicted for the parent compound.

A possibility of the practical usefulness of the integral toxicity prediction is demonstrated in the case study of the antiepileptic drug felbamate mentioned above. Prediction of rat acute toxicity for felbamate (Table 5) makes it possible to suggest that this compound is low toxic. However, the felbamate biotransformation pathway presented in Fig. 3 provides fewer values of felbamate integral toxicity predicted by both methods 3 and 4 (Table 5). According to these estimates, the compound belongs to a more toxic third class of acute intravenous toxicity [31], which corresponds to the experimental results obtained in human but was not observed in rodents.

Fig. 3:

Felbamate metabolism pathway [30].

Prediction results of the integral toxicity obtained for felbamate directly correspond to the cause of its withdrawal [30], [32]. Felbamate was approved in 1993 as an antiepileptic drug. However, it was soon withdrawn due to the cases of aplastic anemia, and hepatic failure resulted in deaths. The manifestation of toxic effects was associated with metabolite 2-phenylpropenal as a reactive electrophile responsible for the toxicities of felbamate observed.