A Liver-Centric Multiscale Modeling Framework for Xenobiotics

We describe a multi-scale, liver-centric in silico modeling framework for acetaminophen pharmacology and metabolism. We focus on a computational model to characterize whole body uptake and clearance, liver transport and phase I and phase II metabolism. We do this by incorporating sub-models that span three scales; Physiologically Based Pharmacokinetic (PBPK) modeling of acetaminophen uptake and distribution at the whole body level, cell and blood flow modeling at the tissue/organ level and metabolism at the sub-cellular level. We have used standard modeling modalities at each of the three scales. In particular, we have used the Systems Biology Markup Language (SBML) to create both the whole-body and sub-cellular scales. Our modeling approach allows us to run the individual sub-models separately and allows us to easily exchange models at a particular scale without the need to extensively rework the sub-models at other scales. In addition, the use of SBML greatly facilitates the inclusion of biological annotations directly in the model code. The model was calibrated using human in vivo data for acetaminophen and its sulfate and glucuronate metabolites. We then carried out extensive parameter sensitivity studies including the pairwise interaction of parameters. We also simulated population variation of exposure and sensitivity to acetaminophen. Our modeling framework can be extended to the prediction of liver toxicity following acetaminophen overdose, or used as a general purpose pharmacokinetic model for xenobiotics.


Biological Components and Processes
A list of the biological components and behaviors (processes), along with ontological references, are given in Table 1. The liver is composed of millions of hepatic lobules. Liver lobules are "plumbed" in parallel and a portion of blood, while transiting the liver once, passes through exactly one lobule. Each hepatic lobule consists of various cell types including hepatocytes, Kupffer cells, stellate cells and endothelial cells and each of these contribute to normal liver function. In our simplified tissue model we focus on the liver's parenchymal cells (hepatocytes) and the two main components of blood (red blood cells and blood plasma). In addition, we have significantly simplified the sinusoidal network of the lobule by assuming a two dimensional, tubular structure lined with hepatocytes. Thus, we only modeled a small section of vascular-tissue volume and extend the result to whole liver by simple multiplication. Coupling of the Hepatic Lobule model to the whole-body PBPK model is done via the concentration of APAP (and metabolites) in the blood as the blood moves from the higher scale (PBPK) to the lower scale (Lobule) and back again. Since concentration is an intensive parameter it is scale insensitive.

Diffusion and Transport of Acetaminophen
In the lobule model we modeled the diffusion and transport using a coarse-grained PDE method. The cells (and pseudo-cell serum portions) are the coarse-graining scale and diffusion and transport is modeled only at that scale. Individual cells are treated as "well stirred" and diffusion within a cell is not explicitly modeled. In the whole-body PBPK model the individual tissue compartments are treated as "well stirred" compartments.

Protein Binding of APAP and Metabolites
We include protein binding of APAP in blood serum characterized by the Fup in the PBPK module. In our model, APAP-F up applies only to transport between blood serum portions and hepatocytes and between blood serum portions and the kideny tubules. We neglect protein binding of APAP within RBCs and hepatocytes. We also assumed that the two APAP metabolites (APAPG and APAPS) do not bind to serum proteins and are not taken up by RBCs.

Subcellular Metabolic Reactions
We include representative common Phase I and Phase II metabolic reactions. For APAP these reactions represent the vast majority of the in vivo metabolism. Many other xenobiotics undergo similar reactions and the sub-cellular metabolic model is at least a starting point for modeling a wide range of chemical entities. The single Phase I reaction modeled represent Cytochrome P450 oxidation of APAP by a pool of enzymes (e.g., Cyp2E1 and Cyp1A2) to give NAPQI. NAPQI reacts rapidly with cellular nucleophiles such as glutathione. The high concentration of cellular glutathione in hepatocytes is generally sufficient to scavenge the NAPQI formed by therapeutic doses of APAP. Phase II metabolism of APAP is known to be quite extensive with both sulfate and glucuronide conjugates being formed. We model these two separately as first order reactions between APAP and unlimited pools of sulfate and glucuronic acid precursors. We are aware that APAP overdose-induced liver necrosis is likely to involve other processes, as well as saturation of the the processes we included, leading to tissue death. However, for the study of therapeutic dose of APAP, we assumed that necrosis-related pathways could be omitted. (1) Parameters that are only dependent on the individual (species, gender, body weight...) being simulated and independent of the compound and (2) Parameters that are only dependent on the compound but independent of the individual. Table 2 lists the allometric scaling factors used to calculate the tissue compartment volumes and volumetric blood flow rates. The scaling equation uses the ratio of the body weight (BW ) (in Kg) to a reference body weight (BW ref ) raised to the 3/4 power for parameters such as cardiac output and Qgfr;

Allometric Scaling of the PBPK model's Compartment Volumes and Blood Flow Rates
Where P i,ref is a standard value at the standard BW . For parameters that don't scale as the 3/4 power, such as the total volume of blood, the scaling equation is; The scaling factors and fraction values used were for humans and are given in Table 2.
The PBPK compartments are of three types; "Circulatory" (arterial and venous compartments), "Organ" (liver, kidney, ...) and "Lumen" (gut lumen, kidney lumen). The transfer equations for movement of compound are generally based on the compartment type and developed based on the transfer out of a compartment. See the main paper for the complete set of transfer ODEs.

PBPK ODEs
The complete set of ODEs for the PBPK model are given in Table 3. In addition, the ODEs can be extracted from the SBML model file or exported from the SBML model file into R, Matlab, etc. code using standard SBML tools (cf. http://sbml.org/Software/SBMLToolbox).

Intercellular Molecular Transfer Model
The molecular transfer model is coded and run in CC3D.

Intercellular Molecular Transfer Rate Equations ODEs
The complete set of ODEs for the Molecular Transfer model for APAP, APAPG and APAPS are given in Table 4,5,6. Superscripts refer to the unique identifiers for each cell. The set of cells C j are the contact neighbors of the cell C i . Subscripts identify the cell type: H = hepatocyte, S = serum portion, R = red blood cell. pbpk F up = Fraction unbound to protein, n i = number of contact neighbors.

Subcellular Reaction Kinetics Model
The metabolic reaction kinetics model is coded and run in SBML in the same way as the PBPK model.

Subcellular Reaction Kinetics ODEs
The complete set of ODEs for the reaction kinetic model are given in Table 7. In addition, the ODEs can be extracted from the SBML model file or exported from the SBML model file into R, Matlab, etc. code using standard SBML tools (cf. http://sbml.org/Software/SBMLToolbox).

Cellular Potts Model in Compucell3D
Cellular Potts Model (CPM) defines a generalized cell as a collection of lattice sites or voxels at location l, acquiring common cell index, σ, in a cell-lattice. Generalized cells refer to either biological cells (hepatocytes, red blood cells) or pseudo cells (blood plasma, blood source cells). Each cell has a cell type Version: March 10, 2016 associated with it. Evolution of cell-lattice configuration in the system is governed by Effective Energy, which is constructed in this model using factors introduced below.

Cell Contact, Volume Constraint
The basic expression of effective energy (H) in CPM has terms contributed by adhesion and volume constraint: The summation of J(τ (σ i ), τ (σ j )) calculates the contact energy between two neighboring cells. J(τ (σ i ), τ (σ j )) represents energy per contact area between the cell σ i with cell type τ (σ i ) located at cell lattice site i and the cell σ j with cell type τ (σ j ) located at cell lattice site j. And calculation excludes the pair of lattice sites belonging to the same cell by multiplication with a factor (1 − δ(σ i , σ j )). The second summation over all generalized cells σ defines a penalty energy for cells with a volume deviated from its target volume V t (σ), and the strength of penalty is described by λ vol (σ).

External Potential
The external potential is incorporated by imposing a directed force on the generalized cells. It's used to achieve persistent cell movement. The additional term in the effective energy (H) caused by the external force is: The summation calculates the contribution of external force to the total effective energy in the system. The term − λ ext (σ( i)) is the force vector for cell σ( i) at cell lattice site i. c is the direction of the pixel copy attempt.

Fields and Flux
We used a course-grained method to model the diffusion processes in the multi-cell sinusoid model. Chemicals are considered uniformly distributed within each generalized cell. Concentration gradients only exists across the interface between different instances of generalized cells. For passive transport, flux is driven by the concentration gradient across the common surface of two cells. For active transport, flux is dependent on the concentration within the source cell.
APAP APAP exists in all generalized cells except BSC (blood source cells). Transport flux of APAP falls into two categories; Passive transport of APAP is modeled between RBC, SP and HEP. Active transport of APAP is modeled for the uptake by HEP (SP → HEP ). Passive transport processes is modeled as a linear ODE, active transport is modeled assuming Michaelis-Menton kinetics. Reaction of APAP to form Phase I and Phase II metabolites are modeled, within each HEP, as Michaelis-Menton kinetics.
APAPG and APAPS We assume that APAPG and APAPS exist only in serum portions (SPs) and this is enforced in the PBPK model by the blood-to-plasma partition coefficient. APAPG and APAPS are assumed to not be taken up by other tissues except the kidney.
GSH, NAPQI and NAPQIGSH We assume that GSH, NAPQI and NAPQIGSH exist only in HEPs and do not transfer between HEPs. GSH participates in a conjugation reaction with NAPQI (to form NAPQIGSH), and GSH is replenished by a relatively slow reaction from unlimited, and unspecified, precursors.

Simulation Parameters and Initial Configuration
In our multiscale model, simulation parameters span several scales. At whole-body level, parameters include blood volumetric flow rate Qs, organ perfusable volume V s, body weight bw, hematocrit hemat, partition coefficient K s and Rs; At sub-cellular level, parameters are mainly reaction rate coefficients Vmax s, Kms and k s. At the multi-cell level, parameters can be divided into chemical-dependent and chemical-independent sub-groups. Chemical-independent parameters include sinusoid geometry, blood flow rate, blood composition and cellular qualities such as volume and adhesion energies between the various cell and pseudo-cell types. Chemical-dependent parameters include diffusion and transport rate constants (or partition coefficients). CC3D parameters include adhesion energy terms J(τ (σ i ), τ (σ j )), volumetric constraint strength λ vol , external potential strength λ ext . Transport parameters are mainly passive transport rate coefficients k PDs, active transport Michaelis-Menton constants Vmax s, Kms. These parameters were estimated by fitting to invivo human data for therapeutic doses of APAP. Geometrical parameters are determined according to the simplified model architecture (see below).

Geometrical Parameters
Our model focuses on a liver lobule section including one sinusoid vessel and two rows of hepatocytes in a 2D plane. The whole lobule was treated as a collection of repeats of this minimal structure. The entire liver is treated as a collection of lobules. The model has a 2D interface but a pseudo-3D setup with thickness allowing estimation of 3D diffusion calculations (e.g., subcellular RK input). In the multi-cell model, the relationship between model pixels and simulated distances is given by the scaling factor px2mm. We assume that each sinusoid vessel has a length of 200µm (approximately 10 hepatocytes), connecting a portal triad to a central vein, and a cross sectional area of 20µm × 4µm. We assume the hepatocytes are cubical with edges of 20µm.

Temporal Parameters
Simulated time is a key component of all three of the sub-models and proper synchronization of time across the sub-models is critical. We take one CC3D "Monte Carlo Step" (MCS) as the minimal temporal length in our model. Every MCS the multicell model is updated. In addition, the CC3D model handles the transport as an innate part of the multicell module and are integrated in CC3D python code. Every 10 MCS, the ODE-based modules are integrated and updated.

Initial Configuration and Condition
Each simulation starts with regularly aligned blood source cells (4 in total) and blood serum portions (4 in each stack) filling in the sinusoid vessel. Hepatocytes are located laterally in contact with the vessel walls (10 for each edge). The key input of the model is the APAP oral dose (in grams) specified in the CC3D script. The dose value is copied to the PBPK module at the start of the run.

Parameter Searching
Parameter searching was carried out by choosing multipliers for each of the 38 parameters in the complete model. The multipliers were chosen from a Log-normal distribution with a mean of 0 and a coefficient of variation of 25% with the restrictions outlined below. The parameter assignment rules for the parameter searching modeling were: Where α is selected from a normal distribution with mean of 1 and standard deviation of 0.25. The choice of α is further constrained by; Constraint C 1 : −2 < α < 2 (P old /1000 < P new < P old · 1000) Constraint C 2 : −0.1 < α < 0.1 (P old /1.26 < P new < P old · 1.26) The parameter scaling and constraints were applied to: Parameters

Population Variability
Population analysis was carried out by choosing multipliers for each of the 38 parameters in the complete model. The multipliers were chosen from a normal distribution with a mean of 1 and a coefficient of variation of 25% with the restrictions outlined below. We used 839 randomly generated parameter sets to represent a population of individuals. The parameter assignment rules for the population variability modeling were: Where β is selected from a normal distribution with mean of 1 and standard deviation of 0.25. The choice of β is further constrained by; Constraint C 1 : 0.75 < β < 1.25 The parameter scaling and constraints were applied to: Parameters not changed : pbpk dose Parameters obeying C 2 : pbpk FupG, and pbpk FupS Parameters obeying C 1 : all other parameters 5 Simulation Code Files

Physiologically based pharmacokinetics module
Refer to file vLiver PBPK based v4 final.xml. Standalone simulations can be done using Jarnac or COPASI. COPASI will complain about the unit definitions, since COPASI doesn't fully support the SBML unit standards, but the model runs correctly. Since a number of the model parameters (compartment volumes and volumetric flows in particular) are assigned via SBML IntialAssignment statements those values will be displayed incorrectly in many user interfaces. When the model is run the correct parameters will be used.

Multi-cell flow and transport module
Refer to the sub-directory "cc3dCode" for reference simulation files. XML files in this directory are not annotated. See the Compucell3D website for tutorials if necessary (http://Compucell3D.org/).

Parameter scanning scripts
Refer to sub-directories log normal, single parameter, pairwise parameter, and population variability for scripts to generate multiple simulation files based on a particular reference simulation parameter set, and relevant submission scripts.

Model Parameters
Detailed lists of parameters for each of the sub-models is given in the tables that follow. Energy and volume parameters for the multi-cell lobule model in CC3D are given in Table 8 and the spatiotemporal parameters in Table 9.

Named Model Parameter Sets
The Table 10 includes a listing of the parameter sets (fixed points) described in the paper. The parameter sets REFSIM and HMPCsim6 both give reasonably good agreement with in vivo human data for a 1.4g oral dose of APAP. Paramter sets LNsim8 and LNsim23 were purposely chosen to represent compounds with ADME characteristics significantly different from APAP.

Model Outputs
The complete list of model outputs for each named parameter set is given in table 11.

Model Reproducibility
While the PBPK and sub-cellular reaction kinetic ODE models are deterministic, the multi-cell lobule module is stochastic, hence the entire model is stochastic. To examine the run-to-run variability caused by the stochasticity we ran 50 replicates of the complete model using the same set of initial conditions and the REFSIM parameter set. As shown in Table 12, the standard deviations for all of the model outputs are parts per billion or less, indicating no significant run-to-run uncertainty or variability caused by the stochasticity of the lobule sub-model. † The resolution of time outputs is limited by the frequency that the multi-scale model writes data to a log file. For these 50 simulations the time resolution in the log file is approximately 5 minutes and any variation of less than that amount is not visible. ‡ See Table 11 for unit definitions.

Sensitivity of model output to model parameters
The sensitivity of the model's output to the models input is given in the following two tables. Table 13 gives the sensitivity of the single model output RMSEsum to the parameter values in the four named parameter sets. Table 14 gives the average sensitivity across all of the model outputs (listed in the first column in Tables 11 or 12) to the parameter values in the four named parameter sets.