Mathematical model of uptake and metabolism of arsenic(III) in human hepatocytes - Incorporation of cellular antioxidant response and threshold-dependent behavior.

BACKGROUND: Arsenic is an environmental pollutant, potent human toxicant, and oxidative stress agent with a multiplicity of health effects associated with both acute and chronic exposures. A semi-mechanistic cellular-level toxicokinetic (TK) model was developed in order to describe the uptake, biotransformation and clearance of arsenical species in human hepatocytes. Notable features of this model are the incorporation of arsenic-glutathione complex formation and a "switch-like" formulation to describe the antioxidant response of hepatocytes to arsenic exposure.
RESULTS: The cellular-level TK model applies mass action kinetics in order to predict the concentrations of trivalent and pentavalent arsenicals in hepatocytes. The model simulates uptake of arsenite (iAsIII) via aquaporin isozymes 9 (AQP9s), glutathione (GSH) conjugation, methylation by arsenic methyltransferase (AS3MT), efflux through multidrug resistant proteins (MRPs) and the induced antioxidant response via thioredoxin reductase (TR) activity. The model was parameterized by optimization of model estimates for arsenite (iAsIII), monomethylated (MMA) and dimethylated (DMA) arsenicals concentrations with time-course experimental data in human hepatocytes for a time span of 48 hours, and dose-response data at 24 hours for a range of arsenite concentrations from 0.1 to 10 μM. Global sensitivity analysis of the model showed that at low doses the transport parameters had a dominant role, whereas at higher doses the biotransformation parameters were the most significant. A parametric comparison of the TK model with an analogous model developed for rat hepatocytes from the literature demonstrated that the biotransformation of arsenite (e.g. GSH conjugation) has a large role in explaining the variation in methylation between rats and humans.
CONCLUSIONS: The cellular-level TK model captures the temporal modes of arsenical accumulation in human hepatocytes. It highlighted the key biological processes that influence arsenic metabolism by explicitly modelling the metabolic network of GSH-adducts formation. The parametric comparison with the TK model developed for rats suggests that the variability in GSH conjugation could have an important role in inter-species variability of arsenical methylation. The TK model can be incorporated into larger-scale physiologically based toxicokinetic (PBTK) models of arsenic for improving the estimates of PBTK model parameters.

Background

Arsenic is a naturally occurring metalloid, abundant in the earth's crust and a component of more than 245 minerals [1]. Exposure to arsenic has been associated with cancers of the liver, bladder, skin and lung [2,3]. Epidemiological studies in Taiwan, Bangladesh and India have reported adverse health effects associated with chronic arsenic exposure including; chronic obstructive pulmonary disease, non-cirrhotic portal fibrosis, hypertension and ischeamic heart disease [4]. The risk of developing serious diseases from chronic exposure to inorganic arsenic in drinking water prompted the US Environmental Protection Agency (EPA) to lower the maximum contamination level (MCL) for arsenic in drinking water to 10 ppb [5].

There are two biologically important arsenic valence states: arsenite (As(OH)3, iAsIII) and arsenate (AsO(OH)3, iAsV). Inorganic arsenic in water is largely in the form of arsenate; it is negatively charged at physiological pH and slowly taken up by cells [6]. Arsenate is rapidly converted to arsenite in vivo [7] which is taken up by cells much more quickly than arsenate [8]. Methylation of arsenicals facilitates their excretion from the cell and therefore was long considered a detoxification process, but recent evidence indicates that monomethylated (MMA) and dimethylated (DMA) arsenicals have many toxic effects including increased oxidative stress [9], chromosomal aberrations (CA), and oxidative DNA damage [10-12]. In hepatocytes, trivalent monomethylated arsenicals (MMAIII) inhibit the activity of thioredoxin reductase (TR), which is a critical antioxidant enzyme controlling the cellular redox balance [13,14].

Uptake and efflux of arsenicals occur primarily through transporter proteins. Uptake of iAsIII in hepatocytes and efflux of MMAIII to blood take place through aquaporin isozymes 9 (AQP9), a family of membrane-spanning proteins that facilitate movement of solutes down their concentration gradient. AQP9 channels are expressed at high concentrations in liver cells and have been shown to transport iAsIII when expressed both in Saccharomyces cerevisiae (yeast) and in Xenopus oocytes [15-17]. Another class of transmembrane proteins facilitating the transport of iAsIII and MMAIII across the cellular membrane of hepatocytes is glucose transporters and especially GLUT2 which is highly expressed in the liver [18,19]. Glutathione conjugated arsenicals are exported to the extracellular space via multidrug resistant proteins (MRPs) and multidrug resistant P-glycoproteins (PGPs)which are ATP-binding cassette (ABC) transporters that export solutes against their concentration gradient [19-21].

Methylation of inorganic arsenic takes place primarily in the liver and specifically in hepatocytes via enzymatic catalysis by arsenic methyltransferase (AS3MT), previously known as Cyt19, producing both mono- and dimethylated arsenicals [22-25]. Two biochemical pathways have been proposed for methylation of arsenates with a key difference in the substrate for AS3MT methylation: (a) a classical process of alternating steps of reduction and oxidative methylation where iAsIII and MMAIII are the substrates and the methylation can happen in the presence or absence of GSH [26] and (b) a process of GSH conjugation and reductive methylation where thiol-bound arsenicals (arsenic triglutathione - ATG, monomethylarsenic diglutathione - MADG) are the substrates [27,28]. GSH has a stimulatory role in both methylation pathways either as a reductant or in direct conjugation with arsenicals [29].

Arsenic activates the redox sensitive transcription factor Nuclear Factor -E2- related factor 2 (Nrf2) causing its increased nuclear translocation and Nrf2 binding to the Antioxidant Response Element (ARE) [30,31]. Arsenic activates Nrf2 in a different manner when compared to other compounds such as sulforaphane (SF) and tert-butylhydroquinone (tBHQ) enhancing the interaction of specific subunits of the E3 ubiquitin ligase [32]. It has been suggested that hepatocytes exhibit an "on-switch" antioxidant behavior when exposed to increasing arsenic doses [33], possibly a result of this Nrf2 activation.

In this study, a cellular-level semi-mechanistic TK model was developed for predicting intra-cellular concentrations of different arsenicals (trivalent and pentavalent) in hepatocytes. Currently, the only published cellular-level TK model for the uptake, biotransformation and efflux of arsenicals is the Easterling et al. [34] model schematically shown in Figure 1; they demonstrated the relative importance of transport processes affecting the accumulation of arsenicals in rat hepatocytes [34]. However, a TK model for humans is needed since the inherent variation of arsenic metabolic capacity across various organisms complicates the extrapolation of rat TK model to humans. Therefore, the mathematical model presented here was developed for humans and parameterized based on data from human hepatocytes.

Figure 1

Schematic depiction of the cellular level TK model for rat hepatocytes. The solid lines that cross the cellular membrane (oval) represent transport processes, while the solid lines within the oval represent biotransformation. The dashed line represents the inhibitory effect of iAsIII on the second methylation reaction (MMA to DMA) (Source: Easterling et al. [34])

Methods

The human TK model applies mass action kinetics in order to predict the concentrations of trivalent and pentavalent arsenicals including arsenite (iAsIII), monomethylated (MMA), and dimethylated arsenicals (DMA) in human hepatocytes. This TK model takes into account processes such as influx, efflux, methylation, oxidation and glutathione conjugation of arsenicals. Moreover, it accounts for induced cellular antioxidant response due to arsenic exposure through a "switch-like" mechanism that alters the model response above a specific threshold concentration [35].

This model has been compared with the Easterling et al. model [34] , in terms of the ability to fit to data of arsenic retention and methylation in human hepatocytes. This comparison aims to highlight the advantages of developing biologically relevant TK models based on data acquired from human cells. Further comparison of these two models in terms of their estimated parameter values aim to study the major intracellular kinetic processes that contribute to the differences in metabolism between humans and rats.

Model Development

The semi-mechanistic TK model describes arsenic transport across the cellular membrane and arsenic metabolism in hepatocytes according to the metabolic reaction cascade proposed by Hayakawa et al. [27] (Figure 2). Figure 3 presents a schematic depiction of the constituents of the TK model that are also explained in Table 1. These constituents include chemical species and enzymes, and the interactions among them. Fundamental assumptions made in the formulation of the TK model are:

Figure 2

A new metabolic pathway of inorganic arsenic biotransformation via arsenic-GSH complexes formation. This pathway includes two separate branches of arsenic biotransformation: MADG → MMAIII → MMAV and MADG → DMAG → DMAIII → DMAV (Source: Hayakaya et al. [27])

Figure 3

Schematic depiction of the cellular level TK model for human hepatocytes showing the major components. Green squares represent extracellular amounts of arsenicals, while red squares represent intracellular amounts. Ovals represent activities of proteins (AQP9, TR, MT1, MT2, MRPa and MRm) and GSH (GS-Pm, GS-Pd). Arrows and hammerheads indicate activation and inhibition respectively.

Table 1

Optimized parameter values of the TK model along with the corresponding process they describe

Notation	Parameter	Value	Units	Description
1	kATG int	0.25	1/min	Rate of production of ATG
2	kMADG int	90	1/min	Rate of production of MADG
3	kDMAG	0.0122	1/min	Rate of production of DMAG
4	koxd	33.254	1/min	Rate of oxidation of DMAIII
5	koxm	0.2375	1/min	Rate of oxidation of MMAIII
6	kiAsIIIint	320	1/min	Rate of production of iAsIII
7	kMMAIIIint	1200	1/min	Rate of production of MMAIII
8	n3	8	-	Hill coefficient of inhibition of MMAIII production
9	Kd3	12.94	μM	Dissociation constant of inhibition of MMAIII
10	fGSHm	0.992	-	Coefficient of inhibition of MADG hydrolysis
11	kgin	704.96	1/min	Rate of GSH production increase
12	kDMAIIIint	0.8472	1/min	Rate of production of DMAIII
13	fGSHd	0.9988	-	Coefficient of inhibition of DMAG hydrolysis
14	Vmaxt2	1.237	pmol/min	Maximal rate of MADG efflux
15	Kmt2	19.47	μM	Half saturation constant of MADG efflux
16	kmin	512.27	1/min	Rate of upregulation of MRP
17	kss	4.26	1/min	Steady state rate of efflux of iAsIII
18	τe	10	min	Time constant of AQP9 inactivation
19	k0	4.2	1/min	Initial rate of efflux of iAsIII
20	kinf	1.613	1/min	Influx of iAsIII
21	kMMAext	0.0006	1/min	Rate of efflux of MMA
22	fm	0.2	-	Coefficient of efflux of MMA
23	Vmaxt1	0.35	pmol/min	Maximal rate of ATG efflux
24	Kmt1	2.3	μM	Half saturation constant of ATG efflux
25	kDMAext	0.002	1/min	Rate of efflux of DMA
26	fd	3	-	Coefficient of efflux of DMA
27	Vmaxm	50	pmol/min	Maximal rate of ATG methylation
28	Kmm	9.32	μM	Half saturation constant of ATG methylation
29	n1	1.22	-	Hill coefficient of ATG methylation
30	Kd1	0.315	μM	Dissociation constant of ATG methylation
31	Kim	1.53	μM	Inhibition constant of ATG methylation
32	Kid	1	μM	Inhibition constant of MADG methylation
33	Vmaxd	8	pmol/min	Maximal rate of MADG methylation
34	Kmd	0.034	μM	Half saturation constant of MADG methylation
35	n2	1.83	-	Hill coefficient of MADG methylation
36	Kd2	2.33	μM	Dissociation constant of MADG methylation
37	kain	1.64	1/min	Rate of AS3MT efficiency increase
38	fA	50	-	Coefficient of MADG methylation inactivation
39	kTR	0.99	1/min	Rate of TR signaling
40	kTRinact	987.13	1/min	Rate of TR inactivation
41	TRC	15	-	Constant of TR inactivation
42	N	1.75	-	Hill coefficient of TR inactivation
43	TR0	9.27	1/min	Steady state value of TR activity

1) Arsenite influx across the cellular membrane pores is governed by their electrochemical potential, and can be described through an ion channel conductance-based formulation.

2) The oxidative stress mediated response of hepatocytes to arsenite exposure exhibits a "switch-like" behavior, and the upregulation of enzyme activities can be described through an approximate step function at a threshold concentration.

3) The methylation reactions are influenced by cooperativity phenomena as well as substrate inhibition, and can be described through a hybrid approach of Hill and Michaelis-Menten kinetics.

4) The GSH-bound hydrolysis and clearance of methylated arsenicals exhibit a threshold-dependent behavior, and can be described using a sigmoidal function.

5) Concentrations of arsenicals are uniform within the hepatocytes as well as the extracellular medium

6) All hepatocytes in the system have identical properties, are uniformly distributed in the medium, and are exposed to the same extracellular concentrations of arsenicals.

Uptake of arsenite by hepatocytes via AQP9s [36] is governed by their electrochemical potential across the cellular membrane (for simplicity we refer to AQP9s as being the ensemble of the activity of both AQP9 and GLUT2 channels). The conductance-based formulation for ion channels proposed in the Hodgkin-Huxley model [37] is used here to describe the regulation of arsenite flux by AQP9s (Equations 1-2). Specifically, the inactivation of AQP9 subunit gates during iAsIII influx is described by Equation 1b, which expresses the increased probability of these gates closing as more transport across the gates occurs.

Where, kinf represents the mass transfer coefficient for influx of arsenite in hepatocytes; represents the rate constant for arsenite conjugation with GSH to form ATG catalyzed by the Glutathione S-Transferase (GST) family of enzymes; is the rate constant for hydrolysis of ATG (reciprocal to conjugation); kSS is the steady state rate constant for efflux of arsenite that is attained at long time periods; k0 is the rate constant describing the basal activity of AQP9; and τe is the time constant governing the regulation of AQP9 gates.

Thioredoxin (Trx) Reductase (TR) is the enzyme that catalyzes the reduction of Trx. Trx is a critical antioxidant protein and an important reductant in the methylation of arsenic by AS3MT [38]. The inactivation of TR by MMAIII leads to signals that account for two different phenomena: the induction of GSH and ABC transporters via a redox sensitive activation of the cellular antioxidant response Nrf2 nuclear receptor pathway [14,39] and the decreased methylation capacity of AS3MT. In this study the main focus is on MRPs as an efflux mechanism of arsenicals since it has been reported that they are regulated by Nrf2 [40,41]. The inactivation is modeled using principles of indirect response model theory [42,43] via the threshold-dependent parameter:. Parameter S depends on the initial exposure concentration of arsenite, (iAsIII)init and a threshold concentration. The value of S is zero when arsenite doses are below the threshold concentration, and gradually increases with greater arsenite doses. The following equations describe this reaction cascade:

where, kTR is the first-order rate constant controlling the activity of TR; TR0 is the baseline activity value of TR; is the first-order rate constant for TR inactivation; TRc is a dimensionless inactivation constant; N is the Hill coefficient for enzyme inactivation from MMAIII; , and are the rate constants governing of the activities of GSH, MRP and AS3MT, , and are the corresponding first-order decay constants.

Methylation reactions of arsenic in the liver have been modeled in published cellular-level and whole-body PBTK models with classical Michaelis-Menten kinetics [34,44-46]. Alternatives include cooperativity models such as the classical Hill-type formulation and the more mechanistic Monod-Wyman-Changeux (MWC) model [47]. In preliminary analyses, these formulations were unable to explain the time course patterns of arsenic retention and methylation in human hepatocytes (results not shown). Therefore, an alternative, non-classical formulation was used in the TK model. In this model, the AS3MT is assumed to exhibit cooperativity and the methylation reaction rate is assumed to exhibit hysteretic sigmoidal behavior as per Frieden [48]. In this formulation, the cooperativity is described by a Hill-type formulation for Vmax that is dependent on the total ATG present in the system. The formulation accounts for the constitutive influence of GSH in the methylation reaction cascade and its role in the increase of Vmax. Moreover this Hill-type formulation for Vmax accounts for a possible colocalization of thiol-containing proteins that interact with GSH (e.g. GSTP1), MRPs and AS3MT in hepatocytes. This colocalization would control not only the production and clearance of ATG but methylation activity as well [49,50].

Previous studies have suggested that exposure of human hepatocytes to elevated doses of iAsIII (0.4 - 4 μΜ) markedly reduced the production of DMAs while at the same time increased the yields of MMAs [51]. Therefore, it is assumed here that the AS3MT inactivation signal (Equation 6) affects only the second methylation reaction rate (Equation 8a).

where, MADGui and DMAGui are the rates of arsenic methylation for the first and second methylation reactions respectively. and are the maximal rates of the first and second methylation reactions respectively; and are the half-saturation constants for the methylation reactions; Kim and Kid are the uncompetitive inhibition constants for the respective reactions n1, n2 are the Hill coefficients; and are the dissociation constants influencing the sigmoidal change in and ; and fA is a coefficient of the second methylation reaction inactivation.

MADG hydrolysis reaction is modeled using a "switch-like" formulation. For doses below the threshold, a Hill-type formulation is used (Equation 9b). Above the threshold concentration, the rates of hydrolysis of MADG and DMAG to MMAIII and DMAIII, respectively, are assumed to be attenuated due to oxidative stress-induced GSH upregulation (Equations 9-10). On the other hand, this non-linear behavior may result in an increase or a decrease of GSH depending on the concentration of iAsIII and the duration of exposure. The non-linear sigmoid function tanh(S) is used here to describe this "switch-like" behavior; this formula has been previously used [52] in a neurocomputational model to describe the non-linear threshold-dependent behavior of neuronal firing rate.

where HDm and HDd are the rates of hydrolysis of MADG and DMAG respectively. is the reaction rate constant for MADG hydrolysis (MMAIII production); is the coefficient of inhibition of MADG hydrolysis; n3 and are the Hill coefficient and dissociation constant, respectively, for the inhibition term; is the reaction rate constant for DMAG hydrolysis (DMAIII production); and is the coefficient of inhibition of DMAG hydrolysis.

Efflux of GSH (or protein) -bound arsenic adducts (ATG, MADG) is assumed to take place via multidrug resistant proteins (MRPs), and is described by classical Michaelis-Menten kinetics. Since MADG is a substrate in the dimethylation reaction (Equation 8a), its efflux rate is assumed to be affected by MRP levels [51].

where, and are the Michaelis constants of the biophysical clearance of ATG and MADG, respectively; and are the rate constants of MMA and DMA clearance, respectively; fm and fd are the dimensionless coefficients of clearance of the respective processes affecting the maximal efflux.

The remaining biotransformation reactions include a series of methylation, glutathione conjugation and oxidation reactions [27] (Equations 17-23).

where is the rate constant of MADG production catalyzed by the Glutathione S-Transferase (GST) family of enzymes; koxm is the rate constant for MMAIII oxidation; kDMAG is the rate constant for DMAG production catalyzed by the Glutathione S-Transferase (GST) family of enzymes and koxd is the rate constant for DMAIII oxidation. It has been suggested for this biotransformation pathway that trivalent arsenicals mostly bound to thiol-containing proteins are conjugated with GSH and methylated in the presence of arsenic methyltransferase (AS3MT) [28]. Therefore, the parameters kATGint, and kDMAG indirectly represent the binding of trivalent arsenicals to thiol-containing proteins in the TK model.

The TK model has been implemented in MATLAB; the system of ODEs comprising the TK model is solved numerically using the stiff solver ode15s. First, the model parameters corresponding to low doses (i.e. below the threshold) were estimated using time course in vitro measurements of arsenicals following exposures to 0.1 μM iAsIII , from Styblo et al. [51]; it was assumed that at this dose, the hepatocytes exhibit no induced antioxidant response. Subsequently, model parameters corresponding to a wider range of doses (i.e. including both low-dose and high-dose behavior) were estimated using dose-response data (for doses ranging from 0.1 - 10 μM) reported by Drobna et al. [53]. The Drobna dataset includes measured concentrations of iAsIII, MMA, and DMA in primary cultured human hepatocytes after 24 h exposure to iAsIII (data for hepatocytes from 8 donors). For this case study, data on hepatocytes from one donor (white female, aged in the 60 s, Donor C) [53] were used, since this donor had similar characteristics to the human donor in the study by Styblo et al. [51] . The deterministic optimization function fmincon was used for parameter estimation in both cases.

Sensitivity analysis

Sensitivity analysis provides estimates of how variation of model's output can be apportioned to different sources of variation in model parameters. This quantity is given by the formula:

where Y denotes model output and P denotes the vector of model parameters; DP and D illustrate the partial and total variance of the model output due to variation in model parameters according to assigned statistical distributions.

The Fourier Amplitude Sensitivity Test (FAST) decomposes the total variance of model output (D) into terms of increasing dimensionality. FAST computes the Total Sensitivity Indices (TSI), which account not only for the variance due to individual parameters (Di), but also estimate the variance due to interaction among parameters (Dij, Dijk, etc.) The total variance for n dimensions is given by

The model parameters were assumed to be normally distributed with a coefficient of variation up to 10%; for some parameters, the coefficient of variation was 1%. Ten thousand (10,000) samples were generated and the normal distributions for all parameters were truncated at 1% and 99% (approximately three standard deviations from the mean value). Three model outputs were selected for the sensitivity analysis: Areas under the Curve (AUCs) of total retention of MMA, DMA and iAsIII in human hepatocytes. The SIMLAB modeling platform [54] was used to perform the global sensitivity analysis.

Comparison with the TK modeling formulation for rat hepatocytes

The TK model presented here was also compared with results from a published model for rat hepatocytes, in order to assess the inter-species differences and the feasibility of direct, cross-species extrapolation. Specifically, the comparison focused on major intracellular processes that influence the different metabolizing rates between these two species. First, the TK modeling formulation of Easterling et al. [34] was parameterized using data from Styblo et al. [51]. Although, a direct comparison is not possible because of major differences in model structures, a subset of parameters was selected for comparison that describe three major biochemical processes which account for similar cellular phenomena in both models. These processes include transport of arsenite across cellular membrane, methylation of arsenic, and biotransformation of AS3MT substrate (iAs in the rat model, and ATG in human TK model). Specific parameters used in comparison include (a) normalized activity of AQP9 (NP), which is defined as the ratio of influx and efflux of arsenite, (b) normalized activity of AS3MT (NM), which is defined as the ratio of the corresponding reaction parameters, and (c) bioavailability of AS3MT substrate (BMRS), which is defined as the ratio of rates of hydrolysis and conjugation of iAs for the human model, and the dissociation constant for the protein binding of arsenic in the rat model. These parameters are specified by Equation 26, and are described in Table 2. Parameters governing the efflux of arsenic were not compared because there is no direct correspondence between the modeling formulations.

Table 2

Comparison of selected parameter values between two hepatocyte-level TK modeling formulations

Kinetic Process	Human-Hepatocyte TK model	Rat-Hepatocyte TK model
NP - iAsIII	0.38	0.38
NM - 1st Reaction (1/min)	3.35	0.02
NM - 2nd Reaction (1/min)	147.06	0.17
Biotransformation - 1st MRS	0.0008	33.33

The two formulations are developed based on data of arsenic uptake and biotransformation in human and rat hepatocytes (NP: Normalized AQP9 activity, NM: Normalized AS3MT activity, MRS: Methylation Reaction Substrate)

In order to facilitate direct comparison, approximate volumes of the cell cultures used in different experiments (per well) were estimated. Hepatocytes were assumed to have a spherical shape with 25 μm diameter for both humans and rats [55]. Cellular volumes (Volc) for human and rat were estimated to be 1.6 and 0.8 μL respectively, based on the number of cells used in human hepatocytes experiments (2*105) and rat hepatocytes experiments (105) [34,51].

Results and Discussion

The semi-mechanistic TK model was parameterized using the fmincon function in MATLAB and time course data of arsenicals in human hepatocytes from Styblo et al. [51]. The parameters are shown in Table 1. The model was able to capture the three distinct modes of the time course patterns corresponding to experimental data (shown in Figure 4, Row 1). In contrast, the only currently existing cellular level TK model for arsenic, from Easterling et al. [34], was parameterized using the same optimization technique and data, but the model was not able to adequately capture these modes (shown in Figure 4, Row 2).

Figure 4

Predicted time course profiles of arsenicals in human hepatocytes using two TK model formulations. The top row shows estimates from the TK model presented here, while the bottom row shows estimates from the TK model formulated using the Easterling et al. [34] approach. The left column shows intracellular levels of iAsIII, total MMA, and total DMA, while the right column shows corresponding extracellular levels. Experimental data are from Styblo et al. [51] (exposures to 0.1 μM arsenic for 2 days)

The time-course estimates from the TK model show that initially (within first minutes of exposure) the rate of influx of AQP9s is substantially higher than the metabolism, thus leading to a fast accumulation of arsenite inside the cells. Then, the influx is reduced, and metabolism increases, thus leading to a slow decrease in arsenite levels (till 8 to 9 hours). During this period, MMA production appears to be the dominant process, as shown by higher levels of MMA compared to DMA, attributable to the high rate of MADG hydrolysis. Subsequently, the arsenite concentrations decrease at a faster rate, the second methylation reaction becomes more significant, and MADG hydrolysis is inhibited (Equation 9).

Figure 5 shows the dose-response profiles estimated by the TK model parameterized using Drobna et al. data [53]. The model explained the dose-response profiles in the data, and captured the significant decrease in DMA amounts at higher arsenite doses. Based on the sensitivity testing the threshold concentration value of 0.1 μΜ was able to adequately explain the arsenicals retention and metabolism, as shown in Figure 6. Threshold values above 0.1 μΜ overestimate the concentration of DMA in hepatocytes by one order of magnitude in the low dose region. On the other hand, threshold values below 0.1 μΜ (e.g. 0.01 μΜ) underestimate the DMA concentration substantially.

Figure 5

Predicted dose-response profiles in human hepatocytes using the cellular level human TK model. The left panel shows total amounts (in hepatocytes and the medium) of iAsIII, MMA and DMA. The right panel shows corresponding intracellular levels. Experimental data are from Drobna et al. [53] for hepatocytes from a 63 year old white female

Figure 6

Sensitivity testing of TK model predictions for varying assumptions of threshold levels. The left panel shows total amounts (in hepatocytes and the medium) of DMA, while the right panel shows intracellular levels of DMA. Experimental data are from Drobna et al. [53] for hepatocytes from a 63 year old white female

Results of the sensitivity analysis showed that the relative contribution of variance of individual TK model parameters varied significantly across different doses of arsenite. As shown in Figure 7A, at low doses (0.1 μΜ), the transport parameter k0 (Parameter 19) contributes most to the variance in intracellular MMA levels. This agrees with Easterling et al. [34], who reported that the transport parameters are the most significant in relation to intracellular concentration of arsenicals. On the other hand, at higher doses, the parameters related to intracellular biotransformation of MMA are the most influential. For 1 μM dose of iAsIII the most significant parameters are koxm (Parameter 5), (Parameter 10), (Parameter 7) and (Parameter 15). The first three parameters (Parameters 5, 10, and 7) directly influence oxidation and glutathione conjugation reactions involving MMAIII, whereas the Michaelis constant (Parameter 15) controls the activity of MRPs that efflux MADG from the cells. At high doses, induced antioxidant response of hepatocytes to arsenic leads to increased production of GSH and MRPs in the cells, leading to higher production of MADG, which can be readily effluxed via membrane-associated proteins. On the other hand, when the oxidation reaction is dominant (koxm; Parameter 5), it results in higher production of MMAV, which becomes accumulated in the cells, leading to overall increase in intracellular MMA levels.

Figure 7

Sensitivity analysis of intracellular MMA, DMA and iAs. Total Sensitivity Indices (TSI) of TK model parameters with respect to intracellular MMA (A), DMA (B) and iAsIII (C) levels for three dose scenarios: 0.1 μM - top, 1 μM - bottom left, 10 μΜ - bottom right

Figure 7B shows sensitivity analysis results for intracellular DMA levels. Similar to results from sensitivity analysis of MMA levels (Figure 7A), the transport parameter k0 (Parameter 19) is the most influential parameter at low doses (0.1 μΜ). For 1 μM dose, DMA production is significantly influenced primarily by (Parameter 13), which affects the rate of DMAG hydrolysis. Furthermore, the oxidation reaction is not important in the case of DMA levels, because DMA transport across the cellular membrane is much faster compared to MMA transport [56]. This is also corroborated by the relative values of corresponding transport parameters for MMA and DMA (fd >> fm,, as shown in Table 1), and by the findings of Styblo et al. [51].

Figure 7C shows sensitivity analysis results for intracellular iAsIII levels. At low doses (0.1 μΜ), the TSIs of most parameters are close to 1, indicating very high contributions. This unusual finding can be attributed to large interaction effects among multiple model parameters on the model output (i.e. binary interactions terms such as Dij,,and tertiary interaction terms such as Dijk). This was verified by computing first-order sensitivity indices, which account for contribution of each individual parameter (Di) to the output variance without taking into account higher-order interactions [57,58]; these indices for all parameters were low (< 0.1) at low doses (results not shown). At higher doses, (Parameter 1) and (Parameter 6) are the most influential. Both these parameters correspond to rate constants in the bidirectional reaction of ATG with iAsIII (glutathione conjugation of arsenite and ATG hydrolysis).

A fundamental hypothesis in this modeling formulation that allows the capturing of the dose-response profiles of arsenic retention and methylation across doses (Figure 5), is the introduction of threshold-dependent non-linear ("switch-like") mechanisms in the metabolic network due to oxidative stress (TR inactivation). This assumption is based on findings that signaling motifs exhibit biological switches under a narrow range of endogenous or exogenous stimuli [59]. This is often described by a Hill equation with a large Hill coefficient (e.g., kinase cascades [60] and nuclear-receptor pathways [61,62]). The large Hill coefficient for inhibition of MADG hydrolysis (Table 1 - Parameter 8) points to potential "switch-like" behavior of the activation of Nrf2 due to arsenic-mediated oxidative stress [32].

The parametric comparison between human and rat hepatocyte TK models for arsenic, presented in Table 2, provides insight into factors that affect arsenic metabolism in hepatocytes. Specifically, the AS3MT activity is found to be not a significant factor. Rats have been reported to be much faster metabolizers of arsenic than humans [51], but based on this study this is attributed to other factors. Specifically, in the rat-hepatocyte TK model (Figure 1), the protein-bound arsenite (p-iAs) is biotransformed to iAs (substrate of methylation reaction) at a much higher rate (five orders of magnitude) compared to the biotransformation of iAsIII to ATG (substrate of methylation reaction) of the human-hepatocyte TK model (Figure 3). This variability could be influenced by polymorphisms related to GSH production in hepatocytes [63] or availability of thiol-containing proteins to interact with AS3MT [28].

This work demonstrates the development of a prototype semi-mechanistic toxicokinetic (TK) model for arsenicals in human hepatocytes introducing features such as cooperativity and "switch-like" antioxidant response. Even though this model is not directly applicable to in vivo systems as a standalone formulation, it can be applied to inform macroscopic metabolism-related parameters in the PBTK model. On the other hand, more experimental studies on arsenicals in human hepatocytes will substantially improve model structure and can help in characterizing inter-individual variability in arsenic metabolism. Currently, the Styblo et al. study [51] is the only study in the authors' knowledge that reports time course profiles of arsenic methylation in human hepatocytes. Furthermore, significant uncertainties exist in experimental data due to the limitations of widely used techniques such as hydride generation-atomic absorption spectroscopy (HG-AAS) and high performance liquid chromatography-inductively coupled plasma-mass spectrometry (HPLC-ICP-MS), where glutathione conjugated arsenic species ATG and MADG have been reported to be degraded to iAsIII and MMAIII during the speciation analysis in the bile of rats exposed to arsenic [64,65].

This cellular-level TK model is based on an arsenic biotransformation pathway where arsenic-GSH adducts (ATG, MADG) are substrates for the respective methylation reactions [27,28]. It should be pointed out that arsenic can be efficiently methylated even in the absence of GSH [49,66], indicating that arsenic-GSH complexes need not be major species in the methylation of arsenic. On the other hand, the explicit consideration of arsenic-GSH complexes allows the description of a hysteresis behavior associated with methylation reactions and the stimulating role of GSH in these processes (Equations 7-8); so, this mechanism has been selected in this model. It should be noted that it is beyond the scope of this manuscript to comparatively evaluate the different arsenic biotransformation mechanisms.

Clearly it is, in principle, possible to incorporate into this cellular level TK model both the oxidative and reductive mechanisms as individual pathways. However, currently available experimental data are not adequate for estimating the relative contributions of each pathway. Development of improved experimental techniques for quantifying binding of arsenicals to GSH and thiol-containing proteins will allow the estimation of the relative contribution of each pathway. Since AS3MT coexists in hepatocytes with a number of competing elements that affect its action, its activity should be determined based on the availability of each of these elements. For instance, to study the effectiveness of the oxidative mechanism, it is possible to knock out GSH biosynthesis in hepatocytes by interfering with the activity of Glutamate-Cysteine Ligase (GCL) [67], and exposing them to various doses of arsenite. Such data sets can be used to estimate AS3MT activities along the two competing reaction pathways; this type of information is necessary in order to extend the mathematical formulation of the model described here to include both competing methylation pathways.

Parameter identification is an important issue in computational biology since most of the models involve more parameters than the available data. The TK model was parameterized using data on total arsenic of three species (iAs, MMA, DMA) and it was able to capture the modes of arsenic retention and methylation in human hepatocytes, but was not able to exactly capture the time-course profiles from the experimental data. In order to reduce the uncertainty associated with this issue, sensitivity analysis and testing were conducted as a means to identify the relative impact of each parameter on model predictions [68]. Additional time-course data (either on intermediate species or under more exposure/dose conditions) can improve model performance, and can help obtain additional mechanistic insights into the dynamics of arsenicals in hepatocytes.

Conclusions

A cellular-level TK model was developed based on a recently proposed pathway of arsenic biotransformation. This model can describe uptake, retention and clearance of arsenicals in human hepatocytes using a semi-mechanistic approach. It highlights the key biological processes that influence arsenic metabolism by explicitly modelling the metabolic network of GSH-adducts formation [27]. Moreover, comparison of the model structure and parameters with a rat-hepatocyte TK model [34] highlights the relative roles of different metabolic reactions in the methylation of arsenic. Ongoing work involves incorporating this cellular-level semi-mechanistic TK model as a module within a whole-body PBTK model of arsenic [44], in order to improve the PBTK model parameterization and its predictions [69].

Authors' contributions

SKS conceived the study as part of his doctoral dissertation, developed and implemented the mathematical model, analyzed the results and drafted the manuscript. CJB and SSI contributed to the development of the mathematical model, analysis of the results and manuscript drafting. PGG conceived the study and supervised the work. All authors read and approved the manuscript.