A systems biology approach to dynamic modeling and inter-subject variability of statin pharmacokinetics in human hepatocytes

Abstract

Background

The individual character of pharmacokinetics is of great importance in the risk assessment of new drug leads in pharmacological research. Amongst others, it is severely influenced by the properties and inter-individual variability of the enzymes and transporters of the drug detoxification system of the liver. Predicting individual drug biotransformation capacity requires quantitative and detailed models.

Results

In this contribution we present the de novo deterministic modeling of atorvastatin biotransformation based on comprehensive published knowledge on involved metabolic and transport pathways as well as physicochemical properties. The model was evaluated on primary human hepatocytes and parameter identifiability analysis was performed under multiple experimental constraints. Dynamic simulations of atorvastatin biotransformation considering the inter-individual variability of the two major involved enzymes CYP3A4 and UGT1A3 based on quantitative protein expression data in a large human liver bank (n = 150) highlighted the variability in the individual biotransformation profiles and therefore also points to the individuality of pharmacokinetics.

Conclusions

A dynamic model for the biotransformation of atorvastatin has been developed using quantitative metabolite measurements in primary human hepatocytes. The model comprises kinetics for transport processes and metabolic enzymes as well as population liver expression data allowing us to assess the impact of inter-individual variability of concentrations of key proteins. Application of computational tools for parameter sensitivity analysis enabled us to considerably improve the validity of the model and to create a consistent framework for precise computer-aided simulations in toxicology.

Background

The discovery and development of new drug entities is strongly handicapped by the circumstance that about 50% of the drug candidates fail in the clinical studies [1]. About one quarter of candidate drugs fail due to toxicity or pharmacokinetic (PK) problems [2], and currently, it is a well known fact, that toxicity is the major cause of attrition in the drug development process [3]. Therefore, it is quite clear that dangerous properties of drug entities have to be revealed very early in the drug evaluation studies [4]. Despite the ever growing effort to apply computational power towards improving the understanding and in-silico prediction of drug pharmacokinetics and response, the overall impact on preclinical safety testing has been modest.

Application of systems biology holds tremendous promise because it aims to understand and quantitatively describe biological phenomena within the framework of the hierarchical levels of metabolic pathways and regulatory networks at the different scales of cells, tissue, organs and ultimately whole organisms [5, 6]. However, despite emerging consensus that such a holistic approach is essential to provide the framework of predictive toxicology, the number of successful case studies is still minuscule [7–9].

Current activities can be grouped into

(1) quantitative structure-activity-relationship (QSAR) methods, computational models based on compound structure and focused on potential interactions of small molecules with major classes of proteins such as drug metabolizing enzymes [10–15], transporters [16] and receptors [16–18]. Also important are physicochemical properties of the drug, for example solubility and permeability that are estimated from the molecular structure [19–22].

(2) in vitro kinetics for prediction of in vivo drug clearance using kinetic data from recombinant cytochromes P450 (CYPs), microsomes and hepatocytes (IVIVE: in vitro-in vivo extrapolations) [23].

(3) physiologically based PK (PBPK) modeling [24–28] which considers the anatomical, physiological and chemical aspects of ADME (absorption, distribution, metabolism and elimination of the drug) [29–31] in multi-compartment models [32].

In addition to these simulations based on mathematical models various computational and bioinformatics approaches are applied to extract information from high throughput data of drug response experiments at cellular, tissue, organ and whole organism level.

A critical assessment of the aforementioned tools, essentially to outline gaps that must be addressed for more reliable predictive simulation-based toxicology, indicates needs for more rigorous network models focusing at systems dynamics beyond kinetics of individual enzymes, consideration of inter-individual variability and systematic investigations of parameter sensitivity and its impact on model verification, discrimination and reduction, to name a few. The first issue is related to the design of the dynamic models for the drug elimination process in the hepatocyte, which should be based on the integration of membrane transport processes for substrates and products as well as phase I and phase II reactions. These models need to be tightly linked to stimulus (dose)-response experiments.

The issue of model parameterization in the context of modeling in toxicology has been already addressed in 1995 by Andersen et al [24]. In addition to the problem of identifiability, particular attention should be given to correlation between parameters, very common in biological systems.

Yet another question of interest concerns the subtle integration of the enormous inter-subject variability in enzymatic phenotypes into the model. This is of outermost importance for predictions in toxicology and also in clinical pharmacology in order to design optimal treatments for individual patients. Consideration of this variability in phenotype should rest on quantitative proteomics and activity data from human liver tissue or hepatocytes representing a statistically significant portion of the population.

The importance of this issue has been extensively addressed in the review of Rostami-Hodjegan and Tucker [23] and discussed in the context of IVIVE (in vitro-in vivo extrapolation) approach. Significant progress has been made by this group through a strategy which is based on the simulation of drug disposition in virtual liver populations. The basic idea of this method is to portray the variability into the r _max /K _M -values estimated from in vitro studies using human liver microsomes, primary or cryopreserved hepatocytes or recombinant expressed enzymes. This is a promising approach; however, improvement of the strategy is imaginable by separation of the two parameters r _max and K _M (maximal enzymatic rate and enzymatic affinity parameter in Michaelis-Menten- and derived kinetics).

This concept, which is a focus of the present article, is driven by the possibility to incorporate separate information from pharmacokinetics and quantitative proteomics as well as future incorporation of the regulatory network responsible for variation of expression level of the enzymes, transporters and receptors. As model drug we chose atorvastatin (AS), one of the most frequently used 3-hydroxymethylglutaryl coenzyme A reductase inhibitors (statins). Statins thereby reduce cholesterol synthesis and they also stimulate the uptake of LDL-cholesterol from the blood. Although they are relatively safe drugs, the lipid-lowering effect of statins is inadequate in some patients, and unpredictable drug-drug interactions can occur, as well as hepatic and extrahepatic adverse effects including hepatotoxicity, myopathy, and rare but severe rhabdomyolysis.

Some aspects of dynamic modeling of statins have been tackled in previous studies, including the deterministic modeling of transport kinetics [33, 34] and the pharmacokinetic modeling of the clearance mechanisms including consideration of unspecific binding effects [35]. In the latter case, a multi-compartmental model approach has been derived on isolated rat hepatocytes. However, the metabolism of AS inside the liver cell has not been described in details, missing the kinetic description of the different metabolite formations, which thus far precludes the consideration of inter-individual variability, apart from the fact that rat hepatocytes had been investigated.

This study describes a deterministic modeling approach of the dynamic biotransformation and transport processes of AS in human hepatocytes considering both kinetics of transport processes as well as intracellular detoxification processes via Phase I and Phase II enzymes. Attention is given to the detailed enzymatic and chemical reactions, including conversions between the acid and lactone form of AS as well as unspecific binding between metabolites and macromolecules like proteins. The additional integration of quantitative protein data of Phase I and Phase II enzymes enables the dynamic analysis of the underlying inter-individual-variability of expression levels of these enzymes.

Methods

Chemicals

Williams Medium E (WME) without phenol-red and without L-glutamine and stabilized L-glutamine were obtained from Pan-Biotech-GmbH (Aidenbach, Germany). Penicillin/Streptomycin and ITS-X were obtained from Invitrogen (Karlsruhe, Germany). Bovine serum albumin (BSA) and dimethyl sulfoxide (DMSO) were purchased from Sigma-Aldrich Chemicals (Taufkirchen, Germany). Dexa Inject was obtained from Jenapharm (Jena, Germany). AS, its metabolites and deuterated standards were purchased from Toronto Research Chemicals Inc. (North York, Canada). PBS was obtained from Invitrogen (Karlsruhe, Germany) and Complete Mini, EDTA-free from Roche Diagnostics (Mannheim, Germany). NaPP, sucrose and LiChrosolv were purchased from Merck (Darmstadt, Germany); acetonitrile from Roth (Karlsruhe, Germany), formic acid from Fluka, (Germany), UGT1A3 monoclonal mouse antibody from Abcam (Cambridge, England). Trypsin was purchased from Promega (Mannheim, Germany). Synthetic peptides were purchased from Sigma-Genosys (Haverhill, UK).

Isolation and cultivation of primary human hepatocytes

Tissue samples from human liver resections were obtained from patients undergoing partial hepatectomy. Experimental procedures were performed according to the guidelines of the charitable state-controlled foundation HTCR (Human Tissue and Cell Research) Regensburg, Germany, and the institutional guidelines for liver resections of tumor patients with primary or secondary liver tumors, Technical University Munich, MRI, Munich, Germany. The use of human hepatocytes for research purposes was approved by the local ethics committees of the Ludwig-Maximilians-University of Munich [36] and the Charité, Humboldt University Berlin [37], Germany, and written informed consent was obtained from all patients. Hepatocytes were cultured on collagen gel precoated 6-well plates at a density of 1.5·10⁶ cells/well. Cells were allowed to attach to the collagen layer. After transport, culture media was disposed and attached cells were cultured 24 h at 37°C in a humidified chamber with 95%/5% air/CO2 in serum-free medium WME, supplemented with albumin (0.1% (v/v)), penicillin/streptomycin (100 U/ml), stabilized L-glutamine (2 mM), dexamethasone dihydrogenphosphate (0.025% (v/v)) and ITS-X (5 mg insulin, 3.35 μg natrium-selenit, 2.75 mg transferrin and 1 mg ethanolamine), further named SFM.

Time-series experiments

Incubation with AS was started by disposing the culture media and cultivation of the attached cells at 37°C in a humidified chamber with 95%/5% air/CO₂ in 2 ml SFM, supplemented with 10 μM AS, 0.1% (v/v) BSA and 0.1% DMSO. At specified time-points, three wells were further treated for the preparation of samples for the measurement of extracellular and intracellular metabolites, respectively. SFM media was collected and 50 μL formic acid and deuterated internal standard was added for the further measurement of extracellular metabolites. Cells were harvested in pre-cooled albumin-free SFM, disrupted by freeze/thaw and ultra-sonification and centrifuged. The supernatant was used for the determination of intracellular metabolites.

For the preparation of samples for the protein measurements, culture medium was disposed and cells were harvested in pre-cooled PBS, supplemented with Complete Mini EDTA-free (1 Tablet/10 ml Buffer). Cell suspensions were centrifuged 5 min (500 g) at 4°C and cell pellets were resuspended in 150 μL NaPP-buffer (0.1 M, pH 7.4), containing 250 mM sucrose and Complete Mini EDTA-free (1 Tablet/10 ml Buffer). Cells were disrupted by ultra-sonification and lyophilized for the analysis of total protein concentration, CYP3A4 and UGT1A3 content.

For cell number determination, cells from two wells were fixed with methanol-acetic acid fixative solution (10 min at 37°C and 4°C) and afterwards nuclei were stained for 15 min with Meyers Hämalaun (Sigma-Aldrich Chemie GmbH, Germany), rinsed with water and air-dried. Stained nuclei were counted in digital images (10 per well) at 40-fold Magnification (ImageJ Image Processing and Analysis Program).

Quantification of atorvastatin and its metabolites

AS and ASL, and their para- (ASpOH, ASLpOH) and ortho-hydroxy-metabolites (ASoOH, ASLoOH), were determined by LC-MS-MS analysis using the respective deuterium labeled analogues as internal standards, essentially as described [38]. HPLC separation was performed at 30°C on a XBridge Shield RP18 column (2.1 × 50 mm, 3.5 μm, Waters) using (A) 1 mM formic acid and (B) acetonitrile as mobile phases at a flow rate of 0.4 ml/min. Gradients were programmed as follows: 63% A for 4 min; linear decrease of A to 60% within 9 min; linear decrease of A to 55% within 2.5 min; 55% A for 1 min; increase of A to 63% in 0.2 min. Equilibration time of the column was 20 min. MS-MS analysis was performed on an Esquire HCT ultra ion trap mass spectrometer (Bruker Daltonics, Bremen, Germany) coupled to an HPLC 1100-System (Agilent, Waldbronn, Germany) consisting of binary pump G1312A, degasser G1379A, well-plate sampler G1367A and column thermostat G1330B. The ionization mode was electrospray (ESI), polarity positive, mass range mode ultrascan, and nitrogen was used as a drying and nebulizer gas. The following parameters were applied: nebulizer 45 psi, dry gas 10 l/min, dry temperature 300°C, capillary 4100 V, scan range 200 - 600 m/z.

Precursor and product ions (m/z) of analytes and internal standards, respectively, were ATV (559 and 440.2; 466.2), [²H₅]ATV (564 and 445.2; 471.2), ATV-L (541.2 and 448.2), [²H₅]ATV-L (546.2 and 453.2), p-OH-ATV (575 and 440.2; 466.2), [²H₅]p-OH-ATV (580 and 445.2; 471.2), p-OH-ATV-L (557 and 448.4), [²H₅]p-OH-ATV-L (562 and 448.4), o-OH-ATV (575 and 466.4), [²H₅]o-OH-ATV (580 and 471.2), o-OH-ATV-L (557 and 448.4), [²H₅]o-OH-ATV-L (562 and 448.4). Sample quantification was possible in a range from 0.5 to 500 pmol.

CYP3A4 and UGT1A3 protein quantification

Protein quantification of CYP3A4 and UGT1A3 in human liver microsomes and relative protein quantification of UGT1A3 in lyophilized samples of primary hepatocytes was performed by immunoblotting as described previously [38, 39].

Cell lysates for absolute quantification analysis of CYP3A4 were prepared from lyophilized human primary hepatocytes by sonification in the presence of glass beads in buffer containing Complete Mini-Protease Inhibitor Cocktail, and following homogenization. In the CYP3A4 quantification assay, three synthetic isotopically labeled peptides (13C/15N amino acid) were used as internal standard for calibration. These isotope labeled standard peptides represented sequence analogues to proteotypic peptides of CYP3A4, which arise from tryptic digestion. After acetone precipitation and resolving of the proteins in 8 M Urea, a definite amount of internal standard peptides was added. The sample mixture was reduced with 5 mM DTT and alkylated with 15 mM iodacetamide in 50 mM ammonium bicarbonate. Subsequently samples were digested with trypsin at 42°C for 4 h (enzyme/substrate ratio of 1:10). Efficiency of tryptic digestion was checked by SDS-PAGE followed by silver staining. The resulting peptides were purified using C18 OMIX^® Tips (Varian, Darmstadt, Germany) according to manufacturer's suggested protocol and separated on a nanoliter-flow Ultimate HPLC system (Dionex, Idstein, Germany). After injection (15 μl), peptides were trapped and desalted on a precolumn (0.3 mm I.D. × 5 mm PepMapTM, Dionex) at a flow rate of 30 μl/min in 0.1% TFA for 6 min. Peptides were transferred to the separation column (75 μm I.D. × 250 mm PepMapTM column, Dionex) and separated in a linear gradient of mobile phase (A: 0.1% formic acid, B: 84% acetonitrile/0.1% formic acid) from 5% B to 35% B over a period of 35 min with a flow of 290 nl/min. The column effluent was continuously directed into the NanoSpray II source of a 4000QTrap mass spectrometer (Applied Biosystems, Foster City, CA, USA). The MS was set up to run a multiple reaction monitoring experiment essentially, as described previously [40], including two to three parent-to-product ion transitions for each internal standard peptide as well as the corresponding transitions of native peptide of CYP3A4. The instrument settings were as follows: ion spray voltage, 3-4 kV; interface heater, 150°C; declustering potential, 50 V; collision energy, peptide specific; entrance potential, 10 V; collision cell exit potential, 10 V. MS data were processed by integrating the appropriate peak areas from extracted ion chromatograms by MultiQuantTM Software (Applied Biosystems). The absolute amount of CYP3A4 protein was calculated from the peak area ratio "internal standard peptide/native peptide". Total protein content of the samples were determined by amino acid analysis (AAA) on a Waters 2695 HPLC system using the AccQ•Tag derivatization method (Waters, Eschborn, Germany), according to manufacturer's instructions.

Identification of metabolic network structure

Numerous metabolic and physicochemical aspects about AS had to be considered in the initial model building. AS exists in two forms, a very lipophilic lactone (ASL) and a comparably hydrophilic hydroxyl-acid (AS). AS is converted enzymatically via an instable intermediate product into ASL, mediated by different UGT isoenzymes [41, 42]. Recent investigations by ourselves and others have shown that the most important contributor to UGT-driven lactonization is UGT1A3, whereas UGT1A1 plays an insignificant role [38, 42]. AS acids and lactones are inter-converted chemically into each other [43]. However, studies have indicated, that the chemical lactonization of AS to ASL can be neglected at physiological pH of 7.4 [43]. Recent studies highlight, that different PON enzymes might also be possible contributors to the lactone hydrolysis and that PON1 is present in liver [44–48].

Both AS and ASL are hydroxylated in human hepatocytes leading to para- and ortho-hydroxy-metabolites, ASpOH, ASoOH, ASLpOH and ASLoOH [49, 50], mainly catalyzed by CYP3A4 [51]. Recent studies have reported, that CYP2C8 and CYP3A5 also hydroxylate AS to a minor extent [50, 52].

Furthermore, AS is transported into the cell via organic anion transport polypeptides (OATP), and recent studies on recombinant systems showed, that OATP1B1 and OATP2B1 contribute to the AS import [53, 54], which are both expressed in the human liver [55, 56]

OATP1B3 is supposed to be also a main contributor to drug transport [54], since it shows high gene expression levels in liver [57], but its importance for AS has not been investigated kinetically so far.

OATP transporters have been reported to be bidirectional facilitated diffusion transporters, independent from ATP and Na+, K+ and H+, but with a possible involvement of reduced glutathione [58, 59]. However, previous investigations on transport mechanisms on rat hepatocytes, showed, that the intracellular concentrations of pitavastatin and other compounds are much higher than outside the cell [33, 34, 60]. Facilitated or passive diffusion do not allow a greater intracellular than extracellular concentration of the parent drug of interest, when the source is the initial extracellular concentration, because it is a concentration-gradient dependent mechanism. Therefore, the import mechanism should be rather considered as active mediated transport, which is not concentration-gradient dependent, rather than as the proposed facilitated diffusion process.

As shown with the recombinantly expressed transporters mentioned above [54], transport of the acidic metabolites, ASpOH and ASoOH, by OATP1B1 was similar to that of AS, and transport of the corresponding lactones was also mediated by this transporter, although at somewhat lower rates. We therefore assumed that the same OATP transporters are responsible for the import of AS and its hydroxylated and lactone metabolites into hepatocytes.

AS and its hydroxylated metabolites, ASpOH and ASoOH, are actively exported out of the hepatocytes into the bile by the ATP-dependent MDR1 transporter [61, 62]. In addition, acidic and lactone form of AS showed inhibitory effects in transport studies of substrates of MDR1 and MRP2 [63–65], pointing to the competitive transport mechanism of these substrates at this proteins. Further, the transporters MRP1, MRP3 and MRP6 are also reported to be responsible for the transport of organic compounds including AS from inside the hepatocytes into the plasma [56, 66].

Passive diffusion might play also an important role in the transfer of AS, ASL and the corresponding metabolites, ASpOH and ASoOH, ASLpOH and ASLoOH, respectively. Since the acidic forms of AS are rather hydrophilic and the lactone forms of AS are rather lipophilic, it can be assumed that passive diffusion plays a more important role for the lactone forms, and the transporter mediated active transport plays a more important role for the acidic forms, as reported earlier for statins [67].

Finally, lipophilic drugs have a high affinity to bind non-specifically to proteins, and previous studies concentrated on the modeling of drug binding in the intracellular and the extracellular space as well as on the surface of the cells [60, 68, 69].

Mathematical modeling

The mathematical model can be described by the system of non-linear ordinary differential equations

(1)

where the change in extracellular, intracellular or unspecific bound metabolite concentration c _j in the extra - or intracellular compartment with V _comp is effected by the conversion or production of contributing chemical or enzymatic reactions or by transport steps r _ij , respectively. The extracellular volume, , equals to the volume of the media used. The intracellular volume, , equals to the total volume of all cells used, and is determined by multiplying the cell number by the volume of a single hepatocyte, estimated to be 14.1 pL by the approximation of a spherical shape with a diameter of 30 μm [70].

Appropriate reaction kinetics r _ij are modeled for the CYP3A4 hydroxylation, the UGT1A3 lactonization, the chemical and enzymatic lactone hydrolysis and intracellular unspecific binding to macromolecules.

Previous studies determined substrate inhibition kinetics of the CYP3A4 mediated hydroxylation of AS on human microsomes [52]. However, inhibition effects contribute severely only at a concentration higher than 100 μM. Furthermore, our model approach considers the competitive nature of alternative substrates, by integrating the CYP3A4 hydroxylation of AS and ASL as reaction kinetics describing the competitive conversion of alternative substrates to alternative products, illustrated for the hydroxylation of AS to ASpOH

(2)

(see Additional file 1 - Derivation of Atorvastatin kinetics at CYP3A4).

The lactonization of AS to ASL is mediated by UGT1A3 enzymes and the reaction is formulated as substrate inhibition kinetics [41].

(3)

The lactone metabolites are either hydrolyzed chemically to the respective acid metabolites [43] inside (c) or outside (m) the cell, or enzymatically by the contribution of PON enzymes inside the cell. Both reactions are described as first order kinetics

(4)

Unspecific binding of intracellular metabolites to macromolecules is formulated as

(5)

with the dissociation coefficient k _dis and the intracellular fraction unbound [71]

(6)

which describes the ratio between the intracellular free concentration to the sum of intracellular free and bound concentration, and , in equilibrium (index eq) (see Additional file 2 - Derivation of kinetics of unspecific protein binding). However, the intracellular free and bound concentrations in equilibrium are not measurable; therefore, the fraction unbound fu _j is set as a parameter to be estimated in the parameter optimization procedure.

Transport steps include active import and export of the metabolites as well as passive diffusion steps. Both active import by OATP1B1 or OATP2B1 and export of AS are described as Michaelis-Menten-kinetics [53, 54]

(7)

whereas the active transport kinetics of the other metabolites are assumed to be of first-order [72].

(8)

Besides the active transport, metabolites undergo passive diffusion through the double-layer lipid-membrane. Passive diffusion, described as

(9)

is driven by the concentration difference between outside and inside the cell, , over the lipid-membrane with thickness d, through all cells with the total surface area A _cells , and controlled by the diffusion coefficient D _j and is comprised as the permeability coefficient P _j .

Optimization procedure for estimation of model parameters

The optimization procedure is based on evolutionary strategies which are implemented with JavaEva (WSI Computer Science Department, Center for Bioinformatics, University of Tübingen, Germany) and a MVA (main vector adaptation) mutation operator [73]. The optimization procedure estimates the parameters in equations (2) to (9) based on the optimization criterion

(10)

where the deviation of calculated and measured concentrations divided by the measurement standard deviation s _j , squared and summed over all metabolites J and all time points N, has to reach a minimum. Additional optimization constraints are fu _AS > fu _ASL , P _ASL > P _AS , P _AS > P _ASOH and P _ASL > P _ASLOH , because the lactones have a higher lipophilicity than the acids and the metabolites are supposed to be more hydrophilic than the corresponding parent lactone or acid drug.

The integration of the differential equations (1) using the reaction kinetics in equations (2) to (9) was performed by the differential algebraic equation solver LIMEX (Konrad-Zuse-Zentrum für Informationstechnik, Berlin) [74].

Relative abundance approach for prediction of r _max -parameters

For pharmacokinetic predictions taking inter-individual variability of CYP3A4 and UGT1A3 expression levels into account, maximal rate parameters are predicted via a relative abundance approach, which is based on the assumption that the maximal rate of the reaction is proportional to the enzyme concentration:

(11)

The maximal velocities of the respective enzyme e, here CYP3A4 or UGT1A3, in the conversion to the product P _j in the liver of interest li is estimated from the respective maximal rate and the enzyme concentration in the reference liver and the enzyme concentration in the liver of interest li.

Computational approach

The mathematical model of AS metabolism was coded in FORTRAN language and linked to the numerical integrator LIMEX, also written in FORTRAN. After compilation to the executable program, optimization was started by the call of JavaEva. The mathematical model is supplemented as SBML-file for review purpose.

Results

Model-setup

The model structure of AS biotransformation in primary human hepatocytes is schematized in Figure 1 and comprises the description of the extracellular, intracellular and unspecifically bound metabolites, as well as corresponding reaction and transport steps and unspecific protein binding.

Figure 1

Scheme of the full model of Atorvastatin metabolism in primary human hepatocytes. The scheme illustrates atorvastatin metabolic and transport pathways in human hepatocytes. The model comprises the active import and export of AS and ASL, and corresponding para- and ortho-hydroxy-metabolites, ASpOH and ASoOH, ASLpOH and ASLoOH, mediated by import proteins, OATP1B1 and OATP2B1 (solid lines), and export proteins, MDR1 and MRP2 (solid lines), as well as passive diffusion steps (dashed lines). AS and ASL are hydroxylated to the corresponding metabolites by CYP3A4, respectively (solid lines). Compound AS is converted via an instable glucuronid-intermediate to ASL mediated by UGT1A3 (solid lines). The lactones are hydrolyzed either chemically or enzymatically by putative PON to the respective acids (solid and dotted line). Acids and lactones are further subject to unspecific binding to macromolecules in human hepatocytes (dashed lines).

Since the calibration and quality samples were prepared in the same media as the experimental metabolite samples and because it can be assumed that the intracellular protein concentration is much higher compared to the outside of the cells, the unspecific protein binding was only considered for the intracellular space.Regarding the unspecific protein binding, it can be assumed that the intracellular protein concentration is much higher compared to the outside of the cells. This assumption is based on the estimation of the ratio of intracellular protein to extracellular media protein concentration. From the total cell protein measurement and the determination of total cell volume a total intracellular protein concentration of about 30 g/l can be estimated. The media includes as sole protein component 0.1% (v/v) albumin, corresponding 1 g/l. Therefore, this assumption is justified.

Because our hepatocyte culture model does not allow to distinguish experimentally import and export, single irreversible, actively mediated transport steps for both directions, as well as passive diffusion steps were considered for each compound, respectively, similar to mechanistic transport models of other studies [33, 34]. Only in the case of AS different transport steps for OATP1B1 and OATP2B1 were implemented in the model, because specific KMs had been determined in previous studies [53, 54].

Model verification

Model verification was performed by quantitative measurements of AS and all metabolites in the extra- and intracellular space during time-series experiments, performed first on a single batch of primary hepatocytes from a single individual 1. Special emphasis was given to the issue of parameter sensitivity and correlation between parameters to improve the quality of the model.

Time series/stimulus-response measurements on primary human hepatocytes

From the metabolite concentrations in primary human hepatocytes of individual 1 (see Additional file 3 - Atorvastatin metabolite concentrations from the time-series experiment on primary human hepatocytes of individual 1) we calculated total recovery at each time point (Table 1). The recovery is calculated from material balance equations and is defined as the sum of intracellular and extracellular metabolite amounts at the respective time-point divided by initial AS amount. As evident from Table 1, the total recovery was close to 100% after ten minutes but decreased at higher time-points. This may be explained by unspecific protein binding leading to the pool of bound metabolites, which increases over time due to intracellular accumulation of metabolites. This result points to the necessity and importance of the implementation of unspecifically bound metabolites in the model.

Table 1 Recovery of atorvastatin metabolites in the time-series experiment on primary human hepatocytes of Individual 1

Parameter estimation

Parameter optimizations were performed with aid of evolutionary algorithms (JavaEva, μ = 8 parent and λ = 4 children) and using the nominal parameter values of the reactions, transport and diffusion steps in equations (2) to (9). The optimization criterion in equation (10) was used, which demands to minimize the difference between experimental concentration data and model simulation. In the optimization procedure certain parameters were fixed to values (Table 2), which were either identified in or assumed from previous investigations. The Michaelis-Menten constants K _M of the CYP3A4 hydroxylation and of the AS transporters OATP1B1 and OATP2B1 were fixed to values determined with recombinant enzymes [50, 53, 54]. The K _M and K _I of the UGT1A3-lactonization were fixed to values determined on human liver microsomes [41]. The rate constant of spontaneous hydrolysis of ASL, k _CR , was estimated from experimental observations [43] and assumed to be same for the hydrolysis of the lactone metabolites, ASLpOH and ASLoOH, respectively. The dissociation rate constant k _dis of unspecific binding was fixed to a value, considered to be very high in a previous study on modeling of protein binding mechanisms [33].

Table 2 Parameters which were fixed in the optimization procedure to literature values

Analysis of the predicted concentration-time-profiles

The model predicted concentration-time-profiles are illustrated together with the measured concentrations (see Additional file 3 - Atorvastatin metabolite concentrations from the time-series experiment on primary human hepatocytes of individual 1) in Figure 2. Both intracellular AS and ASL are converted to the corresponding para- and ortho-hydroxy metabolites, ASpOH and ASoOH, and ASLpOH and ASLoOH, respectively. However, the acidic metabolites, ASpOH and ASoOH, show higher intracellular concentrations than the lactone metabolites, ASLpOH and ASLoOH. The ratio between intracellular AUC _{0-600 min} of the acidic form to that of the lactone form equals 20.1 in case of the para-hydroxy-metabolites and 23.6 in the case of the ortho-hydroxy-metabolites.

Figure 2

Measured concentrations and simulated fits of Atorvastatin compounds, individual 1. Displayed are the measured concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 1, and corresponding simulation profiles (solid lines). First and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH (right).

Accordingly, the extracellular concentrations of the acidic metabolites are higher than that of the lactone metabolites. The ratio between extracellular AUC _{0-600 min} of the acidic form to that of the lactone form equals 54.3 in case of the para-hydroxy-metabolites and 70.8 in the case of the ortho-hydroxy-metabolites.

Further, the results show that the components have higher concentrations in the intracellular space than outside the cell. The ratio of intracellular AUC _{0-600 min} to extracellular AUC _{0-600 min} equals minimally 4.5 in the case of ASoOH and maximally 38.2 in the case of ASLpOH. Similarly, the ratio of intracellular c _max to extracellular c _max equals minimally 2.5 in the case of ASoOH and maximally 34.5 in the case of ASLpOH.

Simultaneous model verification on different individual hepatocyte donors

The process of model verification so far has been applied on a single experiment on primary human hepatocytes of individual 1. In the next step it has to be proven, that the model is also capable to describe different individual metabolic profiles, especially being further able to reflect the inter-individual parameter variability, not only in the parameters r _max of the phase I and II reactions, but also in the transporters. Therefore, the model verification is performed based on AS biotransformation data on primary hepatocytes from three different individuals simultaneously.

Therefore, maximal rate constants r _max of CYP3A4 hydroxylation and UGT1A3 lactonization were predicted for individual 2 and 3 via the relative abundance approach (equation (11)) using the r _{max -}parameters of individual 1 (Table 3) and the protein concentrations observed on the primary human hepatocytes (Table 4). These parameters were then fixed in the optimization. Further, the rate constants k of PON mediated lactone hydrolysis of individual 2 and 3 were estimated in the optimization procedure.

Table 3 Verified parameters of CYP3A4 and UGT1A3 in individual 1

Table 4 Protein concentrations of CYP3A4 and UGT1A3 in primary human hepatocytes

However, the analysis of parameter sensitivity and identifiability (see Additional file 4 - Analysis of parameter sensitivity and following model reduction procedure) showed that the first order kinetics of OATP1B1 mediated import is favored over the zero-order kinetics of OATP2B1 import of AS, and that the first order kinetics are hard to distinguish from first order passive diffusion mechanisms in this evaluation system. Thus, in the following, mechanisms of active transport and passive diffusion are lumped, resulting in the apparent transport rate for import,

(12)

and in the apparent transport rate for export

(13)

described by the product of apparent rate constant κ _im/ex of import or export and extra- or intracellular concentration , respectively (Figure 3). Consequently, the rate constants κ of import and export are set individually, and were to be estimated in the optimization procedure.

Figure 3

Scheme of the model used in the simultaneous model verification. In contrast to the full model illustrated in Figure 1, this model contains lumped transport kinetics, which means, one import and one export rate for each metabolite, respectively. This transport simplification resulted from the local parameter sensitivity analysis (see Additional file 4 - Analysis of parameter sensitivity and following model reduction procedure).

The simultaneous model verification was performed based on the stimulus response data obtained from primary human hepatocytes of individual 1 (see Additional file 3 - Atorvastatin metabolite concentrations from the time-series experiment on primary human hepatocytes of individual 1), individual 2 (see Additional file 5 - Atorvastatin metabolite concentrations from the time-series experiment on primary human hepatocytes of individual 2) and individual 3 (see Additional file 6 - Atorvastatin metabolite concentrations from the time-series experiment on primary human hepatocytes of individual 3), respectively.

The model prediction is in satisfying agreement with experimental data of individual 2 (Figure 4). However, in case of individual 3, the model prediction is comparably poor, because there are relatively high deviations, especially with the intracellular metabolites ASL and ASoOH and extracellular ASpOH (Figure 5, solid line). One reason could be the extended contribution of beta-oxidation in AS metabolism. Beta oxidation at the heptanoic acid side chain is a typical transformation pathway for all statins, but is reported to play only a minor role in humans [75]. However, high activity of beta-oxidation of fatty acids, which is responsible for the supply of ATP for gluconeogenesis in type 2 diabetes mellitus [76, 77], which individual 3 was diagnosed with, may severely influence the AS metabolism. To test this hypothesis, respective reactions with the acid metabolites as substrates were considered in the model verification of individual 3, but the results showed no significant improvement (data not shown). The second reason could be, that CYP3A4 and UGT1A3 protein concentrations of individual 3 used in the estimation of corresponding r _max value via relative abundance approach (equation (11)) differ from the measured mean value (Table 4). Due to the fact that the r _max parameters of CYP3A4 and UGT1A3 show high sensitivities in the parameter sensitivity analysis (see Additional file 4 - Analysis of parameter sensitivity and following model reduction procedure), a variation in the values would have a high impact on the metabolic profiles. Therefore, in a further optimization step, the CYP3A4 protein concentrations are allowed to vary in the interval of mean ± standard deviation (Table 4), where the standard deviation of UGT1A3 protein concentration is assumed to be 30% of the mean value. Notably, an improvement in the model prediction on individual 3 could be achieved in the case of the intracellular metabolites AS and ASL and extracellular AS (Figure 5, dashed lines). But there are still major deviations in case of intracellular ASL and ASoOH and extracellular ASpOH, which could not be explained any further, but shows, that there must be some other reactions or effects in the system, which have not been discovered yet.

Figure 4

Measured concentrations and simulated fits of Atorvastatin compounds, individual 2. Displayed are the measured concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 2, and corresponding simulation profiles (solid lines). First and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH (right).

Figure 5

Measured concentrations and simulated fits of Atorvastatin compounds, individual 3. Displayed are the measured concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 3, and corresponding simulation profiles (solid lines and dashed lines). First and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH (right). The dashed lines show the additional consideration of putative beta-oxidation of atorvastatin acids, AS, ASpOH and ASoOH, and of adaptation of CYP3A4 and UGT1A3 protein concentrations in the model fit.

The estimated model parameters for individual 1, 2 and 3, summarized in Table 5, display the proposed parameter variability of enzyme mediated reactions as well as of the transport steps. However, the rate constants of ASLoOH import and export of individual 2 and of ASL import and export of individual 3 are considerably higher than the other transport rate constants, which points to strong linear dependencies of these parameters.

Table 5 Parameter Variability in the simultaneous model verification on individual atorvastatin metabolism data.

The mathematical model of AS metabolism in human hepatocytes of individual 1 is supplemented as SBML-file (see Additional file 7 - Model of atorvastatin metabolism in primary human hepatocytes of individual 1, and BioModels database).

Dynamic analysis of inter-individual CYP3A4 and UGT1A3 expression level variability

Based on the model version of optimized parameters, gained in the simultaneous model fit, the effect of inter-individual variability of CYP3A4 and UGT1A3 protein expression levels was investigated by linking the protein expression data of 150 liver samples (Figure 6) via the described relative abundance approach in equation (11), using individual 1 as reference. UGT1A3 and CYP3A4 protein concentrations of individual 1 were converted from based on total protein amount to based on microsomal protein amount by multiplying with the factor 0.22, determined in human liver homogenates and corresponding microsome fractions via Bradford test [78].

Figure 6

Individual variability of CYP3A4 and UGT1A3 protein concentration level. Displayed is the distribution of CYP3A4 (top) and UGT1A3 (bottom) protein concentration data of individual human liver microsomes (n = 150). Protein concentrations are normalized to minimal value, respectively.

Simulations were performed with an initial extracellular AS concentration of 50 (pmol ml^-1), which is in the range of physiological plasma concentrations [54, 79], over a time period of 1200 min.

The most important question that arises from this dynamic analysis was, how the metabolic profiles of the intracellular metabolites AS, ASpOH and ASoOH are influenced by this variability, since they are considered to be the active drugs, which inhibit HMGCoA-reductase [80]. Therefore, AUC, c _max and t(c _max ) of the concentration-time-profiles of either AS alone or the sum of concentration of AS, ASpOH and ASoOH were calculated for each liver sample over a time period of 1200 min and the distributions over all liver samples were evaluated, respectively. Finally, appropriate probability density functions are fitted to the distributions (Figure 7). The probability density function characteristics are summarized in Table 6. Obviously, there are differences between the examination of only AS or the sum of concentrations of AS, ASpOH and ASoOH. AUC, c _max and t(c _max ) have lower values in case of AS alone compared to the sum of all acidic metabolites. The population mean of AUC is 75051 (pmol ml^-1 min) in the case of AS and 203617 (pmol ml^-1 min) in the case of the sum of AS, ASpOH and ASoOH. Also, c_max is lower in the case of AS alone, 201 (pmol ml^-1), compared to the sum of the acidic metabolites, 366 (pmol ml^-1). Further, the maximal concentration appears at a shorter time point, 48 min, in the case of AS alone, compared to the time point, 100 min, of the sum of AS, ASpOH and ASoOH. The results are quite explainable, since ASpOH and ASoOH are the hydroxylated products of AS and therefore their maximal concentrations event at a delayed time-point compared to AS.

Figure 7

Distributions from dynamic analysis of the Inter-individual variability in Atorvastatin metabolism. Distributions (bars) of AUC (top), c_max (middle) and of t(c_max) (bottom), and fitted probability density functions of simulated profiles of AS (solid lines) and of the sum of AS, ASpOH and ASoOH (dashed lines), respectively. Distributions are predicted by the dynamic simulation of Atorvastatin metabolism, implementing individual CYP3A4 and UGT1A3 protein concentration data (Figure 6).

Table 6 Probability density functions fitted to inter-individual variability distributions

Further, except for the AUC of AS, the probability density functions show a very narrow shape, regarding the relative standard deviation, when comparing to the sample standard deviations of the underlying distribution of CYP3A4 and UGT1A3, respectively. The probability density functions fitted to CYP3A4 and UGT1A3 distributions have a relative standard deviation of 259% and 137%, respectively, whereas the relative standard deviations of the probability density functions of AUC, c_max and t(c_max) are lower than 50%, except for the probability density function of AUC of intracellular AS, which equals 123%.

Discussion

The deterministic modeling of drug metabolism has a major advantage compared to traditional pharmacokinetic models. The distinction between metabolism and elimination inside the hepatocytes and the detailed description of the biotransformation network structure allows the observation of the reaction kinetics and the transport mechanisms involved. Our modeling considered the major acid and lactone metabolites in the extracellular and intracellular space as well as appropriate reaction kinetics and transport steps. An experimental limitation in primary hepatocytes concerns the lower limit of quantification of each metabolite. Therefore, a relatively high initial concentration of AS of about 10 μM had to be chosen in order to permit measurement of the described major metabolites. Interestingly, the intracellular concentration of the parent drug AS was higher than the extracellular concentration. This indicates, that either the OATP mediated import does rather follow an active uptake mechanism than the proposed facilitated diffusion mechanism [58, 59], or so far unknown transporters may be involved in the transport steps. By any means it justifies the chosen implementation of kinetics for the active mediated transport beside the passive diffusion mechanism. In contrast to concentrations determined in plasma in clinical studies [54, 79, 81], extracellular concentrations of the acidic metabolites were much higher compared to the lactones. A possible explanation of these discrepancies may be that the in-vitro investigation on isolated primary human hepatocytes was performed over rather short time intervals (hours), whereas analyses in humans usually cover longer periods (days). A second explanation could be the relatively high initial concentration of AS of about 10 μM used in this study, compared to plasma concentrations observed to be lower than 200 (pmol ml^-1). Furthermore, we found that the recovery of metabolites decreases over the time. This confirms on the one hand the contribution of unspecific binding to macromolecules, most severely in the intracellular space, as observed previously [35], but on the other hand could also be contributed to a certain extent by the unspecific binding to the collagen layer or the plates [34], especially of highly lipophilic ASL and lactone metabolites. However, due to the washing procedure preliminary to the cell harvesting and disruption procedure, this effect could not be observed in the chosen experimental set-up.

The parameter identification was considered satisfactory, as the simulation profiles fitted well to the corresponding measured metabolite concentration data. But, it was questionable if the optimized parameters could be considered to be sensitive and identifiable. Therefore, a local parameter sensitivity and identifiability analysis based on the Fisher-Information-Matrix was performed. The results showed that the full model of AS metabolism is not identifiable. Consequently a model reduction procedure was set-up and could be applied successfully, indicating that the parameter identification difficulties are caused by non-identifiable transporter steps, displayed by low parameter sensitivities and high parameter errors and correlations, which were then reduced in the model. However, this outcome does not necessarily mean that the non-identifiable transport steps are not present or not used in human hepatocytes, but rather implies, that they cannot be distinguished from the remaining transport steps in the model verification. Thus, the remaining transport steps in the model capture the probable superimposed contribution of several transport mechanisms in-vivo.

Unfortunately, the local identifiability of the final model version reveals some uncertainties caused by remaining high correlations present in the transport steps, as well as in the intracellular reaction network. An additional problem attributed to the modeling and parameter optimization is related to the fact that several parameters have lumped characteristics, because they are used for more than one compound, like the fraction unbound factor, which is set to be the same for all acidic or lactone AS compounds, respectively. A possible solution would be to set individual parameters for each compound. However, then the parameter identifiability difficulties would be even more severe. This is due to the fact, that the experimental observation on primary human hepatocytes is constrained strongly on several limitations. For example, the range of initial concentrations is strongly limited to the quantification of metabolites of interest. Further, the number of data points is strongly limited to the available cell number of primary hepatocytes, isolated from surgery removal. Finally, also the individual character of human cells due to individual genetic and medical background precludes an exact reproduction of experiments and the reproducibility of experimental results. The remaining parameter uncertainty seems to be crucial for the prediction of the population probability on the first sight. Actually, this is not the case. Since the parameter errors and the correlations indicate which parameter changes in what extent do not influence the computed time courses, the variability in the respective parameter error range would not influence the probability distributions. In general, the question if parameter errors and insensitivities are sufficiently small can only be answered when knowing the requirements of the pharmacological application as drug, which was not given in this study.

The simultaneous model verification based on primary human hepatocytes from different individuals illustrates the individual character of drug metabolism in the liver. Further, it indicates the inter-individual variability of rate parameters not only of the phase I and II enzymes, but also of the transport enzymes. In case of patients with type 2 diabetes, AS metabolism is probably influenced strongly by beta-oxidation and further so far unknown effects. Therefore, the individual model verification might also be used to reveal undesired pharmacokinetic effects. However, this analysis tool should be checked carefully on other drug systems in the future, too.

The individual nature of the AS metabolism was investigated in this study in the domain of inter-individual variability of CYP3A4 and UGT1A3 protein expression levels in human hepatocytes, by performing dynamic analysis on the verified model, individualized by implementation of a comprehensive set of protein data from a liver bank. The results show, that the inter-individual variability strongly affects the biotransformation behavior and therefore reveals the individual character of AS metabolism. Consequently, this individuality in the pharmacokinetics of AS also points out the individual character of pharmacodynamics at the drug target, namely the HMGCoA-reductase in human hepatocytes.

However, this subject-variability can also be expected to be present in the transport protein expression, as previously shown for OATP-C [82]. Therefore, this variability should be taken into account in the future by implementing corresponding population expression data in the dynamic analysis. Another source of variation is the well-known polymorphism in CYP- and UGT- enzymes and likely transport proteins, which very probably causes inter-individual variations of the catalytic activity and the substrate affinity K _M . Thus, further efforts should also consider this individual difference, leading to the improvement of the model prediction.

To resolve the contributing transport mechanisms in more detail, additional information should be implemented from transport investigations on recombinant systems, as described previously [72, 83], which enable the differentiated identification of both basolateral and apical transporters of AS and the corresponding transport kinetics. Further, also the estimation of unspecific drug protein binding could be improved in the future by using radiolabeled compounds in the experiments and modeling approaches as described previously [60, 69].

Conclusions

We believe that the results of our simulations provide strong arguments for rigorous dynamic modeling of drug biotransformation at the cellular level embedded in a systems biology approach. By resolving the detailed metabolic network structure with metabolites and catalyzing enzymes, we investigated the dynamic variation of atorvastatin metabolism affected by the inter-individual variability of expression levels of phase I and phase II enzymes.

In contrast to experimental investigations on recombinant systems or tissue fractions of hepatocytes, like microsomes, the investigation on primary human hepatocytes enabled the holistic and most realistic in-vitro observation of drug biotransformation, because it is possible to observe the coupled contribution of metabolism and transport to the entire processes. De novo of this study, we identified intracellular concentration profiles of atorvastatin metabolites in primary human hepatocytes in a time-series approach.

Such an approach is essential for integration of further-reaching issues, such as drug-drug interactions, impact of regulation networks linked to nuclear receptors and particularly to quantitatively account for subject-variability. The integration of this variability caused by genetic or environmental variations is crucial for predictive pharmacokinetic modeling.

Such a rigorous modeling approach critically depends upon a tight link between experimental observations and model design, simulation and verification. While the results are promising, some limitations in the parameter identification were still encountered. In the long term, these open problems can only be solved by stronger links to other research areas, such as pharmacogenetics, characterization of transporters, etc. On the whole, our contention is that the problem of parameter identifiability is an indispensable ingredient of model verification. Systematic investigations - if possible linked to optimal experimental design - can greatly strengthen the credibility of the models.

However, we present a model that goes much further. The domain of application does not remain the system behavior for which it has been elaborated. The model provides the possibility to link further modules such as gene regulation, drug target metabolism and present the important links to be implemented into the PBPK environment. Finally, the model structure used in this study should be considered as a module to be integrated into the framework of multi-scale whole body modeling and simulations necessary to tackle the drug disposition in patient populations.