In-vitro to in-vivo acetaminophen hepatotoxicity extrapolation using classical schemes, pharmaco-dynamic models and a multiscale spatial-temporal liver twin

In vitro to in vivo extrapolation represents a critical challenge in toxicology. In this paper we explore extrapolation strategies for acetaminophen (APAP) based on mechanistic models, comparing classical homogeneous compartment pharmaco-dynamic (PD) models and a multiscale digital twin model resolving liver microarchitecture at cellular resolution. The models integrate consensus detoxification reactions in each individual hepatocyte. We study the consequences of the two model types on the extrapolation and show in which cases these models perform better than the classical extrapolation strategy that is based either on the maximal drug concentration (Cmax) or the area under the pharmaco-kinetic curve (AUC) of the drug blood concentration.


INTRODUCTION
: In vitro -in vitro extrapolation approaches. In vitro, toxicity is measured by exposing cell populations of interest to certain concentrations Cvitro of drug (here: APAP) and measuring the fraction of cell death with time (I). The toxicity at a certain time T (often: 24h, or 3d) is considered as the time at which the toxicity is compared (II). A common consensus for extrapolation to in vivo is that the in vitro concentration value Cvitro is identified by the Cmax-value (III) of the corresponding pharmaco-kinetic (PK) curve for the drug in the blood (IV). The PK curve depends on the administered drug dose, so that by identifying Cmax and Cvitro, the dose value associated with that Cmax-value can be associated with the toxicity value associated with the corresponding Cvitro value (V). (Alternatively, to the Cvitro / Cmax values, the area under the curve (AUC) is used). The model-based strategy mimics drug toxicity in a model by simulating the process of drug detoxification in vitro by the cells, whereby the toxicity pathway is explicitly represented in each cell and eventually integrated with a PK model (IV) and a compartment model of the organs of interest (VI). The latter may represent organ microarchitecture, here the liver lobule (VI-C) as repetitive minimal tissue unit or consider a well-mixed compartment (VI-B), both integrating an intracellular PD model (VI-A). Simulations with (VI-B, C) directly yield the in vivo toxicity (VII).

3
As direct risk-assessment of chemicals in human is, of course, prohibited, risk-evaluation usually is based on animal experiments. For cosmetics, animal experiments are forbidden in EU since 03/2013 1 . Large ongoing research programs focus on the ultimate goal to deliver testing strategies that enable animal-free risk assessment of chemicals (Godoy et al., 2013;Albrecht et al., 2019;Cherianidou et al., 2022). However, despite of large efforts in quantitative risk evaluation by in vitro and/or in silico methods, identification of no-observed-adverse-effect levels (NOAEL) or lowest-observed-adverseeffect levels (LOAEL) by alternative methods is presently not possible with satisfactory accuracy (Leist et al., 2017). Therefore, replacement of animal experiments has so far not been possible in most fields of regulatory toxicology. However, long-term, extrapolation of adverse drug effects to humans from in vitro experiments constitutes a key objective (Brecklinghaus et al., 2022).
In vitro testing usually begins with the establishment of a concentration-response relationship in cell cultures or in more complex 3D culture systems ( Fig. 1(I)). Typical responses studied in relation to concentration of a test compound are cytotoxicity or biomarkers linked to adverse effects. A next challenge is in vitro to in in vivo extrapolation (Sachinidis et al., 2019). Here, an in vitro concentration-response curve is usually translated to a dose-response curve, e.g. of organ toxicity, in vivo (Fig 1(V)). The easiest way of in vitro -in vivo extrapolation is to compare toxicity at in vivo relevant concentration ranges. An accepted procedure is the quantitative comparison based on the plasma peak concentration (Cmax ; Cmax) or the area under the curve (AUC) of a test compound (Fig. 1(IV), Fig. 1(III)); named "direct Physiologically-Based-Pharmacokinetic (PBPK) based extrapolation'' from now on. This procedure often fails because: (1) In vitro, concentrations of test compounds in the culture medium are either relatively constant during the exposure period or decrease only slowly, due to the much larger volume of culture medium compared to the cell fraction. In contrast, faster concentrations changes occur in vivo, where compounds are eliminated by zero or first order kinetics unless the compound is infused or delivered over a long period of time, which is impractical for many drugs. (2) In vitro / in vivo-differences in the enzyme activities; and (3) the spatial organization of cells impacting the spatial-temporal transport and detoxification profiles.
In this paper, we studied several computational model (CM)-based extrapolation approaches from in vitro to in vivo hepatotoxicity of APAP, that gradually take each of these three factors into account.
APAP is a well-suited, clinically relevant toxic model substance to study in vitroin vivo extrapolation strategies. Relevant, as overdosing of APAP is the main reason for acute liver failure in several countries (Olson et al., 2017). An APAP overdose is cytotoxic for hepatocytes expressing cytochrome P450 (CYP) enzymes; mainly Cyp2e1, Cyp1a2 and Cyp3a4 (McGill and Jaeschke, 2013). These are localized close to the central vein of the liver lobules (Ghallab et al., 2019b), which constitute the repetitive anatomical and functional units of the liver. Liver lobules have a complex architecture facilitating exchange of metabolites between blood and hepatocytes, the parenchyma of the liver. Some chemical reactions are only executed in specific spatial regions of a lobule, a phenomenon called zonation (Gebhardt, 1992;Bartl et al., 2015). In animal models, here mice, it is possible to dose-dependently determine the fraction of dead hepatocytes induced by a test compound in vivo; this can be done under conditions where blood is sampled time-dependently from relevant sites, e.g. the portal vein as liver inflow, the liver vein as liver outflow or mixed venous blood, as internal exposure monitoring (Sezgin et al., 2018). Moreover, hepatocytes can be isolated from the same mouse strain and incubated concentration-dependently in vitro (Ghallab et al., 2022).
Here we confront a classical PB/PK/PD model (CL-model) approach, assuming a well-mixed compartment composed of CYP-positive-and CYP-negative hepatocytes (indicated by the yellow/blue box in Fig. 1(VI-B)), and a spatial-temporal PB/PK/PD-model (ST-model) resolving liver microarchitecture (indicated by the lobule in Fig. 1(VI-C)), with experimental data on APAPbased hepatotoxicity in mice. The well "mixed-mixed" model approximates the population of hepatocytes as a stirred container filled with cells. The ST-model represents each individual hepatocyte as basic model unit inside a realistic sinusoidal network (Hoehme et al., 2010;Drasdo et al., 2014) representing a digital twin model of APAP toxicity in liver microarchitecture, hence naturally capturing the effect of liver zonation (which is not the case for the well-mixed model). The APAP detoxification model is executed inside each individual hepatocyte in both, the well-mixed and the ST-model ( Fig. 1(VI-A)). We tested several model variants for the CL-model and different model-based extrapolation strategies (summarized in SFig. 1). Many a-priori plausible strategies did not function. The final strategy provided us with a parameter set that gives a valid prediction. However, the parameter identification procedure in that strategy (explained below) provides severaldifferent -parameter sets leading to almost the same agreement to the experimental data sets that were used for the parameter identification procedure, but not all parameter sets lead to an in vivo toxicity prediction within the experimental uncertainty. The failures are informative and indicate that limitation of parameter ranges by proper experiments is indispensable for reproducibility of a model-based extrapolation. For successful parameter sets the extrapolation outperforms the classical extrapolation strategy based on Cmax and AUC, which we performed as reference ( Fig. 1(IV)). The difference between the CL-model and ST-model for our simulations was small, however, the advantage of the ST-model is that it permits to directly include architectural distortions as for example fibrotic zones in liver fibrosis and cirrhosis, which is not possible in the CL-model.
The successful extrapolation simulation pursued in this work was composed of the following steps ( (1) As input data for the calibration of the in-vivo model (CL-model or ST-model) hepatotoxicity data from the in vitro measurements ( Fig. 1(I), Fig. 2D-F), and the APAP pharmaco-kinetics ( Fig.1(IV'), Fig. 2I) were determined. More specifically, this data serves to calibrate the model parameters of a model that composed of an intracellular PD model coupled to a supra-cellular PK model composed of different body compartments ( Fig. 3C.1).
(2) An intracellular PD model ( Fig. 1(VI-A)) and an organ compartment (physiologically-based) PB/PK model were set up, whereby for the latter the liver was represented by two distinct approaches: (3a) A well-mixed liver model whereby all cells were considered as independently fed by the APAP blood concentration ( Fig. 1(VI-B); Fig. 3C.1). The model integrating the intracellular PD-model into that well-mixed liver model is a classical (CL-)model (denoted as CL-model 3 in Fig. 3).
(3b) A spatial-temporal (ST-) resolved liver micro-architectural multiscale model, where the intracellular PD model was executed in each individual cell of a virtual lobule ( Fig. 1(VI-C)), referred to as ST-model ( Fig. 3D.1).
(4) Within a single fit procedure ( Fig. 3C.1), the parameters of the intracellular PD model and the parameters of the PB/PK model (together forming CL-model 3, cf. Fig. 3C.1) were fitted together subject to the following conditions: (i) the parameters of the intracellular PD models were fitted to the in vitro data ( Fig. 1(I)). (ii) For the simultaneous fit of the (coupled) PD/PB/PK model (composed of intracellular PD-model and the well-mixed PB/PK model, cf. Fig. 3C.1) to the PK data ( Fig. 1(IV')), the enzyme activities for CYP2E1 and CYP1A2 and the GSH levels obtained from the fit to the in vitro data (i) were corrected by the ratios found from an experimental determination of enzyme activities from in vivo hepatocytes. I.e., despite the fits to the in vitro toxicity data and to the in vivo PK-data were executed simultaneously, the parameters of the intracellular PD-model used in both sub-steps (in vitro / in vivo) were not the same but differed by the levels of the aforementioned three factors (CYP2E1, CYP1A2, GSH). However, the model structure of the intracellular PD-model as well as the other parameters were the same.
Once a suitable fit of the parameters of CL-model 3 was obtained, the so parameterized CL-model 3 composed of parameter sets {intracell-PD, PB/PK, cf. Fig. 3C.1} was executed to predict the in vivo toxicity ( Fig. 3C.2), and the simulation results were compared to experimental data ( Fig. 1(VII)).
(5) The same parameter set was then chosen to execute the spatial temporal (ST-)model ( Fig. 3D.1), whereby the micro-architecturally resolved model required one additional parameter, the liver flow rate. We found that the CL-model and the ST-model performed similarly for the same parameters used. The extrapolated in vivo toxicity results partially performed very well. However, for the fit of the intracellular PD-model, different parameter combinations performed equally well with regard to the fitness criterion, but they did not with regard to the in vivo toxicity prediction. This indicates that the parameter landscape is rough so does not guarantee convergence to the same set of intracellular model parameters.
From the simulation results it would be expected that experiments narrowing the parameter ranges of the intracellular model are likely to result reproducibly in valid in vivo extrapolations that outperform the classical extrapolation schemes based on Cmax and AUC.
Below, first an overview of the tested strategies and models is given, explaining which of the strategies failed, and in which way they failed, before finally explaining the model that succeeded. The success case confirms that a CL-model as well as a ST-model is a-priori capable to simultaneously explain in vitro and in vivo hepatotoxicity, so the structure is a priori suited. However, validation of the model would require to further constrain its parameters, which was out of the feasibility of this work. The ST-model is presented here as a proof of a concept of a spatially resolved multi-scale model of in vivo hepatotoxicity. Such models are expected to capture key components of hepatotoxicity (Holzhütter et al., 2012;Schwen et al., 2016;Ho and Zhang, 2020).

Animal experiments
Eight-ten-week-old male C57Bl/6N mice, 20-25 grams body weight were used (Janvier labs, France). The mice were fed ad libitum with the Ssniff R/M-H, 10 mm standard diet (Ssniff, Soest, Germany) and housed at controlled ambient temperature of 25• C with 12 h day, 12 h night cycles. All experiments were approved by the local animal ethics committee.

Induction of liver injury by APAP.
In order to test the dose response of APAP in livers of mice, various doses (89 up to 500 mg/kg) were administered, where APAP was dissolved in warm PBS and injected intraperitoneally (i.p.) into overnight fasted mice, with an application volume of 20 ml/kg (Schneider et al., 2021b). On day one after APAP injection, liver tissue samples were collected and processed for histopathology and immunohistochemistry analyses as previously described (Campos et al., 2020;Holland et al., 2022) ). Three mice were used for each tested dose. All experiments with mice were approved by the local authorities.

Pharmacokinetic analysis of acetaminophen.
For pharmacokinetic analysis of APAP, a dose of 450 mg/kg APAP was administered i.p.. Blood samples were collected in a time-dependent manner after APAP injection (0,5,15,30,45,60,120,240 and 480 minutes) from the portal vein 'representing 75% of the liver inflow', from the heart 'representing 25% of the liver inflow' and from the hepatic vein 'representing liver outflow' as described in Ghallab et al. (2016). After blood collection plasma was immediately separated by centrifugation at 13,000 rpm for 10 minutes and stored at -80°C until analysis. APAP concentrations were determined by LC-MS/MS as previously described (Sezgin et al. 2018).

Hematoxylin and eosin staining
Hematoxylin and eosin (H&E) staining was performed in 5 µm-thick formalin-fixed paraffinembedded liver tissue sections as described in Ghallab et al. (2021). Representative images were then acquired with a bright field microscope using Cell˄F software (Olympus, Hamburg, Germany) and from these the dead cell areas determined.

CYP2E1 and CYP1A immunostaining
Immunostaining of CYP2E1 and CYP1A were performed in 5 µm-thick cryo-sections using antibodies against CYP2E1 (Cat. No. MFO-100, Stressgen, Victoria, BC, Canada) and CYP1A (a gift from Dr. R. Wolf, Biochemical Research Centre, University of Dundee, Dundee, UK). In order to detect antibody binding, the tissue sections were incubated with a horseradish peroxidaseconjugated secondary antibody (Cat. No. P0217, DakoCytomation Denmark A/S, Glostrop, Denmark). Following washing steps, antibody binding was visualized by covering the tissue sections with AEC+ high sensitivity substrate chromogen (Dako, USA) for 10 minutes. The AEC-stained sections were then preserved by mounting with an aqueous mounting media (Schenk et al., 2017). Quantification of the CYP positive area was done in whole slide scans as previously described (Ghallab et al., 2019b) (see also SI section 2.10).

Measurement of CYP450 enzyme activity in isolated liver microsomes
Liver microsomes were prepared as described before (Lang et al., 2001). Microsomal incubation mixtures contained 50 μg of microsomal protein in 0.1 M sodium phosphate buffer (pH 7.4), NADPH-generating system (5 mM MgCl2, 4 mM glucose 6-phosphate, 0.5 mM NADP+, and 4.0 U/ml glucose 6-phosphate dehydrogenase), and 50 µM chlorzoxazone (CYP2E1 substrate) and phenacetin (CYP1A2 substrate) in a final volume of 100 μl. Samples were preincubated in a water bath for 5 min at 37 ºC and the reaction started by addition of the NADPH-generating system. Reactions were stopped after 15 min by adding 10 μl of 250 mM formic acid and 10 µl of internal standard solution containing a mixture of deuterium labelled analogues of the analytes and cooling on ice. The supernatant was analyzed by LC-MS/MS as described previously (Feidt et al., 2010).

Isolation and cultivation of primary mouse hepatocytes
Primary mouse hepatocytes were isolated from the livers of male C57Bl/6N mice, 8-10-week-old, (Janvier labs, France) as previously described (Godoy et al., 2016). Briefly, the liver was perfused through the vena cava for 15 minutes with EGTA buffer at 37°C. Subsequently, the liver was perfused with collagenase buffer for 5-7 minutes. The liver was then excised and dissociated in a suspension buffer. Following a filtration step, the cell suspension was centrifuged for 5 minutes at 50 g. The cells were re-suspended in 10 ml suspension buffer. The cell viability was checked by trypan blue exclusion. In order to prepare confluent hepatocyte culture, the cells were seeded in six-well plates at a density of 800,000 cells per well in Williams E medium supplemented with 10% fetal calf serum (FCS), 2mM L-glutamine, 100 units/ml penicillin, 0.1 mg/ml streptomycin, 10 μg/ml gentamycin (PAN Biotech, Aidenbach, Germany) and 100nM dexamethasone (Sigma-Aldrich, Munich, Germany). The cells were allowed to attach for 2h at 37°C, 5% CO2. Subsequently, the floating cells were washed out. In order to prepare sandwich cultures a second layer of collagen-1 was added and allowed to polymerize for 30 minutes. For all subsequent cultivations normal media was used without FCS.

Measurement of CYP450 enzymes activity in cultivated hepatocytes
In order to measure the activities of CYP2E1 and CYP1A2 in isolated cells, primary mouse hepatocytes were incubated with 50 µM chlorzoxazone and phenacetin directly after isolation (fresh hepatocytes) or at different time intervals (0h, 2h, 4h, 1 day, 3 days and 7 days) after attachment of the cells. After 60 and 120 minutes of incubation, 50 µl of the cell culture supernatant was taken and the reaction was stopped by adding 5 µl of formic acid (250 mM). After mixing, the samples were stored at -20 °C until analysis. After thawing samples were spiked with internal standard mixture and analyzed by LC-MS/MS as described previously (Feidt et al., 2010).

Cytotoxicity of acetaminophen in vitro
To check APAP-induced cytotoxicity in cultivated primary mouse hepatocytes (sandwich cultures), the cells were incubated with various concentrations of APAP (0, 0.25, 0.5, 1, 2, 4, 8, 16 and 32 mM). Incubations periods started 2h after plating of hepatocytes. This early exposure was chosen because primary mouse hepatocytes are known to show a particularly fast decrease of expression of phase I and phase II metabolizing enzymes even when cultivated in organotypic 3D cultures (Godoy et al., 2016). 24 h later, the cells were incubated with propidium iodide (PI) 1.5 mM (Thermo Scientific, MA, USA) diluted 1:500 in cultivation media for 5 minutes at 37°C. Moreover, 3D cultivated mouse hepatocytes were incubated for 5 min or 3 h with APAP followed by washout and analysis of cell morphology and PI uptake at 24 h to understand the influence of the exposure period to APAP (SFig 1A, B). Already incubation for only 3h led to an increase in PI-positive hepatocytes compared to controls for all tested APAP concentrations (SFig 1C). PI positive nuclei were visualized using a combined fluorescence/phase contrast microscope (Nicon, Dusseldorf, Germany). The number of dead cells was manually quantified based on PI positive nuclei and cell morphology (SFig. 1B). The representative result of PI-stained hepatocytes demonstrates that the technique allows a clear differentiation between PI positive and negative hepatocytes ( Fig 2B). Five representative images were quantified per well. The results are from three independent experiments (with different mice), and each experiment was performed in triplicate.
A second set of experiments was performed for the 4mM APAP-concentration at 5 different time points: 1, 4, 8, 18, 24 hours. This data set was rescaled to match the 24h data set (see SI Fig. 1). The scaling was necessary as the 24h-toxicity and the dead-cell kinetics were measured using hepatocytes from two different mice, each mouse having different levels of GSH prior to starvation that were not known.

Measurement of GSH concentration
The glutathione content in liver tissue homogenate was measured using LC-MS/MS as previously described (Sezgin et al., 2018).

Mathematical modeling
Pharmacokinetic (PK) model The PK model considers two compartments: the peritoneal cavity into which APAP is injected, and the blood compartment, into which APAP diffuses from the peritoneal cavity. Clearance of APAP in the blood is modeled by a first-order decay term. Only a fraction of the dose injected intraperitoneally is assumed to reach the blood compartment according to a limited bio-availability. Parameters were fitted for the experimentally measured blood kinetics observed for the dose of 300 mg/kg of body weight and one data point measured at 24h each for 56 mg/kg, 167mg/kg, and 500mg/kg through minimization of the log-likelihood function using the CMA-ES algorithm (Hansen and Ostermeier, 1996;Hansen, 2006). In a next step, prediction of the concentration-time curve for the doses 89, 158, 281, 375 and 500 mg/kg was done in MATLAB using the estimated parameters. Equations are presented in the supplementary material. This model was used to perform the AUC and Cmax strategies, and for CL-model 1 (SFig. 1-A.2).

In vitro-PD compartment model of APAP
The in vitro metabolism of APAP was described using an ODE-based pharmacodynamic model. In brief, APAP diffuses passively through the cell membrane as this constitutes the main transport mechanism (Prescott, 1980). Once in the cell, it is metabolized by UGT and SULT to non-toxic metabolites and by CYP2E1 and CYP1A2, the two enzymes involved in APAP metabolism in mice (Snawder et al., 1994;Lee et al., 1996;Genter et al., 1998) to the toxic metabolite NAPQI. When the level of GSH, an intracellular antioxidant, is sufficiently high, NAPQI is detoxified by GST. If GSH levels get depleted through the detoxification reaction, NAPQI forms adducts in the mitochondrial proteins leading to the formation of ROS (reactive oxygen species), which get amplified through activation and translocation of JNK to the outer mitochondrial membrane. This also leads to further depletion of GSH. High levels of ROS trigger the MPT (mitochondrial permeability transition), leading to the impairment of ATP production. In our model, a cell is considered dead when its ATP level, initially around 5000 µM, falls below 100 µM. This value was chosen arbitrarily low relatively to the initial value (2% of the initial value) but selected simulations varying this parameter indicated that the parameter has only a negligible effect.
Given a cellular input exposure (concentration of APAP outside the cell over time) in the culture medium, the fraction of dead cells over time is obtained by simulating an ensemble of 3251 cells, each cell differing from the other one reflecting the cell-to-cell variability in the population (Furusawa et al., 2005;Sigal et al., 2006;Spencer et al., 2009). The population size of 3251 has been chosen to match it with the number of cells considered in the spatial temporal model of a virtual lobule.
Cell-to-cell variability is modeled by varying all parameters that depend on enzyme levels (i.e maximal velocities and production rates). The distribution of the varying parameters in the cell population is assumed to be log-normal around the mean parameter value, with a certain coefficient of variation CV (standard deviation over the mean). The mean values and the CVchosen to be the same for all values -are fit parameters. The maximal velocities of CYP2E1 and CYP1A2 is varied in the cell population according to the experimentally measured distribution (see SI section 1.10).
It was verified that the results do not depend critically on the number of cells chosen. Using 160 cells instead of 3251 cells gave for the same seed an only very slightly different parameter distribution such that the fraction of death cells differs by only ~ ±5% from those where 3251 cells have been used.
The equations are presented in the supplementary material. They were implemented in MATLAB. Physiologically relevant parameter ranges were determined from literature for each parameter and parameters having an influence on the simulation were estimated using the in vitro fraction of dead cells after a 24h exposure to 0.25mM, 0.5mM, 1mM, 2mM, 4mM or 8mM APAP.

In vivo-PD compartment model of APAP
The in vitro PD model was modified in the following way to represent the in vivo settings: i) The input APAP concentration outside the cells was replaced by the concentration of APAP in blood through two main strategies: (1) the pharmacokinetic model of APAP is plugged as an input concentration profile for 3215 cells representing a piece of a liver lobule (used in CL-model 1: Fig. 3(A.1, A2) and CL-model 2: Fig. 3(B1, B2) and (2) as obtained from a multi-compartment model (CL-model 3, Fig. 3(C.1)). The latter model considered a peritoneum compartment into which APAP was injected. From there, APAP is then released to a blood compartment from which the kidneys and liver eliminated it. The kidney effect was modeled by a first-order elimination term in the blood compartment. The liver uptake is modeled as the sum of the uptake by 3215 cells within a well-mixed liver compartment ( Fig.  3(C1)), and later upscaled. The 3125 cells were chosen to match the cell population size of the spatial-temporal (ST-) model ( Fig. 3(D.1)) below. For upscaling, this sum is multiplied with the estimated number of equivalent liver lobule pieces, i.e. the ratio of the total number of cells in the liver and the number of cells in the spatial-temporal lobule (see supplementary material for more details, section). ii) The CYP enzyme activities ( ) were modified to capture the higher enzyme activity in vivo compared to the in vitro case (based on experimental data, see section 2.5 and 2.7 . iii) The initial GSH concentration was modified to two times lower initial concentrations compared to the in vitro case (based on experimental data, see SFig. 4). In vivo-PD spatialtemporal model of APAP The cells uptake in the space-free in-vivo-PD model of APAP (sect. 2.10 above) was replaced by a spatial lobule of 3215 cells, representing each individual hepatocyte within liver lobule tissue microarchitecture. The in vivo-PD model of APAP (from 2.10 above) was now solved in each hepatocyte. The spatial lobule represents a full lobule in (x, y)-direction with a height of 10 hepatocytes in z-direction. APAP enters the sinusoidal (liver capillary) network of the model liver lobule with the blood by the portal veins and hepatic artery, that are not distinguished in the model for the blood flow. The hepatocytes were assumed to take up APAP by the same mechanism as for the well-mixed liver model (sect. 1.11.3), whereby the blood concentration of APAP in the sinusoid at the position of each individual hepatocyte was used. The blood left the liver by the central vein hence representing the lobule outlet. The liver outlet is then connected back to the blood compartment. The entire liver compartment was modeled as parallel arrangement of identical liver lobules such that each sub-volume of blood entering the liver would pass only one liver lobule. Blood flow, transport in the blood vasculature and intra-cellular reactions were computed in only one lobule that was considered as "representative" -whose contribution was multiplied by the total number of liver lobules. Because the liver is now considered as a spatial compartment, the coupling to the blood compartment involves the liver flow rate as an additional parameter. This parameter was calibrated manually around experimental data (see SI section 1.6.4).
For the detailed micro-architecture of the lobule in the computer simulation, statistically representative liver lobule architectures obtained from three-dimensional volume reconstructions of confocal laser scanning micrographs where chosen (Hoehme et al., 2010;Hammad et al., 2014).
Simulations directly in the 3D volume reconstructions were not feasible as the 3D volume data sets did not represent the portal veins, while determination of the flow boundary conditions require an entire liver lobule in (x,y)-plane. The construction of representative liver lobules described in (Hoehme et al., 2010;Drasdo et al., 2014) circumvents this problem. It displays portal triads, the sinusoidal network, the hepatocytes, and the central vein.
Within the sinusoidal network, steady-state blood flow is assumed, and calculated in each sinusoid from the total entering flow via the portal vein and hepatic artery. Poiseuille-flow is considered, which relates the volume flow rate to the pressure difference along each sinusoid of the sinusoidal network. To account for the dependence of the effective blood viscosity with the sinusoid diameter (Fahraeus-Lindqvist effect) an empirical effective blood viscosity model was chosen (Boissier et al., 2021).
The transport of APAP with the blood flow, and its uptake by hepatocytes, is simulated by partial differential equations mimicking the advective transport of molecules within the blood, and the transport of APAP from the blood into and out of each hepatocyte. The in vivo-PD model is executed in each individual hepatocyte of the liver lobule with the hepatocytes' kinetic parameters varying due to cell-to-cell variability, now in space. (The equations and further details are provided in the supplementary material.)

Fitting procedure
A similar fitting procedure has been applied to the in vitro and in vivo data.
As the volume of the solution in vitro was large such the APAP concentration only changed negligibly over the measurement period of 24 hrs. Hence, simulation of the cytochrome P450negative cells was unnecessary, and the fitting procedure could be focused on the CYP-positive cells only, which contained the NAPQI, SULT and the UGT conversion reactions.
To identify the values of the model parameters to the in vitro data, the following strategy was performed. First, the initial values of the 34 parameters and the ranges over which they can vary were chosen according to literature values. In the next step a simple parameter sensitivity on these 34 parameters was performed to identify those parameters, that influence the hepatotoxicity. Each parameter was separately varied by several orders of magnitude (divided and multiplied by 100) for one fixed concentration of APAP, while the other parameters were kept constant.From this analysis, 17 parameters were determined to only have a negligible impact on the hepatotoxicity, quantified by the simulated fraction of dead hepatocytes, and were thus fixed to their initial values. To calibrate the remaining 17 sensitive parameters, the log-likelihood function using all the available data of in vitro toxicity was maximized using the CMA-ES algorithm (a global convergence algorithm for non-linear functions). The finally computed dead cell fraction used the entire cell population size i.e., CYPpositive-and negative cells into account in the denominator. (More details and equations can be found in supplementary information.) For the in vivo model fit, no sensitivity analysis to separate sensitive and non-sensitive parameters was performed as the number of parameters (5) was small, and the ranges were better known. However, as the APAP concentration varies markedly in the blood, now both, the cytochrome P450 positive and negative cells were taken into account, whereby the enzyme activities of the UGT and SULT-pathway were assumed to be the same in all hepatocytes.

Establishment of in vitro and in vivo data of APAP hepatotoxicity
In a first step we determined in vitro and in vivo APAP hepatotoxicity (corresponding to Fig. 1(I), (V)).
In vitro: Primary hepatocytes were isolated from male C57BL6/N mice and cultivated in a 3D configuration between two soft gel collagen layers (SFig. 1). Six different concentrations ranging between 0.25 mM and 32 mM APAP were tested and cytotoxicity was analyzed by propidium iodide (PI) staining and cell morphology after a 24 h exposure period ( Fig. 2A-C, SFig. 1C, D). For the 4mM dose, cytotoxicity was analyzed at 5 different time points (Fig. 2C). The fraction of dead cells increases to almost its saturation value after 10 hours. Analysis of hepatocytes of three different mice resulted in a concentration-dependent increase of cytotoxicity, where approximately 40% of all hepatocytes lost their ability to exclude PI at concentrations of 4 mM and 8 mM APAP (Fig 2C). At even higher, concentrations (16 and 32 mM, not shown in Fig. 2C) a sharp increase of the concentration response curve was obtained, leading to cell death of all (i.e. not only CYP-positive) hepatocytes.

In vivo:
To determine the APAP-induced hepatotoxicity in vivo (Fig. 2D-F; SFig. 2), the same mouse strain (C57BL6/N, male) was studied that has been used to generate the in vitro data with cultivated hepatocytes. APAP was intraperitoneally administered at five doses ranging between 89 and 500 mg/kg body weight and livers were analyzed 24 h after administration (SFig 2A). The result demonstrates a dose-dependent increase of the pericentral dead cell area (SFig 2B). For doses up to 375 mg/kg the dead cell area corresponds approximately to the pericentral fraction of hepatocytes that express CYP2E1 (SFig 2C). Next, the pericentral dead cell area was quantified in relation to the total area of the liver tissue to calculate the fraction of dead cells in vivo. Mean values and standard deviations of three mice per dose showed a dose-dependent increase in the fraction of dead hepatocytes, with 89 mg APAP/kg b.w. representing the NOAEL and 158 mg APAP/kg b.w. the LOAEL (Fig 2F). These in vivo-hepatotoxicity data serve as reference for comparison to all in vitroin vivo extrapolation strategies.
In the next step, the classical extrapolation strategy based on AUC and Cmax has been studied.

Direct PBPK based in vitro to in vivo extrapolation: equivalent exposure does not mean equivalent toxicity
For the classical extrapolation scheme (Fig. 1(I)-(V)), based on AUC and Cmax (Figs. 2G, H), first the pharmacokinetics of APAP in the blood plasma was investigated by a combined experimentalmodeling strategy ( Fig. 1(IV), Fig. 2I).
For determination of the dose-dependent drug concentration at the site of action, four doses (56, 167, 300 and 450 mg APAP/kg b. w.) were administered to male C57BL6/N mice and blood was taken from the portal vein to determine APAP plasma concentrations. The data was collected 30 minutes after injection for all doses, and over additional time points to determine the temporal kinetics for the 300 mg/kg dose. A simple pharmacokinetic (PK) model (Methods, SI) was then fitted to the data and used to simulate concentration time curves for all APAP doses that have been studied for hepatotoxicity in vivo (Fig 2I). The concentration (C(t))-curves show a sharp increase after administration of APAP up to a maximum value Cmax that increases with the dose, and a first-order decay after that maximum. The integral under the simulated C(t)-curves corresponds to the AUC. Hence the C(t)-curves permit to attribute a unique Cmax and an AUC value to each dose.
For the AUC approximation, a concentration-equivalent is defined by performing a time integration over the C(t) curve, generally by = ∫ ( ) 0 with t=0 denoting the injection (administration) time point of the drug, T the time at which the toxicity is measured, x {vitro, vivo} identifies whether the in vitro or in vivo values are used. AUCx/T again has the unit of a concentration. Hence by equivalently equalizing either AUCvitro = AUCvivo or the concentration values defined by AUCvitro/T= AUCvivo /T, this procedure allows plotting the fraction of dead cells in vitro and in vivo for equivalent exposures (Fig 2G). AUCvivo /T denotes the average APAP blood concentration to which the hepatocytes are exposed to during the observation period T in vivo, so that the aforementioned approximation assumes that this average determines the hepatotoxicity. The administered in vitro APAP concentration corresponds to Cvitro(t=0), which is generally a good estimator for Cvitro(t) in the time interval [0, T] for T=24h if the volume of the culture medium is much larger than that of the volume of the hepatocyte population i.e., Cvitro(tT)  Cvitro(0).
For equivalent exposures based on AUC the fraction of dead (PI positive) hepatocytes shows a large deviation of the predicted from the measured value in vivo (Fig. 2G).
As an alternative to the AUC, the Cmax values of the PK-curve (C(t)) has been approximated with the drug concentrations administered in vitro (Cvitro  Cmax; Fig. 1(II, III)).
The intention for this approach is the assumption that the maximum APAP blood concentration Cmax is a good estimator for the APAP concentration used in the in vitro experiment, Cvitro(0), and that this maximum concentration determines the hepatotoxicity. Hence, again each administrated dose could be associated with a dead cell fraction through its corresponding drug concentration value, permitting to assign a hepatotoxicity value to each dose (Fig. 2H).
However, the comparison of dead cells based on Cmax did not result in a good agreement of in vitro and in vivo data (Fig 2H). One difference is the lower slope of the concentration effect curve in vitro. Moreover, the LOAEL occurs at a slightly lower concentration in vitro than in vivo, which contrasts the AUC-based comparison.
In conclusion, equivalent exposures in vivo and in vitro based on Cmax and AUC do not lead to equivalent toxicity. A second problem of this approach is the large difference of the results obtained by the AUC-and the Cmax-based extrapolation: a-priori, one cannot know if basing the extrapolation on Cmax or AUC leads to similar predictions, or if one or the other leads to a better prediction. Thus, the 'direct PBPK model extrapolation' fails to predict the correct in vivo toxicity. This calls for a refinement of the models used for extrapolation.

Setting up the intra-cellular detoxification (PD) model in vitro and in vivo
An alternative approach to the above-described direct PK-based extrapolation is a pharmacodynamic model (PD) mimicking the effect of APAP on each cell of a population (Fig 3A.1, Fig. 4A). The processes considered include diffusion of APAP between extracellular space and cytosol, and the key processes of APAP metabolism.
The final in vivo PD model consists of an intracellular PD-model module, and a model module specifying how the APAP feeds into the cell (Fig. 3A.1, Fig. 4A). The intracellular processes considered are the detoxification by UGT and SULT, and the activating metabolism by CYP2E1 and CYP1A2 (Fig. 4A). CYP enzymes lead to production of NAPQI which later binds to proteins. This triggers production of reactive oxygen species (ROS) modulated by JNK. After ROS accumulate to high levels, the cell reaches membrane permeability transition which impacts ATP production eventually leading to cell death (Snawder et al., 1994;Lee et al., 1996;Genter et al., 1998). The processes, which don not have any impact on the ATP levels (downstream processes) are not modeled here in order to reduce the complexity of the model.
Each cell in the model contains the same ODE system. Briefly, for all parameters of the ODE model a range of physiologically relevant values was used based on published data. To account for the cellto-cell variability the distribution of the enzyme activity parameters in the cell population was assumed to be log-normal with a coefficient of variation (CV, ratio of standard deviation to the mean), which is part of the estimated parameters (cf. section 2.10.2).
Not all cells are CYP-positive as confirmed by immunostaining of CYP enzymes which displays zonation around central veins (SFig. 5). After an image processing step, 50.9% of cells were considered CYP2E1 positive and 60.5% of cells are CYP1A2 positive. Cells that are CYP2E1positive are also CYP1A2-positive, but not necessarily vice-versa. A gradient of activity is assumed from highest CYP activity close to the central vein to lowest CYP activity in the distant region from the central vein. The gradient values were estimated based on the analysis of images stained for the cytochrome P450 enzymes assuming that the gradient of the gray values translate into a CYP-activity gradient. Because the cells considered in the in vitro experiments are harvested from the same animals used for the in vivo toxicity analysis, we assume the same fraction of cells expressing the gradient. The CYP enzymes' Vmax -values were thus not subjected to the cell-to-cell variability according to CV.
Because the medium volume is much larger than the total volume of all cells together (which corresponds to the in vitro reaction volume), the concentration in the medium does not vary significantly. The APAP concentration Cvitro(t) in the medium corresponds to a slightly linearly decreasing function that was estimated from experiments as an input function of each cell APAP metabolic model.
A simulation of the in vitro setting thus consists in solving the ODEs involved in 3215 cells (number of cells involved in a piece of virtual lobule as described later) with a given initial concentration of APAP in the medium over time. The cells which have an ATP level below 100 µM are considered as dead. Fig. 1(VI). For the explanation, see text. Green text indicates those parameters that were fitted in the respective modeling step. The prime in B.2 and C.1 indicate parameters used for the in vivo simulation.

Figure 3: Sketch of model-based extrapolation strategies, detailing
The parameters of the PD model were calibrated based on data of the in vitro experiments in primary mouse hepatocytes. An excellent agreement between experimental data and PD model could be obtained (Fig. 4B). However, due to the insufficient number of data points and the relatively wide range of physiologically compatible parameter values, no unique parameter set could be identified. Several fits of equivalent agreement quantified by their very similar standard error, were performed using different random seeds, which led to largely different parameter sets. This is expected to impact on the extrapolation to the in vivo toxicity.

Figure 4: Non-coupled PK/PD model strategy. (A) Scheme of the APAP metabolic model where the PK model is plugged as an input, (B) best fit on in vitro toxicity data over 40 different random seeds (shown is a particular seed, label 13), extrapolation to in vivo fraction of dead cells after correcting CYP enzyme activity and GSH concentration with (C) same permeability in vitro and in vivo (D) and with lower permeability in vivo.
In the next step the parameter sets obtained from the fit of the intracellular model to in vitro experimental toxicity data were used to extrapolate to in vivo toxicity, considering potential in vitro / in vivo differences of the intracellular parameters and in the exposure of hepatocytes to APAP (Fig.  3A).
In a first cycle of model extrapolations, we kept the same parameters for the Cytochrome P450 enzymes and GSH in vivo as in vitro but extrapolations totally failed (Cellière, 2016). As from those model simulations CYP2E1 and CYP1A2 activities (the two major APAP activating enzymes in mice) have been found to have the strongest impact on the model-based in vitroin vivo extrapolation, they were analyzed in cultivated hepatocytes and freshly isolated hepatocytes, whose activity is known to be similar to those in liver tissue (Godoy et al., 2016) (Godoy et al., 2016) (SFigs.  3, 4). These analyses permitted to quantify the decrease of both enzymes in culture compared to the in vivo situation.
Compared to hepatocytes used for the in vitro toxicity experiments, the in vivo hepatocytes' CYP activities was experimentally found to be lower by a factor of 3.3 for CYP2E1 and 1.8 for CYP1A2 in vitro than in vivo (SFig. 3), while the GSH concentration in vivo was determined to be 2 times lower in vivo than in vitro (see SFig. 4).
Concerning the difference in APAP exposure in vivo and in vitro, our first strategy consisted in computing the concentration of APAP in the liver from the previously introduced PK model assuming that the number of APAP-molecules taken up from the blood by the cells during one passage of blood through the liver is negligible compared to the number of APAP-molecules in the blood, and that each hepatocyte in vivo is exposed to approximately the same APAP blood concentration. The approach is simple and computationally cheap because one does not need to represent explicitly the uptake of cells from the blood and because the PK-model and the intracellular PD-model can be run independently of each other. The procedure in this case is: (i) fit the intracellular parameters for the PD model with in vitro toxicity data (Fig. 3A.1), (ii) adjust the enzyme activities and metabolite concentrations, (iii) fit the extracellular parameters of the PK model on the APAP blood concentration data (Fig. 3A.2). Once the parameters are calibrated, the in vivo PD model is simulated for each toxicity dose and compared to data (Fig. 3A3).
This approach significantly overestimated the fraction of dead cells (Fig. 4C) but yielded good extrapolation results if the permeability for APAP in vivo was assumed to be significantly lower than in vitro (Fig. 4D). This may be justified as in vivo APAP in the blood has to pass the space of Disse which contains extracellular matrix, such as collagens. However, the assumption that the uptake of APAP molecules is negligible compared to the number of APAP molecules remaining in the blood turned out to be violated for all parameters. When considering a lower permeability in vivo, the hypothesis is verified but the total hepatocytes uptake is way too low and thus not realistic either (less than 1%) (Prescott and Wright, 1973;Dai et al., 2006;Malfatti et al., 2020).
For this reason, we dropped this hypothesis and moved to another approach that explicitly models the interaction between the evolution of the blood concentration and the uptake of the cells. Coupled PK/PD model scheme with assumed storage from ECM, (H) fit of the in vitro-PD model on in vitro toxicity data for random seed 25 that lead to lower permeability, (H) fit of the PK/PD in vivo model extracellular parameters on APAP blood in vivo concentration data using random seed 25 and assumed ECM storage, (I) extrapolation to in vivo toxicity, (J), fit of the PK/PD in vivo model with additional ECM storage hypothesis, (K) and with extracellular parameters and effective permeability in vivo on APAP blood in vivo concentration data for random seed 25, (K) extrapolation to in vivo toxicity (L).

Stepwise, independent fit in vitro, then fit in vivo blood concentration extracellular parameters with a coupled PK/PD model approach fails (CL-model 2)
The second model-based strategy integrated the PK and PD in one single PK/PD model that simultaneously mimics the APAP-PK in the blood, the elimination of APAP from the blood compartment, and the intracellular detoxification pathway reactions ( Fig. 5A; Fig. 3B.2). For this purpose, the transport of injected APAP into the blood from the peritoneal cavity as well as its elimination by the liver itself and of all other possible clearance pathways, as for example removal by the kidney, was computed. The elimination by non-liver sources was lumped into one first-order elimination term. The removal of APAP by the liver was modeled by summing the uptake of APAP by all cells of the liver. In practice one sums over the number of cells in a representative slice of the liver lobule of about 10 hepatocytes in height. Then for each cell of the liver lobule slice the set of ODEs representing the intracellular PD model (Fig. 4A) is solved for each time step. This model thus corresponds to a multi-scale model with a well-mixed blood compartment coupled to individual space-free cells (Fig. 5A).
The parameters of the PK/PD are determined in a multi-step procedure ( Fig. 3B.1, B.2).
(1) Firstly, that set for the parameters of the intracellular model is chosen as a starting parameter set that generated the best fit of the APAP in vitro-PD model on in vitro toxicity data (Fig. 5B, Fig. 3B.1). (2) Secondly, the activities of CYP2E1 and CYP1A2, as the one of GSH were replaced according to the factors determined in section 2.3. (3) Thirdly, the so modified intracellular model was integrated with the blood and peritoneal cavity compartments to fit the extracellular parameters such that the PK data could be captured (Fig. 5C, Fig. 3B.2). An extra constraint was added to represent elimination by other factors than liver to represent 15% of the total bioavailable dose (Prescott and Wright, 1973;Dai et al., 2006;Malfatti et al., 2020) Using this second strategy, the fit to the APAP blood concentration (Fig. 5C) misses out the lowest concentration data and still deviates markedly from the profile over time for the 300 mg/kg dose. Moreover, the APAP concentration decreased too quickly with time compared to the data in Fig. 5C. Finally, and most importantly, the fraction of dead cells is largely overpredicted for small APAP doses (Fig. 5D). We hypothesized that the permeability in vivo that is determined by both the hepatocyte membrane and the crossing of the space of Disse, may be smaller than in vitro (where only the hepatocyte membrane needs to be traversed), which would lower the APAP amount passing into the hepatocyte per time unit, so may reduce cell death in particular for small doses (McPhail et al., 1993). Permeability in vivo can thus be considered as an effective permeability.
To test this hypothesis now the extra-cellular parameters including the effective permeability were fitted. While this improved the fit quality it was still not enough to fully capture the blood concentration data (Fig. 5E). In particular, the profile over time is closer to the data but the concentration for the lowest concentration data is missed out. Interestingly, the permeability re-fitted led to a lower value of the predicted in vivo toxicity, despite the fraction of dead cells in vivo was still overestimated (Fig. 5F). (Notice, that refitting of the permeability to the in vivo toxicity data after fitting of all other parameters as in Fig. 5C would modify the fit to the PK-data, hence provides no option.) Since the simulated blood concentration overall fit to the PK data was not great, neither in Fig. 5C, nor in Fig. 5E, a storage mechanism was introduced representing the hypothesis that APAP is transported with a delay into the hepatocytes as it would be the case if it is first stored in the space of Disse (e.g. by adhering to the ECM) and then released from the ECM for transport into the hepatocyte with a delay (Fig. 5G). As lowering the permeability improved the fit quality at the previous step, we chose here the set of intracellular parameters fitted on in vitro data for which the permeability was the lowest to test whether this could eliminate the overprediction of in vivo toxicity at small doses (Fig. 5H). The storage mechanism improves significantly the capturing of the blood profile over time (Fig. 5I). The lowest concentration data is also captured accurately (Fig. 5I). However, even these two steps, (1) adding a delay in the transport of APAP from the blood into the hepatocyte, (2) choosing the parameter set with the smallest permeability, did not lead to a better toxicity prediction (Fig. 5J).
In a last attempt, letting the permeability in vivo be a fit parameter (as in the transition from Fig. 5C to Fig. 5I), slightly improves the fit to the in vivo APAP blood concentration (Fig. 5K) but the model could still not predict the in vivo fraction of dead cells (Fig. 5L). This is expected as we intentionally chose a random seed parameter set for which the in vitro permeability was already low.
In conclusion, fitting the in vitro data in a first step (Fig. 3B.1), then replacing in vitro cytochrome P450-and GSH activities in that model by their in vivo counterparts, and finally fitting a model integrating the so adapted intracellular model with a blood pharmacokinetics (PK) model to the PK data in a last step (Fig. 3B.2), failed to yield reasonable predictions of the in vivo toxicity (Fig. 3B.3).
Hence, we asked the question if the integrated PK/PD model if it would be fitted to both the in vitro and PK data in a single step, could yield a valid in vivo prediction. In this last fitting strategy (Fig. 3C.1, C.2; Fig. 1(VI-A, B)), the parameters of the intracellular APAP detoxification model were fitted to the in vitro toxicity data, and in the same step i.e., simultaneously, the extracellular parameters of the full model, composed of the intracellular in vivo model that is obtained by modifying the enzyme activities of CYP2E1, CYP1A2 and glutathione in the in vitro model (that is fitted in the same step) according to the experimentally measured factors and the in vivo compartments, are fitted to the in vivo blood (PK) concentration data ( Fig. 3C.1, C.2). The fitness function to minimize was the sum of both fitness function previously used for both data sets individually. For simplicity the storage mechanism of Fig. 5E was not included. The model is able to capture well the data in vitro (Fig. 6B) along with in vivo blood concentration data (Fig. 6C). When extrapolating to in vivo, the model was able to predict the in vivo toxicity data with sufficient accuracy (Fig. 6D; Fig. 3C.3). However, for a different random generator value during the fit procedure (for example, starting the fit procedure from different points in the parameter space) an overestimation of the damage at small doses was observed (SFig. 6D), even though the quality of fit to the in vitro and PK data did not change (SFig. 6B, C). Nevertheless, that extrapolation still performed better than in the other computational model-based extrapolation strategies.

Simultaneous parameter determination in vitro and in vivo blood
These results indicate that simultaneously fitting the intracellular parameters to the in vitro toxicity and extracellular parameters within a model integrating the intracellular model and the relevant body compartments to the PK data, is able to generate a very good agreement to that data on one hand, and a valid extrapolation to the in vivo toxicity data on the other hand. This means that the mechanisms of the model considered here are enough to explain the data available. However, redoing the fit also yielded parameters sets for which the in vivo hepatotoxicity prediction failed to be accurate enough although in vitro and PK data were well captured. This indicates that in order to take a final decision on whether the studied set of mechanisms is de facto explaining the hepatotoxicity further information is necessary. This can for example be through more narrow parameter ranges or additional in vitro toxicity measurements, either determining the hepatotoxicity in vitro for additional doses of APAP or measuring at more time points.
In the liver lobule, the chemical reaction volumes are given by the hepatocytes in a specific spatial organization due to liver lobule zonation. Downstream hepatocytes see lower APAP concentrations that upstream hepatocytes. In order to study in how far this spatial organization may modify the detoxification of APAP, finally the effect of microarchitecture on the APAP detoxification was studied within a spatial-temporal micro-architecture liver lobule model. In the final step, a multi-level hepatic lobule model was created that considers the complex lobular microarchitecture to include detoxification effects emerging from the spatial organization of hepatocytes within the liver lobule (cf. Fig.1(VI-A, C), Fig. 3D). Such an approach has the fundamental advantage, that it could directly represent disease-related architectural distortions as they occur in fibrosis, steatosis or cirrhosis (Ghallab et al., 2019b(Ghallab et al., , 2021. E.g., in periportal fibrosis, cytochrome P450-negative hepatocytes are partially replaced by scar tissue (Ghallab et al., 2019a), while in septal fibrosis, fibrotic streaks connect central veins of neighboring lobules hence reduce the number of CYP-positive hepatocytes (Ghallab et al., 2019b).

Considering lobular microarchitecture by a multi-level virtual liver lobule (ST-model)
As a proof of concept, the simulations have been performed in a virtual liver lobule generated by sampling from statistical distributions of geometric parameters defined to characterize 3D lobule architecture in confocal laser scanning micrographs of healthy mice (Hoehme et al., 2010) (cf. Fig.1 (VI-C))). The geometric parameter distributions were obtained by analysis of a few tens of liver lobules that could not be directly used for the simulations as none of them displayed an entire lobule. The virtual lobule permits to construct an entire lobule in (x, y)-direction displaying all portal triads ( Fig. 7A). This ensured that flow and pressure boundary conditions could be uniquely determined. The representative lobule model represents each individual hepatocyte within the lobule as well as the sinusoidal (the capillary) network. Blood flow is modeled by Poiseuille flow, APAP transport within the lobule by a partial differential equation, and APAP metabolism by the ODE-based PD model that is solved within each individual hepatocyte (SI). The same parameters for the ODE model were used as obtained in the previous section. Cytochrome P450-enzymes were zoned as experimentally observed (cf. SFig. 5).
This virtual spatial liver lobule replaces the set of independent, well-mixed hepatocytes of a classical (compartment) model (see definition above). The APAP injected intraperitoneally enters the blood compartment, which feeds the portal veins of the liver lobule and is at the same time fed via the central vein of the liver lobule (Fig. 7A). Hence, the concentration at the inlet of the liver lobule is set according to the concentration in the blood compartment, while the blood compartment receives inputs from the peritoneum and the central veins. In this approach, the liver blood flow rate had to be included in the model to couple the blood compartment to the liver lobule. This is an additional parameter, which was adjusted in the narrow range of values found in published references (see Supplementary Material, SFig. 7). All other parameters were chosen from the well-mixed model CLmodel 3. Once adjusted, the full multiscale ST-model reproduces with similar quality as the wellmixed model the blood concentration data and captures closely the fraction of dead cells in vivo data ( Fig. 7B, C).
In the spatial-temporal (ST)-model, APAP enters the periphery of each liver lobule via the blood of the portal vein and the hepatic artery. From there it passes through the sinusoidal network along the hepatocytes until it drains into the central vein. The periportal hepatocytes do not express CYP2E1 or CYP1A2 and therefore CYP450-mediated hepatotoxicity does not take place in this part of the liver lobule. Nevertheless, glucuronidation and sulfation reactions of APAP take place in the periportal hepatocytes, and the formed phase 2 metabolites can no longer contribute to pericentral toxicity. Therefore, periportal hepatocytes already eliminate a fraction of APAP before it reaches the pericentral region without generating cytotoxicity. This effect is automated accounted for in the STmodel but lacking in the compartmental PD model that lumped all hepatocytes with the same enzyme equipment into one well-stirred compartment (Figs. 1(VI B, C), 5A, 6A). As a consequence, the actual APAP concentration to which the CYP positive cells are exposed to, is smaller in the STmodel (Fig. 7) than in the well-mixed CL-model 3 (Fig. 6).
Moreover, the multilevel PD spatial-temporal (ST-)model gives the opportunity to visualize APAP concentrations and toxic consequences three-dimensionally in a virtual liver lobule. Stills of the 24 h after administration of 89, 281 and 500 mg APAP/kg b. w illustrate the dose-dependency of pericentral cell death and allow a comparison to the histological slides of the corresponding doses as those exemplarily presented in SFig. 5.
Because the extra flow rate parameter introduces an additional degree of freedom, even if chosen from a small parameter range, a direct comparison with the well-mixed liver approach on previous section is difficult. However, different from the well-mixed model, the parameters of the ST-model were not calibrated with the in vitro and the PK-data. Hence, the parameters have not been optimized for this model, so it is likely that a parameter set may exist for which the agreement to the data is better than that in Fig. 7. What is hindering this fit is the long simulation time with the ST-model, which makes a direct parameter identification impossible with the current generation of computers. It is possible to circumvent this difficulty by an iterative parameter identification strategy that would include iterations between data, well-mixed and ST-model. As the difference in the agreement between ST-model and well-mixed model where so small, this iterative procedure was not applied here and will be presented for another application elsewhere. Finally, we believe that this spatial representation could be of major interest in case diseases with architectural modifications such as fibrosis or cirrhosis.

DISCUSSION
An important question in toxicological research is how accurate in vitro toxicity data can be extrapolated to the toxicity in vivo. Guided by this question, we studied a number of extrapolation strategies for APAP-overdose induced liver injury, starting with the classical extrapolation strategy, followed by mathematical model-based extrapolation strategies varying the complexity of the model and the parameter identification strategy. APAP was chosen as its toxicity mechanism was considered as relatively well understood, so mathematical models could be built based on a scheme of toxicity mechanisms that was believed to largely reflect consensus, in particular involving NAPQI-induced hepatocyte death. In so far, the model results may also inform about in how far the consensus mechanisms permit to explain in vivo toxicity.
To permit evaluation of the accuracy of the extrapolation strategy, first hepatotoxicity was experimentally determined over a wide range of concentrations in vitro and doses in vivo, complemented by measurements of the drug pharmacokinetics (PK) in the blood, activity determinations of the enzymes CYP2E1, CYP1A2, and glutathione concentration, respectively. The in vitro measurements use cultivated primary hepatocytes of the same mouse strain that was used in vivo.
A number of different approaches (some in different variants) was then studied with regard to their ability to predict hepatotoxicity in vivo from in vitro hepatotoxicity data: 1. Str1 (Fig. 1(I-V)): The classical extrapolation scheme based on AUC and Cmax of the blood PK data. 2. Str2 (Fig. 3-A): A three-step model approach, in which first the parameters of a PD model were determined by fitting the PD model to in vitro toxicity data. Then, in a second step, those enzyme activities measured to change in vitro from their in vivo values were modified in the PD model and the so modified model was in a third step used to predict hepatotoxicity in vivo by receiving its input by a blood pharmaco-kinetic model for APAP. 3. Str3 (Fig. 3B): A three-step approach, modifying only step 3 of the previous approach (Str2.) by embedding the modified PD model into a multicompartment in vivo model such that the pharmacokinetics was now a model outcome. The intracellular parameters were first fitted on the in vitro toxicity data and then the extracellular parameters were fitted on the in vivo kinetic data. 4. Str4 (Fig. 1(VI-B), Fig. 3C): A two-step approach using the same model structure as in the previous approach (Str3) but now fitting simultaneously the in vitro and in vivo model parameters to the in vitro hepatotoxicity data and the PK data. The hepatotoxicity of the model with the so determined parameters corresponds to the in vivo hepatotoxicity prediction. 5. Str5 (Fig. 1(VI-C) Fig. 3D): A spatial-temporal model resolving liver microarchitecture, using the same parameters as determined in strategy (Str4.) plus a liver flow rate parameter The classical approach (Str1) based on AUC and Cmax of the blood PK-data, complemented by a population PK model, failed to yield an accurate hepatotoxicity prediction. The predicted hepatotoxicity using the AUC is largely deviating from the experimental hepatotoxicity data. Using Cmax the predictions are overall markedly better, but the hepatotoxicity is overestimated for small APAP doses and moderately underestimated for large APAP doses. This points to the difficulty that because the AUC and Cmax-based extrapolation strategies do not yield the same prediction, one would have to know a-priori (which is not the case), which of the two characteristics, AUC or Cmax, should be chosen.
Hence, we studied whether pharmacodynamics (PD) models of the classical well-mixed type (Str2-4) or spatio-temporal (ST; Str5) model resolving micro-architecture may help to overcome the problems of the classical extrapolation strategy (all schematized in Fig 3).
The created model of the PD in a single hepatocyte reflects the consensus mechanisms of APAP hepatotoxicity including conversion of APAP to NAPQI catalyzed by CYP2E1 and CYP1A2, and conversion of APAP by SULT, and UGT. The underlying structure of the intracellular (cell: hepatocyte) PD model was part of each modeling/fitting strategy Str2-Str5.
In strategies (Str2: Fig. 3A.1, Str3: Fig. 3B.1) in a first step the parameters of the intra-cellular PD model were directly fitted to the in vitro data taking into account the fraction of hepatocytes in which the cytochrome P450 enzymes CYP2E1 and CYP1A2 are expressed. The intracellular PD model simulates the mechanism of action of APAP, including its activating and inactivating metabolism, generation of reactive oxygen species, mitochondrial toxicity, ATP depletion and cell death.
In a next step, for all strategies, the CYP2E1 and CYP1A2 activities as well as the GSH initial concentrations are of the intracellular PD model were adapted to their in vivofound values to obtain an "in vivo intracellular PD-model", before coupling the latter to a PK module (Str2: Fig.4; Str3: Fig.  5; Str. 4: Fig. 6; Str5: Fig. 7). Discrepancies as between Cmax and AUC for the classical extrapolation scheme (Str1) cannot occur anymore. A second main advantage of basing the toxicity extrapolation on a model that uses the same intracellular model structure in vitro and in vivo is that its mechanistic nature permits to easily incorporate the differences between in vitro and in vivo. This is particularly helpful when equivalent exposure in vivo and in vitro does not lead to similar toxicity. This is the case when cells in vitro differ, e.g. in their metabolic activities from corresponding cells in a tissue in vivo.
The coupling of the intracellular PD model and the PK model has been done in three different ways for Str2, Str. 3, Str4/5 (Str5 is a spatially resolved version of Str4).
For Str2 (Fig. 4, Fig. 3A) the PK model was assumed to remain unaffected by the uptake of APAP by the hepatocytes, which would be a reasonable assumption if the removed APAP per passage of blood through the liver would be negligible compared to the APAP blood concentration. This assumption was violated for most parameter sets that were able to yield a reasonable fit to in vitro toxicity and no extrapolation to in vivo toxicity data can be achieved. When considering a lower effective permeability in vivo, the uptake to cells is too low (< 1%)and thus not realisticalthough excellent extrapolation to in vivo toxicity data can be achieved.
To avoid such a violation of self-consistency, in Str3 the PK was mimicked by a mathematical compartment model representing peritoneum, blood and hepatocyte compartment instead of considering the PK-curve as input to the in vivo intracellular PD model (Fig. 3B). I.e., the PK becomes an output of the model integrating physiological body compartments and the intracellular detoxification pathway. In this case, the toxicity at small APAP doses was systematically, partially even largely, overestimated (Fig. 5).
Hence in a further strategy (Str4, Fig. 3C) we fitted the parameters of the intracellular in vitro PD model (to the in vitro toxicity data) and the parameters of the compartment model simultaneously (to the APAP blood concentration data, Fig. 6), which for several parameter sets yielded subsequently good in vivo hepatotoxicity predictions, however, with varying accuracy. This accuracyand reproducibilitymay be improved by additional sets of experiments as outlined below.
In a final step we studied the influence of tissue micro-architecture on detoxification within a spatial temporal model (Str5, Fig. 3D), which readily takes liver zonation into account. When entering the periportal fraction of the liver lobule in that model, APAP first passes hepatocytes that do not express the enzymes CYP2E1 and CYP1A2, which downstream form hepatotoxic products, but detoxifying UDP-glycuronosyl-transferases and sulfotransferases, which form products that do not harm the cell. Therefore, a fraction of APAP was detoxified by sulfation or glucuronidation before it reached the pericentral fraction of cytochrome P450 expressing hepatocytes. The sequential removal of APAP at the path of the blood through the liver from the portal tract towards the central vein led to small removals of APAP by cells not expressing CYP-enzymes with only minimal reduction in hepatotoxicity (Fig. 7).
While the approach can be very useful in case of strong architectural distortions or loss of fenestration in chronic liver diseases, for acute DILI has only a minor influence compared to the well-mixed in vivo PD model (Fig. 6).
Recently, (Ghallab et al., 2022) have identified a new mechanism that contributes to APAP-induced liver injury in vivo. APAP overdose causes a breach of the blood-bile barrier and leakage of bile acids into the sinusoidal blood, which are transported again to the pericentral hepatocytes via the sinusoidal uptake carriers. This results in accumulation of bile acids in the pericentral hepatocytes above toxic thresholds, and thereby aggravates APAP hepatotoxicity. Interestingly, blocking the sinusoidal bile acid uptake carriers strongly ameliorates APAP-induced hepatocyte death (Ghallab et al., 2022). This discovery was made possible using functional intravital imaging techniques which allows imaging at subcellular resolution (Hassan, 2016;Schneider et al., 2021a;Vartak et al., 2021). This novel mechanism was not considered in our modeling approach but can prospectively be integrated into the spatial temporal model. However, during the establishment of the modeling strategies shown in this paper, we hadbesides the mechanisms detailed in this work -introduced further black-box mechanisms that act dose-dependent to search for additional potential mechanisms, but even with those a better agreement of data and model than depicted in Figs. 6 and 7 and a higher reproducibility could not be obtained. However, integrating the bile canaliculi network in the STmodel as a direct approach may provide further insight.
In summary, the explanation for the data obtained was the better, the more the structure of the model resembled the real liver micro-architecture and physiology. The final model and fitting strategy permit a reasonable prediction of in vivo hepatotoxicity from in vitro toxicity measurements. However, the difficulty to reproduce the parameter sets upon multiple fit repetitions indicate that (1) the parameter ranges known for the enzyme activities of all reactions may still be too large, or / and (2) the number of data points may be still too small to ensure that the fitting algorithm always finds the same parameter set for those parameters that are identifiable. This would require a large set of extra experiments, which may perhaps be better resolved within a community effort. To finally decide on the model, accuracy is fundamental, otherwise it is not possible to take a final decision on the contribution of the CYP-catalyzed mechanism via NAPQI to hepatotoxicity for APAP overdose. Once this accuracy in the parameters is reached, a sufficiently accurate prediction of NOAEL and LOAEL following our modeling strategy Str4 and/or Str5 seems within reach.  The toxicity at a certain time T (often: 24h, or 3d) is considered as the time at which the toxicity is compared (II). A common consensus for extrapolation to in vivo is that the in vitro concentration value Cvitro is identified by the Cmax-value (III) of the corresponding pharmaco-kinetic (PK) curve for the drug in the blood (IV). The PK curve depends on the administered drug dose, so that by identifying Cmax and Cvitro, the dose value associated with that Cmax-value can be associated with the toxicity value associated with the corresponding Cvitro value (V). (Alternatively, to the Cvitro / Cmax values, the area under the curve (AUC) is used). The modelbased strategy mimics drug toxicity in a model by simulating the process of drug detoxification in vitro by the cells, whereby the toxicity pathway is explicitly represented in each cell and eventually integrated with a PK model (IV) and a compartment model of the organs of interest (VI). The latter may represent organ microarchitecture, here the liver lobule (VI-C) as repetitive minimal tissue unit or consider a well-mixed compartment (VI-B), both integrating an intracellular PD model (VI-A).
Simulations with (VI-B, C) directly yield the in vivo toxicity (VII).