A modeling approach to evaluate the balance between bioactivation and detoxification of MeIQx in human hepatocytes

Compound concentrations in human hepatocytes in response to MeIQx initial concentrations

The data set used for this present work has been previously published (Langouët et al., 2001). Briefly, human hepatocytes isolated from human liver biopsies were treated with low concentration of MeIQx (1, 10 and 50 µM) not to induce cytotoxic effects (Langouët et al., 2001). The culture media were retrieved at 24 h and the metabolites were purified and quantified by mass spectrometry. Figure 2 shows the percentage of total metabolite quantity for the ten following compounds: MeIQx, MeIQx-N²-Gl, MeIQx-N²-SO₃H, N-desmethyl-7-oxo-MeIQx, 7-oxo-MeIQx, IQx-8-COOH, 8-CH₂-OH-IQx, HON-MeIQx-N²-Gl, HONH-MeIQx, and Potential-Genotoxic-Compound. To describe the metabolism of MeIQx with a numerical model based on ordinary differential equations (ODE), we first converted each compound percentage in molar concentration by using a proportional relationship: ${\left[\text{compound}\right]}_{t=24}=\frac{1}{100}{\left[\text{MeIQx}\right]}_{t=0}\ast {\left(\text{percentage of total compound concentration}\right)}_{t=24}$.

In order to decrease the final ODE model variability carried out by too large a number of parameters, we operated several reductions based on the following hypotheses:

1. It is assumed that 8-CH₂-OH-IQx has a transitory effect because it was not detected at t = 24 h. Consequently, the pathway that hydroxylates MeIQx into 8-CH₂-OH-IQx and IQx-8-COOH were merged in one generic reaction which produces a generic compound named “C-Hydroxy-MeIQx”. This generic compound gathers all carbon C8 hydroxylated compounds (Fig. 3).

2. Assuming that the formation of N-desmethyl-7-oxo-MeIQx and 7-oxo-MeIQx is mediated by a same enzyme, probably xanthine dehydrogenase, we introduced a novel generic compound “oxo-MeIQx” which gathers the compounds with an oxo chemical group. This generic compound oxo-MeIQx is produced by a single generic reaction with MeIQx substrate.

3. It is assumed that HONH-MeIQx has a transitory effect because it was not detected at t = 24 h.

4. In order to model a closed system, we wish to assume that all MeIQx present at t = 0 was either recovered or transformed and distributed between the various metabolites of the system. To that purpose, we introduced a last variable in the system, named “Potential-Genotoxic-Compound”. Its concentration was estimated by summing up all concentrations at t = 24 h to the initial dose of MeIQx. Notice, however, that according to the literature, the main non-measured metabolite in the system is Ester-O-NH-MeIQx. Therefore, the variable ”Potential-Genotoxic-Compound” corresponds to the Ester-O-NH-MeIQx and the potential unidentified compounds.

As illustrated in the reduced model (Fig. 3), the concentrations of the new generic compounds were computed as follows.

${\displaystyle \left[\text{C-Hydroxy-MeIQx}\right]=\left[\text{IQx-8-COOH}\right]+\left[{\text{8-CH}}_2\text{-OH-IQx}\right]}$ ${\displaystyle \left[\text{oxo-MeIQx}\right]=\left[\text{7-oxo-MeIQx}\right]+\left[\text{N-desmethyl-7-oxo-MeIQx}\right].}$

Together, we obtained an exhaustive estimation of the concentrations of the different compounds involved in the numerical model at time t = 24 h: MeIQx, C-Hydroxy-MeIQx, MeIQx-N²-Gl, MeIQx-N²-SO3H, oxo-MeIQx, HONH-MeIQx, HON-MeIQx-N²-Gl and Potential-Genotoxic-Compound. Quantities for each compounds are depicted in Table 1 except HONH-MeIQx that is assumed to be zero because of its transitory effect.

As a main output of the model, we aimed at comparing the bioactivation pathways with the detoxification pathways according to MeIQx concentration. To that purpose, we introduced the ratio $\frac{\text{[Bioactivation]}}{\text{[Detoxification]}}$ based on the two variables [Bioactivation] and [Detoxification] computed as follows.

${\displaystyle \left[\text{Detoxification}\right]=\left[{\text{MeIQx-N}}^2\text{-Gl}\right]+\left[{\text{MeIQx-N}}^2{\text{-SO}}_3\text{H}\right]+\left[\text{oxo-MeIQx}\right]+\left[\text{C-Hydroxy-MeIQx}\right]+\left[{\text{HON-MeIQx-N}}^2\text{-Gl}\right]}$ ${\displaystyle \left[\text{Bioactivation}\right]=\left[\text{Potential-Genotoxic-Compound}\right].}$

Non parameterized numerical model for MeIQx metabolism

According to the data set and the simplifications depicted in the Method section, we built a ODE numerical model to describe MeIQx metabolism. The model contained eight variables: MeIQx, C-Hydroxy-MeIQx, MeIQx-N²-Gl, MeIQx-N²-SO3H, oxo-MeIQx, HONH-MeIQx, HON-MeIQx-N²-Gl and Potential-Genotoxic-Compound. For the three reactions regulated by an enzyme which was not involved in any other reaction, the concentration of the product was mainly assumed to be produced with a Michaelis-Menten dynamics. Each Michaelis-Menten function was parameterized by a catalytic constant K_m and a maximum rate V_m, leading to six parameters: Vm_{MeIQx−N²−SO3H}, Vm_oxo−MeIQx, Vm_PotGenotoxC, Km_{MeIQx−N²−SO3H}, Km_oxo−MeIQx, Km_PotGenotoxC.

A main characteristic of the model is that both CYP1A2 and UGT catalyze pairs of reactions leading to competition for the enzyme resources. First, CYP1A2 metabolizes its substrate MeIQx into two products P1 = C-Hydroxy-MeIQx and P2 = HONH-MeIQx. Second, UGT metabolizes both the substrates S1 = MeIQx and S2 = MeIQx-N²-Gl into the products P1 = MeIQx-N²-Gl and P2 = HON-MeIQx-N²-Gl, respectively.

In both cases, an enzyme E catalyzes two reactions R1 and R2. In that context, we introduced two parameters, named distribution coefficients, α (for the CYP1A2 competition) and β (for the UGT competition) to take into account the proportion of total enzyme concentration used in the reactions:

1. α[CYP1A2] and (1 − α)[CYP1A2] denotes the quantity of CYP1A2 involved in the catalyzing MeIQx into C-Hydroxy-MeIQx and HONH-MeIQx, respectively.

2. β[UGT] and (1 − β)[UGT]) denotes the quantity of UGT involved in the catalyzing MeIQx into MeIQx-N²-Gl and of HONOH-MeIQx into HON-MeIQx-N²-Gl, respectively.

Notice that the underlying assumption is that the total enzyme amount of CYP1A2 and UGT is limited. This assumption is supported by former analyses of enzymatic activities estimated for CYP1A and 1B using ethoxyresorufin as a substrate. These experiments evidenced that enzymatic activities showed relatively constant values during the culture (see Table 2 in Langouët et al., 2001). The corresponding generic ODE models are depicted in Fig. 4. In addition to the distribution coefficients α and β, modeling competitions required to introduce other parameters: four catalytic constants (k₁, k₂, k₃, k₄), the Michealis-Menten constants (Km₁, Km₂, Km₃, Km₄), and the total concentration of CYP1A2 and UGT.

To conclude, the model (Fig. 3) reports the dynamics of eight variables according to 16 fixed parameters (catalytic constants, maximal rates, and CYP1A2 and UGT total concentrations) and two additional parameters named distribution coefficients (α, and β) which may vary according to the experimental setup.

Modeling the enzyme competition process according to the initial MeIQx concentration

To evaluate the dynamics of CYP and UGT enzyme distributions between the two metabolic pathways, we implemented different models: a classical Michaelis-Menten model testing the saturation hypothesis and a dynamic model in which the amount of enzyme involved in a pathway can change dynamically in function of the substrate concentration. The main difference between the two hypotheses is related to the dependency with respect to the dose of the coefficients α and β modelling the enzyme distribution coefficients of CYP1A2 and UGT.

• The saturation hypothesis assumes that both the distribution coefficients α and β are fixed with respect to the concentration of [MeIQx].

• The dose-dependent hypothesis assumes that both the distribution coefficients α, β depend on the concentration of [MeIQx] at any time point, that is, α = φ(a₁, b₁, θ₁, n₁, [MeIQx]) and β = φ(a₂, b₂, θ₂, n₂, [MeIQx]). To model these relations of dependency, we used as generic sigmoid function $\phi \left(a,b,\theta ,n,S\right)=a\frac{{\theta }^n}{{\theta }^n+S^n}+b$, where a and b (0 ≤ a, b ≤ 1) are the boundaries of the sigmoids; n = 5 describes the slope of the sigmoid; θ is the inflection point of the curve.

Model fitting, parameters estimation, and a posteriori filtering

The parameters of all ODEs numerical models were selected such that the model predictions optimally fitted with the data set from the experiment HL-1 (Table 1) at t = 24 h.

For this, we denoted by t₁ = 0 and t₂ = 24 h the two time-points introduced in Table 1. We also introduced d₁ = 1 µM, d₂ = 10 µM and d₃ = 50 µM the initial concentrations of MeIQx used in experiments. We also denoted by C₁, …, C₈ the compounds that include 7 compounds from Table 1 and HONH-MeIQx that is not detected. In this setting, we denoted by $y_{k,l,i}^{obs}$ the measured concentration of each compound C_k at time t_l in HL-1 hepatocytes in presence of initial MeIQx concentration d_i. Together, we obtained 48 experimental data points.

All models describe the variations of the eight variables C₁, …, C₈ according to time. Each model is parameterized by a set of parameters θ, including (i) the 16 fixed parameters of the model and (ii) the parameters allowing to specify the values of the distribution parameters α and β according to the saturation and dose-dependent hypotheses. The prediction concentration of a compound C_k at time t_l in HL-1 hepatocytes in presence of initial MeIQx concentration d_i is denoted by y_k(θ, t_l, d_i).

The PottersWheel toolbox for Matlab software (Maiwald & Timmer, 2008) was used to sample the family of all possible parameters set θ. The main steps of the method are as follows. First, a random set of parameters is selected. From this set of parameters, a trust-region algorithm is applied iteratively to decrease the following mean-square distance (Raue et al., 2009), where σ_k,l,i’s depict an internal error linked to the parameter estimation process. ${\displaystyle {\chi }^2\left(\theta \right)=\sum _{k=1}^8\sum _{l=1}^2\sum _{i=1}^3{\left(\frac{y_{k,l,i}^{obs}-y_k\left(\theta ,t_l,d_i\right)}{{\sigma }_{k,l,i}}\right)}^2.}$

When the trust-region algorithm has reached a local minimum, the PottersWheel toolbox reports a p-value which allows evaluating the goodness of the selected set of parameters. The set of parameters was considered as relevant if the p-value was equal to 1 (most stringent case).

The added value of PottersWheel is to iterate this procedure in order to avoid bias introduced by local minima and to explore the full space of set of parameters. To that goal, the former set of parameters identified by the PottersWheel toolbox was perturbed enough to allow the exploration of new regions of the search space. The trust-region algorithm was again applied iteratively to reach another local best-fit set of parameters which is also stored as a relevant solution when its associated p-value equals 1. This procedure was iterated 1,000 times for both the saturation and dose-dependency hypothesis models. To improve the performance of the algorithm, data were transposed in a logarithmic scale prior to applying the fitting procedure. All the parameterized values used to run PottersWheel are detailed in Table S2.

Therefore, the training methods enabled the minimization of the distance between the 48 experimental data points $y_{k,l,i}^{obs}$ and the parameterized model predictions y_k(θ, t_l, d_i) and resulted in a family of sets of parameters θ.

In a second step, the family of reported suitable sets of parameters θ was filtered in order to eliminate false positive results. This was done by simulating the behavior of the family of parameterized numerical models with the “radau5” solver associated to the default thresholds Absolute Tolerance =1 × 10⁻¹⁰ and Relative Tolerance =1 × 10⁻⁶. Based on this procedure, we retained only the models with the following characteristics: (i) the formation of all variables follows a monotone dynamics, (ii) In order to limit the variability of models on each single data point, we constrained the species that are the most variable with respect to their measurements, e.g., C-Hydroxy-MeIQx and [Detoxification] at initial dose d₁ = 1 µM, as follows: ${\displaystyle \left|\frac{y_{\text{C-Hydroxy-MeIQx}}\left(\theta ,t_2,d_1\right)-y_{\text{C-Hydroxy-MeIQx},2,1}^{obs}}{y_{\text{C-Hydroxy-MeIQx},2,1}^{obs}}\right|\le 0.10,-0.05\le \frac{\left[\text{Detoxification}\right]\left(\theta ,t_2,d_1\right)-{\left[\text{Detoxification}\right]}_{2,1}^{obs}}{{\left[\text{Detoxification}\right]}_{2,1}^{obs}}\le 0.2;}$ where ${\left[\text{Detoxification}\right]}_{2,1}^{obs}$ denotes the value of detoxification at time t₂ = 24 h for the initial dose d₁.

Finally, the Detoxification pathway, the formation of Potential-Genotoxic-Compound, and consequently the bioactivation of MeIQx were estimated by simulating the final family of parameterized models for MeIQx initial concentrations equal to 0.05, 0.1, 0.5, 1, 2, 5, 10, 15, 30, 50, 75 and 100 µM during 120 h, in order to evaluate the balance between the two pathways.

Model-fitting suggests that the dose-dependance hypothesis is more relevant than the saturation hypothesis

As detailed in Material and Methods, in order to study the metabolism of MeIQx in humans, we set up a scheme of numerical models using ordinary differential equation (ODE) based on the data from the experiment using human hepatocytes HL-1 reported in Langouët et al. (2001). Two hypotheses were introduced to model the enzyme competition process according to the initial MeIQx concentration and tested by first searching for the sets of parameters which best explain the data and then testing their relevance with individual data.

As shown in Table 2, the saturation hypothesis (line 1) cannot be explained with the family of parameters identified by the fitting procedure. Indeed, all models reported by PottersWheel toolbox procedure had a p-value lower than 1, which is too low to be considered as relevant. An explanation can be found in the study of local Chi²-scores involved in the computation of the p-value shown in Table 2. Indeed, four local Chi² scores are being found outside the deviation area. In particular, if the distribution parameters α and β are fixed, the local score for MeIQx-N²-Gl and C-Hydroxy-MeIQx is too distant from biological data, respectively, to retain the saturation hypothesis.

On the opposite, the sets of parameters found for the dose-dependent hypothesis (line 2) can reproduce and explain the biological data from Hl-1 experiment. The parameter searching procedure resulted in a family of 260 sets of parameters (also called “fit”) which all equivalently best-fit the data and have a p-value equal to 1. Their filtration according to biological criteria resulted in a family of 73 sets of parameters which equivalently explained the HL1 data set and are based on the saturation hypothesis. These results suggest that the distribution of CYP1A2 and UGT between their respective pathways dynamically depends on the concentration of MeIQx. Additionally our data demonstrate that modeling such dependency with sigmoid function is appropriate describing the dynamics of biological processes.

In the following, the 73 corresponding numerical models were used to predict MeIQx metabolism activation at various concentrations. This resulted in 73 prediction curves, each corresponding to a set of parameters.

Model validation based on literature datasets

Among the detoxification metabolites, the major one is the C-Hydroxy-MeIQx derivatives (Langouët et al., 2001). We first explored the behavior of C-Hydroxy-MeIQx and HONH-MeIQx-N²-Gl for short time and low concentrations, corresponding to the usual human exposures. For this purpose, we used different values of parameters corresponding to the 73 valid numerical models.

For MeIQx initial concentrations between 0.05 to 0.5 µM, the formation of C-Hydroxy-MeIQx reaches a plateau (steady-state) before 7 h of treatment (Fig. 5). At this time, for any initial dose less than 0.5 µM, the predicted quantity of C-Hydroxy-MeIQx varies from 30% to 70% of the initial dose. This is consistent with data from Gu et al. (2010) reporting that C-Hydroxy-MeIQx accounted for 35.8% to 73.1% (see Table S1) of the dose of MeIQx measured in urine 10 h after ingestion of 920 ng of MeIQx by human volunteers.

Similarly, when the initial dose of MeIQx varies from 0.5 µM to 1.5 µM, the predicted quantities of HONH-MeIQx-N²-Gl reaches a plateau between 7 h and 12 h of treatment. During this time interval, the predicted quantity varies from 6% and 16% of the initial dose. This is consistent with data from (Stillwell et al., 1999) reporting that HONH-MeIQx-N²-Gl accounted for 2.2 to 17.1% of the dose of MeIQx measured in urine at 12 h for an exposition between 1,500 and 3,000 ng (that is 0.7 to 1.4 µM) (Fig. 6). These comparisons provide an indirect validation of our modelling approach and highlight the high level of variability in the measured response of the system, which is properly reflected by the choice of a family of ODE models instead of a single one.

Evaluation of the balance between bioactivation and detoxification

As shown in Fig. 3, the global MeIQx metabolism can be decomposed into detoxificated and bioactivated pathways. The detoxification pathways include C-hydroxy-MeIQx together with all phase II conjugates (green compounds in Fig. 7). In regards, bioactivation corresponds to the formation of Potential-Genotoxic-Compound, that is, Ester-O-NH-MeIQx and the potential unidentified metabolites (red compound in Fig. 7). In order to evaluate MeIQx bioactivation according to its initial concentration, we used our 73 mathematical models to simulate the range of detoxification and bioactivation pathways of MeIQx from 0 to 180 h, when all models have reached a steady-state. For this purpose, we again used different values of parameters corresponding to the 73 valid numerical models.

Figure 7 illustrates the distribution of MeIQx at steady-state into bioactivation and detoxification products expressed either in percentage (Fig. 7A) or in concentrations (Fig. 7B, µM). Our results show that the detoxification pathway always predominates although the balance between bioactivation and detoxification changes in function of MeIQx concentration. More precisely, the percentage of MeIQx transformed into bioactivated pathway increases when the initial concentration of MeIQx switches from 0.05 µM and 10 µM. In contrast, for concentrations greater than 10 µM, the balance decreases progressively when the MeIQx dose increases.

In addition, as depicted in the Table contained in Fig. 7B, we noticed that the duration required to transform 50% of MeIQx corresponds to 10% to 25% of the duration required to reach a complete degradation of MeIQx. For instance, for an initial dose of 0.05 µM of MeIQx, the complete metabolism takes 32 h although only 2.5 h is needed to metabolize 50% of the MeIQx dose and 9 h are needed to metabolize 90% of the MeIQx dose (Fig. 7B).

Evaluation of the Bioactivation/Detoxification ratio

In order to figure out the modulation of the ratio between bioactivation versus detoxification pathways, we plotted the evolution of this ratio between 0.05 and 100 µM at 6 h, 24 h, 72 h and 120 h in Fig. 8. This clearly shows that, whenever the considered time point, the ratio between bioactivation and detoxication pathways is not linear but mainly concentration dependent with a maximum concentration around 10 µM. This confirms the observations of Fig. 7. In addition, although the curve profiles and their associated values are quite similar, subtle variations can be observed. For example, we noticed that the peak shape is more acute at 6 h than at 120 h. Finally, these simulations suggest that the family of 73 ODE models have similar behaviors when the dose of MeIQx is greater than 100 µM, since the variability of the predictions for such dose is very low.

Simulations of the distribution coefficients dependency to MeIQx suggest that both CYP1A2 and UGT regulation occur at low concentrations of MeIQx

To deepen the regulation mechanism, we plotted the variation of the CYP1A2 distribution coefficient α (Fig. 9A) and the UGT distribution coefficient β (Fig. 9B) with respect to the MeIQx concentration. Let us recall that all the 73 parameterized numerical models which explain the best the transformation of MeIQx into several compounds were built by introducing a MeIQx-dependent distribution coefficient for CYP1A2 and UGT modeled with a sigmoid function. The evolution of this dependency with respect to the concentration of MeIQx is depicted in Fig. 9. Our simulations suggest that the distribution coefficient for CYP1A2 ranges from 0.005 at high concentrations of MeIQx to 0.05 at low concentrations of MeIQx. This suggests that the distribution of CYP1A2 always favors the bioactivation pathway, but that for MeIQx concentrations lower than 0.5 µM, the relative quantity of CYP1A2 involved in detoxification in multiplied by a factor 10. A possible interpretation is that the low quantity of CYP1A2 allows the system to favor detoxification pathway although this pathway is kinetically disadvantaged. This is consistent with experimental measurements of enzymatic activities by Turesky et al. (2001) evidencing that the detoxification pathway kinetics is slower than the bioactivation pathway.

Similarly, our simulations suggest that the UGT distribution coefficient highly favors the direct transformation of MeIQx in glucoronic compounds with respect to the transformation of intermediary compounds HONH-MeIQx into glucuronic compounds. Although, this tendency is leverage at high dose of MeIQx, with an inflection point at 1.25 µM. This suggests that the system needs to transform very toxic intermediary metabolites HONH-MeIQx when their concentration is too high.

Together, the simulation of the distribution coefficients according to MeIQx concentration suggest that the following scenario of regulation holds. For initial dose greater than 2 µM, CYP1A2 is mainly involved in the catalysis of the bioactivation pathway whereas UGT in involved in the transformation of intermediary compounds HONH-MeIQx. When the concentration of MeIQx reaches a threshold on 1.25 µM, UGT is redistributed to favor a direct transformation of MeIQx into glucoronic compounds. Finally, when the concentration of MeIQx reaches 0.75 µM, CYP1A2 favors the catalysis of the most energetic consuming detoxification pathways to gain in efficiency. Although the thresholds (inflection points) are clearly patient dependent, it should be noted that these phenomena occur in a range of MeIQx corresponding to potential real exposure.

Influence of CYP1A2 on balance between bioactivation and detoxification

The PotterWheel matlab toolbox was used to perform a 2D sensibility analysis of the influence of each parameter over the system. As depicted in Fig. S1, this analysis highlighted that the total concentration of CYP1A2 is a key actor of bioactivation products formation. In order to illustrate and understand better the influence of this parameter, the parameter value of the total concentration of CYP1A2 was modified in the 73 models (multiplication by 10, and division by 10). The perturbed models were simulated in order to compute the balance between bioactivation and detoxification. The results of these simulations is depicted in Fig. 10. By decreasing the concentration of CYP1A2 by a factor of 10, we observe that that the detoxification pathway largely predominates, so that the bioactivation pathway is minor and accounts for less than 10% of the total metabolism. In addition, the complete MeIQx metabolism is not reached at 180 h for a concentration of 100 µM. On the contrary, by increasing the concentration of CYP1A2 by a factor of 10, we observe that that the bioactivation predominates for high concentration. At low concentration, the detoxification still predominates despite the increasing of the importance of the bioactivation pathway. This suggests that, together with the concentration of MeIQx, CYP1A2 is a major actor in the regulation of the balance between bioactivation and detoxification, and that both its concentration and its operating mode are crucial to modulate the balance.