Bottom-Up Physiologically-Based Biokinetic Modelling as an Alternative to Animal Testing

There is a growing need to develop alternatives to animal testing to derive biokinetic data for evaluating both efficacy and safety of chemicals. One such alternative is bottom-up physiologically-based biokinetic (PBK) modeling, which requires only in vitro data. The primary objective of this study was to develop and validate bottom-up PBK models of 3 HMG-CoA reductase inhibitors: rosuvastatin, fluvastatin, and pitavastatin. Bottom-up PBK models were built using the Simcyp® simulator by incorporating in vitro transporter and metabolism data (Vmax, Jmax, Km, CLint) obtained from the literature and proteomics-based scaling factors to account for differences in transporter expression between in vitro systems and in vivo organs. Simulations were performed for single intravenous, single oral, and multiple oral doses of these chemicals. The results showed that our bottom-up models predicted systemic exposure (AUC0h-t), maximum plasma concentration (Cmax), plasma clearance, and time to reach Cmax (Tmax) within two-fold of the observed data, with the exception of parameters associated with multiple oral pitavastatin dosing and single oral fluvastatin dosing. Additional middle-out simulations were performed using animal distribution data to inform tissue-to-plasma equilibrium distribution ratios for rosuvastatin and pitavastatin. This improved the predicted plasma-concentration time profiles but did not significantly alter the predicted biokinetic parameters. Our study demonstrates that quantitative proteomics-based mechanistic in vitro-to-in vivo extrapolation (IVIVE) can account for downregulation of transporters in culture and predict whole organ clearance without empirical scaling. Hence, bottom-up PBK modeling incorporating mechanistic IVIVE could be a viable alternative to animal testing in predicting human biokinetics. This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 International license (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution and reproduction in any medium, provided the original work is appropriately cited.


Introduction
The replacement, refinement, and reduction of animal use in research (3R) was first established in 1959 (Russell and Burch, 1959), with growing efforts in recent times to establish alternatives to animal testing in the risk assessment of xenobiotics (Chapman et al., 2013;Paini et al., 2017).In 2006, the REACH Regulation implemented in the EU demanded that testing on animals be done only as a last resort (EU, 2006).In 2009, an international association of validation bodies and agencies representing the USA, EU, Japan, and Canada adopted a memorandum to cooperate on reducing animal toxicity testing (ICATM; EURL ECVAM, 2009).Finally, under the EU Cosmetics Regulation (EU, 2009;Creton et al., 2009), a complete ban on testing of cosmetic products and ingredients in animals in the EU was imposed in 2013.Despite intensive efforts, the current strate-2 Methods

General PBK model structure and assumptions
All PBK models were developed using a population-based absorption, distribution, metabolism, and excretion (ADME) Simcyp ® simulator (Simcyp Ltd, version 17 release 1, Sheffield, UK).The design qualification and the quality assurance framework of Simcyp has been previously described (Jamei et al., 2013).An in-depth description of the sub-models in Simcyp has been previously published (Neuhoff et al., 2013;Jamei et al., 2014).
Briefly, the advanced dissolution, absorption, and metabolism (ADAM) model (Jamei et al., 2009) was used to describe the intestinal absorption process.Distribution was described using the full PBK distribution model to account for permeability-limited distribution in several organs (Jamei et al., 2014).We assumed that perfusion-limited distribution describes the distribution of the statin to all organs in the body except the intestine, liver, and kidney; and that there is instantaneous distribution of the statin between the intracellular and extracellular space.Where suitable in vitro data indicated that a statin was the substrate of one or more transporters, the relevant organ was converted to a permeability-limited distribution model.For the intestine, the permeability-limited distribution model was used for pitavastatin only.For the liver, the permeability-limited liver (PerL) model (Jamei et al., 2014) was used for all 3 statins.For the kidney, a permeabilitylimited mechanistic kidney model (Mech KiM) (Steffansen et al., 2013) was used for rosuvastatin and pitavastatin.This permitted us to incorporate transporter kinetics when modeling the distribution and disposition of the statins.The volume of distribution (V ss ) was estimated using tissue-to-plasma equilibrium distribution ratios (K p ) calculated using the method of Rodgers and Rowland (2007).We also assumed that for statins that undergo transporter-mediated biliary excretion, enterohepatic recirculation (EHC) occurs and that 100% of the biliary excreted drug is available for reabsorption in the intestine.Renal clearance was estimated from kidney transporter kinetics using the Mech KiM model.

Model development
Bottom-up models were constructed using in vitro experimental data obtained from the literature or predicted using quantitative structure-activity relationship (QSAR) models.Where more than one source of in vitro data was available for a certain parameter, an iterative simulate-and-refine approach was used to select the optimal value that provides the biokinetic parameter within the two-fold criterion described subsequently.For transporter and metabolism kinetics, two forms of in vitro data were available, i.e., intrinsic clearance (CL int ) or maximal transport/metabolic rate (J max /V max ) and Michaelis-Menten constant (K m ) val-in vivo data.The role of PBK modeling as an integrated testing strategy to assess xenobiotic exposure is increasingly recognized (Hartung et al., 2013;Tan et al., 2018).The major challenges for the bottom-up approach are the considerable uncertainties due to suboptimal predictions via IVIVE, especially for transporter kinetics (Margolskee et al., 2017;Shebley et al., 2018).Earlier studies reported that IVIVE using data obtained from human hepatocyte (HHEP) assays underestimated hepatic intrinsic clearance (Kotani et al., 2011;Menochet et al., 2012).It was proposed that this might be due to a reduced expression of hepatic transporters in isolated hepatocytes, as compared to that in the liver (Vildhede et al., 2015).We realized that quantitative proteomics could be used to account for the differential expression of metabolic enzymes and transporters between in vitro systems and in vivo organs and could be applied as mechanistic scaling factors to improve IVIVE in bottom-up PBK modeling.Hence, the primary objective of this study was to develop a bottom-up PBK model using mechanistic scaling factors based on quantitative proteomics measurements of transporter and metabolic enzyme expression levels for IVIVE to predict the in vivo exposure of xenobiotics in humans.This effort is part of an Organisation for Economic Co-operation and Development (OECD) case study to provide guidance on PBK modeling as an alternative to animal testing in the regulatory assessment of chemicals.
We chose three 3-hydroxy-3-methylglutaryl-CoA (HMG-CoA) reductase inhibitors: rosuvastatin, fluvastatin, and pitavastatin as substantial in vitro and in vivo data is available on them for our bottom-up PBK modeling.Importantly, the statins are differentially eliminated: rosuvastatin's elimination is transporter-dependent, fluvastatin's elimination is metabolism-dependent, while pitavastatin's elimination is dependent on both processes.Using only in vitro data from the literature, bottom-up models were built to simulate intravenous (IV) single, and repeated oral administration of each of the statins.While PBK models of these three statins have been previously published (Bosgra et al., 2014;Jamei et al., 2014;Zhang, 2015;Vildhede et al., 2016;Duan et al., 2017;Wang et al., 2017;Mitra et al., 2018), these were developed based on empirical scaling and the middle-out approach.By leveraging recent developments in transporter kinetics and quantitative proteomics that were not available at the time of development of the published models, we demonstrate the utility of mechanistic scaling factors for IVIVE via bottom-up PBK modeling.Additionally, we recognize that animal biokinetic data may be available or routinely generated in some settings and that it is advantageous to incorporate all available data including in vivo animal data to improve model accuracy.To demonstrate this middle-out approach option, we have used data from animal distribution studies to improve the prediction of the distribution profiles of rosuvastatin and pitavastatin.
1 doi:10.14573/altex.1812051sues.Wherever possible, J max /V max and K m were used as they allow modeling of saturation of enzyme/transporter kinetics.CL int was used when appropriate J max /V max and K m values were unavailable.Additionally, middle-out simulations were performed with data from distribution studies conducted with Sprague Dawley rats to optimize the prediction of K p for rosuvastatin and pitavastatin (Kimata et al., 1998;Nezasa et al., 2002).

Rosuvastatin model parameters
Rosuvastatin's physicochemical characteristics, ADME parameters, and source are provided in Table S1 1 .Absorption of rosuvastatin is suggested to be mediated by various uptake and efflux transporters found on the apical membrane of the enterocyte (Huang et al., 2006;Varma et al., 2011).However, suitable in vitro data defining the transporter kinetics of rosuvastatin in the enterocyte were unavailable.Therefore, we used an apparent permeability value from the colorectal adenocarcinoma (Caco-2) cell transwell assay to define its passive and active permeability across the intestinal wall.
In contrast, the role of canalicular efflux transporters is less well established.Using Madin-Darby canine kidney (MDCK II) cells expressing P-glycoprotein (P-gp) and membrane vesicles from mammalian cells expressing breast cancer resistant protein (BCRP), Huang et al. (2006) proposed that BCRP, but not P-gp, is important for the canalicular efflux of rosuvastatin.However, upon using double-transfected MDCK II cells, Kitamaru et al. (2008) demonstrated that P-gp is involved.His findings corroborated the results of Li et al. (2011) based on transporter knockdown Caco-2 cells.To circumvent the lack of precise, transport-er-specific data for canalicular efflux, we used an overall intrinsic clearance (CL int ) from sandwich-cultured human hepatocyte (SCHH) experiments to define the biliary excretion of rosuvastatin (Jones et al., 2012), as the SCHH model is able to delineate between overall sinusoidal uptake and canalicular efflux.
For renal transport, it has been found that organic anion transporter 3 (OAT3) is likely to be the main transporter responsible for the uptake of rosuvastatin into renal proximal tubular cells (PTC) at the basolateral membrane (Windass et al., 2007).Hence, transporter kinetics (J max and K m ) obtained from OAT3-expressing oocytes were used (Windass et al., 2007).On the other hand, BCRP, MRP2, and MRP4 are the transporters responsible for the efflux of rosuvastatin from the renal PTC at the apical membrane (Verhulst et al., 2008).Using a mixed culture of human PTC and distal tubular cells (DTC), the uptake of rosuvastatin at the basolateral membrane was found to be the rate-limiting process of renal clearance (Verhulst et al., 2008).Hence, we assumed that uptake transporter kinetics (OAT3) were equal to the efflux transporter kinetics in PTC cells.As BCRP and MRP2 are not included in Simcyp, the efflux kinetics were assigned to MRP4 in the model.

Fluvastatin model parameters
Fluvastatin's physicochemical characteristics, drug parameters, and source are provided in Table S2 1 .Compared to rosuvastatin, the passive permeability of fluvastatin in the Caco-2 assay is high (16.236x 10 -6 cm/s) (Lindahl et al., 2004).While transporters could be involved in the intestinal uptake of fluvastatin, insufficient in vitro data is available in the literature.
In contrast to rosuvastatin, the major route of elimination is metabolism, involving mainly CYP2C9 (largest contribution), CYP3A4, and CYP2C8, while CYP2D6 and CYP1A1 play minor roles (Fischer et al., 1999).The use of CL int values from recombinant human CYP450 (rhCYP) permits simultaneous modelling of enzyme activities in multiple organs (e.g., intestine, liver, and kidney).However, the sum of CL int from rhCYP does not correspond to that obtained from human liver microsomes (HLM).Therefore, CL int values obtained from rhCYP (Fischer et al., 1999) were further scaled to CL int from HLM experiments (Watanabe et al., 2010) using CYP abundance data in HLM (Nakamura et al., 2016).OATP1B1, OATP1B3, OATP2B1, NTCP, and MRP2 were found to be involved in the hepatic transport of fluvastatin (Kopplow et al., 2005;Noé et al., 2007;Bi et al., 2013;Ellis et al., 2013).Due to the lack of sufficient good-quality in vitro data for each specific transporter, CL int from SCHH experiments was used to inform the sinusoidal uptake and canalicular efflux in the PerL model (Jones et al., 2012;Izumi et al., 2018).Renal clearance was set as zero due to the lack of in vitro renal transporter data.
Following the latest modeling recommendations (Shebley et al., 2018), we evaluated our model by overlaying our simulations on the corresponding observed clinical data.The biokinetic parameters obtained from our simulations, such as area under the plasma concentration-time curve (AUC 0h-t ), maximum plasma concentration (C max ), plasma clearance (CL), and time to reach C max (T max ), were compared against those derived from clinical studies.The goodness-of-fit between simulated and observed biokinetic parameters was judged based on a two-fold difference success criterion.

Derivation of IVIVE scaling factors IVIVE for hepatic transporters
Typically, in vitro J max is scaled to in vivo J max using a relative expression or activity scalar to account for differences in transporter expression or activity between the in vitro expression system and primary HHEP (Hirano et al., 2004;Jamei et al., 2014).This is important as J max is affected by transporter protein expression and its scaling allows for accurate IVIVE.However, this fails to account for the decreased expression of transporters in isolated HHEP compared to human liver tissue that has been observed (Vildhede et al., 2014), and thus an additional scaling factor is required.Therefore, we defined a new relative expression factor (REF) as follows: (1) where ( 2) or ( 3) and (4) SF1 accounts for the difference in abundance or activity of the j th transporter between the isolated HHEP and the expression system, while SF2 accounts for the difference in abundance of the j th transporter between hepatocytes in liver tissue and in isolated HHEP (Fig. 1).
The units for transporter abundance in the liver and in HHEP are identical (pmol/mg membrane protein).Hence, SF2 is dimensionless.Expression of uptake and efflux transporters in pooled samples of fresh liver tissue and isolated hepatocyte membrane were obtained from the literature (Ohtsuki et al., 2011;Vildhede et al., 2015Vildhede et al., , 2016)).For SF1, depending on the expression system used, the units of transporter abundances differed and additional scaling factors were required.If membrane vesicles were used: (5) where transporter j abundance HHEP is the abundance of the j th transporter in HHEP (pmol/mg protein), transporter j abundance vesicles is the abundance of the j th transporter in the membrane vesicles (pmol/mg vesicular protein), f membrane is the frac-
Lactonization of pitavastatin by UGT1A3 and UGT2B7 is the major metabolic pathway of pitavastatin (Fujino et al., 2003).It is a reversible process, where ring-opening occurs non-enzymatically (Yamada et al., 2003;Fujino et al., 2004a).Metabolic data of the lactonization of pitavastatin by supersomes derived from insect cells infected with baculovirus-expressing human UGT1A3 and UGT2B7 were included in our model (Schirris et al., 2015).Pitavastatin is also metabolized by CYP2C9 and CYP2C8.Due to the lack of appropriate rhCYP metabolism data, CL int values from metabolism by HLM were used.To allow for the modelling of metabolism in organs other than the liver (intestine), we allocated the CL int value to CYP2C9, the major CYP450 enzyme involved (Fujino et al., 2003).This is because the model allows CYP450-specific kinetic data to be incorporated in multiple organs where that enzyme is expressed, while HLM kinetic data is confined to the liver.
Similar to rosuvastatin, active secretion of pitavastatin in the kidneys is mediated by OAT3 (Fujino et al., 2005).In vitro data from experiments conducted with kidney slices were used to inform the CL int for active renal uptake (Watanabe et al., 2011).As the K m for OAT3 (3.3 µM) was previously determined (Fujino et al., 2005), J max of OAT3 was back-calculated and included in our model (Results S1 1 ).Similar to our approach for rosuvastatin, the same kinetic parameters used for OAT3 were applied to MRP4 efflux along the apical membrane as basolateral uptake is the rate-limiting process in renal clearance (Watanabe et al., 2011).

Bottom-up model simulation and evaluation
Bottom-up simulations were performed for single IV and oral as well as repeated oral administration of all three statins.Clinical trial parameters, population characteristics, dose, and regimen were matched to the clinical studies.Details of the clinical studies are summarized in Table S4 1 .For each simulation performed, 10 trials were conducted.Dissolution characteristics of Crestor ® tablets were included for oral dosing simulations of rosuvastatin (Tab.S5 and S6 1 ).A solution formulation was used for fluvastatin and pitavastatin simulations.where HPPGL represents the number of hepatocytes per gram of liver and LiverWt is the liver weight of an individual (g).Post-scaling, in vivo J max has a unit of pmol/min.If Equation 6was used to obtain SF1, in vitro J max should be expressed as pmol/min/10 6 cells, while SF1 and SF2 are dimensionless.If Equation 5, 7, or 8 was used, in vitro J max should be expressed as pmol/min/mg protein, while SF1 has a unit of mg protein/10 6 HHEP and SF2 is dimensionless.If necessary, in vitro J max must be converted to the correct units using the total amount of protein per million HEK-293 cells (worked example in Results S2 1 ).Correct application of the scaling factors would result in the final unit of in vivo J max being pmol/min.Equation 9 was performed directly within Simcyp.

IVIVE for intestinal transporters
Similar to hepatic transporters, there is a different expression of transporters between the expression system and enterocytes.In our model, only the Caco-2 cell system was used to inform the kinetics of intestinal uptake transporters.To account for expression differences, Equation 10 was used: (10) where transporter n abundance Jejunal enterocytes, the abundance of the n th transporter in the jejunal enterocyte (pmol/mg protein), was compared with transporter n abundance Caco-2 cells, the abundance of the n th transporter in the Caco-2 cells (pmol/mg protein).Finally, to scale the J max of an uptake transporter from Caco-2 cells to in vivo intestinal CL int, the following equations were used2 .Firstly, J max and K m were converted to in vitro transporter-mediated appar-tion of membrane protein among total protein in a hepatocyte, f protein is the amount of total protein per million HHEP (mg protein/10 6 cells), and f inverted is the fraction of membrane vesicles that are active (inverted).Both f membrane and f inverted are assumed to be one-third and cancel out, while f protein is calculated to be 0.3235 mg protein/10 6 HHEP using data from Vildhede et al. (2016).When transfected HEK-293 cells were used, SF1 was obtained from Equation 6, 7, or 8 depending on whether relative expression ( 6), ( 7), or relative activity (8) was considered: (6) or ( 7) or ( 8) where transporter j abundance HHEP is the abundance of the j th transporter in HHEP (pmol/10 6 HHEP), transporter j abundance HEK-293 cells is the abundance of the j th transporter in the transfected HEK-293 cells normalized to total number of HEK-293 cells (pmol/10 6 HEK-293 cells), transporter j abundance HEK-293 protein is the abundance of the j th transporter in the transfected HEK-293 cells normalized to the total amount of protein (pmol/mg protein), CL int,HHEP is the intrinsic transporter clearance in HHEP (µL/min/10 6 HHEP) and CL int,HEK-293 is the intrinsic transporter clearance in HEK-293 cells (µL/min/mg protein).
Using the scaling factors obtained from Equations 1 to 8, in vivo J max was obtained from the following equation:  1) scaling from the in vitro expression system (membrane vesicles or HEK-293 cells) to isolated HHEP systems (SCHH or suspension hepatocytes) using SF1 and (2) scaling from isolated HHEP systems to hepatocytes in liver tissue using SF2.
resents the unbound intrinsic metabolic clearances for i metabolic pathways mediated by j isoforms of rhCYP (µL/min/pmol CYP j ), and CYP j abundance HLM represents the abundance of the j th CYP in HLM (pmol CYP j /mg protein).A single ISEF was calculated for fluvastatin metabolism and applied to CYP3A4, 2C9 and 2C8.For rosuvastatin and pitavastatin, default ISEF values from the Simcyp database were used.Finally, application of ISEF to scale in vitro metabolic kinetics to in vivo levels was performed in Simcyp and the equation is as follows: (16) Where CL int represents the in vivo intrinsic metabolic clearance in the liver (L/h), V max,i (rhCYP j ) represents the maximal metabolic rate for i metabolic pathways mediated by j isoforms of rhCYP (pmol/min/pmol CYP), K m,i (rhCYP j ) represents the Michaelis-Menten constant for i metabolic pathways mediated by j isoforms of rhCYP (µM) and MPPGL represents the amount of microsomal protein per gram of liver.

Optimization of Kp and middle-out simulations
Initial bottom-up simulations using the method of Rodgers and Rowland (2007) to estimate K p resulted in poor predictions of the V ss and suboptimal distribution profiles for rosuvastatin and pitavastatin.Therefore, we optimized organ K p values with data obtained from distribution experiments conducted with Sprague Dawley rats (Kimata et al., 1998;Nezasa et al., 2002).Subsequently, middle-out simulations were performed to compare the results with the bottom-up simulations to evaluate the degree to which simulations are improved by judicious incorporation of in vivo data.
This was performed for rosuvastatin and pitavastatin but not for fluvastatin, as the tissue concentration (based on radioactivity) could not be accurately attributed to unchanged fluvastatin owing to its extensive metabolism.Conversely, rosuvastatin undergoes minimal metabolism, while we assumed pitavastatin metabolism was confined to the liver and gut.Hence, K p values were calculated for both compounds by using the ratio of drug concentration between a particular tissue and plasma.Where multiple K p values were computed because of multiple sampling times, we chose to use the largest K p value.

REF for the hepatic transport of rosuvastatin, fluvastatin, and pitavastatin
The exact values, derivation, and references used for SF1, SF2, and REF are summarized in Tables S7-S13 1 .For uptake transporters expressed in transfected HEK-293, SF1 varied from 1.571-121 for rosuvastatin (Tab.S7 1 ) and 0.155-1.790for pitavastatin (Tab.S8 1 ).For the same transporter, a large variation in the relative expression of uptake transporters in HEK-293 cells versus hepatocytes was observed, leading to highly variable SF1 values.Among the uptake transporters, SF1 for ent permeability for each i th segment of the intestine: (11) where P app,trans,n,i represents apparent permeability contributed by the n th transporter in the i th segment of the intestine (10 -6 cm/s), J max is the maximal transport rate of the n th transporter (pmol/min), A represents the transwell surface area (cm 2 ), K m is the Michaelis-Menten constant (µM), fu inc represents the unbound fraction in the incubation system and C lumen,i represents the concentration of substrate in the lumen of the i th segment (µM).If J max is normalized to the total amount of protein (we refer to as J max '), an additional unit conversion was applied to express it as pmol/min: (12) where J max ' (pmol/min/mg protein) is multiplied by the total protein per transwell filter (mg) to obtain J max .
Next, in vitro P app,trans,n,i is converted to in vivo effective active permeability via a regression equation obtained from Loc-I-Gut jejunal permeability and Caco-2 studies (Sun et al., 2002): (13) where P eff,trans,n,i represents in vivo effective active permeability contributed by the n th transporter in the i th segment (10 -4 cm/s), while the regression coefficients are A = 0.939 and B = -0.8787(default values in Simcyp).Finally, P eff,trans,n,i is converted to uptake clearance for each n th transporter in the i th segment: (14) where CL int,T,i represents the uptake clearance for each i th segment of the intestine (L/h), SA i represents the surface area of the i th segment (m 2 ), RSA n represents the relative segmental abundance of the n th transporter compared to that of the jejunum I segment.In ADAM, the intestines are segregated into 8 segments: duodenum I, jejunum I and II, ileum I, II, III and IV, and colon.Each segment has its own CL int , depending on the particular transporter expression level and surface area within each segment.Equations 11, 13, and 14 were performed within Simcyp.

IVIVE for metabolism
To account for differences in the intrinsic activity of enzymes between rhCYP and that measured in liver microsomes, we adapted the intersystem extrapolation factor (ISEF) from Procter et al. ( 2004) to obtain Equation 15: (15) where CL int,u (HLM) represents the unbound intrinsic metabolic clearance in HLM (µL/min/mg protein), CL int,i,u (rhCYP j ) rep-OATP2B1 exhibited the greatest variability, with a value of (for rosuvastatin) from Bosgra et al. ( 2014) compared to 0.155 (for pitavastatin) from Hirano et al. (2006).In other words, there is significant inter-lab variability in the degree of transporter expression in transfected cells.In contrast, for efflux transporters expressed in membrane vesicles, a clear overexpression in the membrane vesicles compared to the isolated hepatocytes was observed, with SF1 ranging from 0.000327 to 0.0105 (Tab.S7, S8 1 ).
Separately, the values of SF2 varied from 1.230 to 5.114, indicating a consistent underexpression of both uptake and efflux transporters in isolated hepatocytes compared to liver tissue (Tab.S9, S10 1 ).Upon multiplying SF1 by SF2, REF varied from 0.190-148.830for uptake transporters and 0.00052 to 0.028 for efflux transporters (Tab.S11-S13 1 ), indicating that the directionality of adjustment of J max values is transporter-and expression system-specific.

REF for the intestinal transport of pitavastatin
In contrast to the underexpression of hepatic uptake transporters in the in vitro systems used, Caco-2 cells exhibited an overexpression of OATP2B1 as compared with jejunal epithelial cells.A REF value of 0.0625 was obtained from the same study from which intestinal OATP2B1 kinetic data was obtained (Ölander et al., 2016).Detailed calculations are shown in Results S3 1 .

ISEF for the metabolism of fluvastatin
ISEF for CYP3A4, 2C9, and 2C8 calculated using Equation indicated that metabolic clearance required additional scaling by 2.444-fold.Detailed calculations are shown in Results S4 1 .

Effect of scaling factors on bottom-up IV simulations
To demonstrate the importance of the ISEF and REF (SF1 and SF2), IV infusions of rosuvastatin and pitavastatin were simulated with and without scaling factors.In the absence of scaling factors, simulated AUC 0h-t exhibited a fold difference of 10.25 and 5.35 compared to the observed values for rosuvastatin and pitavastatin, respectively (Tab.1).Significant underestimation of the CL was observed as well.Upon applying the relevant scaling factors, simulated biokinetic parameters improved to within two-fold of the observed values for both simulations (Tab.2).As   ( Martin et al., 2002Martin et al., , 2003)).The simulated biokinetic parameters with their geometric means overlaid with the observed data and the fold difference are listed in Table 2.For all 3 dosing scenarios, the C max , T max , AUC 0h-t , AUC 0h-∞ , and CL were within 1.5-fold of the observed data.However, the model was unable to recapitulate the triphasic decline observed in intravenous dosing of rosuvastatin (Fig. 3A).

Bottom-up simulations of single and multiple dose biokinetics of fluvastatin
As shown in Table 2, the predicted C max , T max , AUC 0h-t , AUC 0h-∞ , and CL for the 3 dosing scenarios fell within two-fold expected, the simulated plasma concentration-time profiles improved upon application of the scaling factors (Fig. 2).

Bottom-up simulations of single and multiple dose biokinetics of rosuvastatin
Having established that the derivation and application of mechanistic scaling factors permitted reasonable IVIVE for transporter-dependent, metabolism-dependent and mixed-mode statins, we proceeded to simulate various dosing regimens for each compound.Figure 3 shows the simulated plasma concentration-time profile of rosuvastatin after a single 8 mg IV infusion, a single 20 mg oral dose, and a 10 mg once daily dosing for 14 days  eters within 2-fold of the observed data (FDA, 2012) (Tab.2).However, the 5-day repeated dosing was sub-optimal.Simulated parameters fell outside the two-fold criterion, with C max and AUC 0h-t being underpredicted while CL and T max were overpredicted.This is in contrast to the predicted concentration-time profile, where the terminal gradients for all 3 dosing regimens were poorly predicted.The predicted terminal gradient was much steeper than the observed clinical data.

Effect of rat Kp on rosuvastatin and pitavastatin simulations
To investigate the importance of rat K p to our models, the same simulations mentioned previously were performed with and without the use of rat K p data for rosuvastatin and pitavastatin.Initially, the bottom-up model underestimated V ss with a of the observed parameters (FDA, 2012).Only the T max for the 10 mg single oral dosing was beyond the two-fold criteria with a fold difference of 2.35 versus the observed value.Similarly, the plasma concentration-time profile (Fig. 4) shows that the simulations are in close agreement with the observed profile.Furthermore, the nonlinear biokinetics observed during repeated 40 mg oral dosing of fluvastatin due to saturable first-pass elimination were recapitulated in our model (Tse et al., 1992).

Bottom-up simulations of single and multiple dose biokinetics of pitavastatin
The plasma concentration-time profiles of pitavastatin after a 2 mg IV infusion over 1 hour, 2 mg single oral dosing and 4 mg once daily dosing for 5 days are shown in Figure 5. Simulated biokinetic parameters for the IV infusion and single oral dose were well predicted, with the simulated biokinetic param- of disposition processes remains a major barrier to its routine use in risk assessment (Wambaugh et al., 2018).Building upon the work by Lundquist, Vildhede and colleagues, we show that the key to greater accuracy in IVIVE of transporters lies in the use of quantitative proteomics to account for differential expression of transporters in in vitro systems compared with isolated hepatocytes and liver tissue (Kimoto et al., 2012;Lundquist et al., 2014;Vildhede et al., 2015Vildhede et al., , 2018)).The inclusion of an additional scaling factor (SF2) that accounted for this downregulation permitted predictions of CL, AUC 0h-t , and C max to within two-fold of the observed data after an IV infusion or single oral administration of the 3 compounds.This is notable for the following reasons: Firstly, this demonstrates the ability of bottom-up PBK modeling to accurately recapitulate biokinetics in humans when comprehensive in vitro data and accurate scaling factors are available.Secondly, our results confirm our previous hypothesis that the suboptimal IVIVE for transporter-mediated elimination is due to the downregulation of transporters when hepatocytes are isolated from liver tissue.Thirdly, we illustrate that this mechanistic scaling approach is sufficiently robust to accommodate data generated from a variety of in vitro systems across different laboratories, measured and reported using different approaches.Our promising findings fuel further simulations using compounds with different biokinetic profiles to confirm their robustness.
In our derivation of SF2, data of transporter protein expression in liver and hepatocytes were obtained from three publications (Ohtsuki et al., 2011;Vildhede et al., 2015Vildhede et al., , 2016)).However, there was an inter-laboratory fold difference of 1.8-90 in the predicted value of 0.12 l/kg and 0.091 l/kg compared to the observed values of 1.73 l/kg and 1.90 l/kg for rosuvastatin and pitavastatin.Simulated values were underpredicted by more than tenfold compared to the observed values.Thus, optimizations of K p using data from rat distribution studies were performed and their references are shown in Tables S14-S16 1 .This led to an improved value of 0.62 l/kg and 0.45 l/kg for rosuvastatin and pitavastatin.
However, biokinetic parameters with and without using rat K p were relatively similar and fell within two-fold of the observed data even without the use of rat K p (Tab. 3).This is in contrast to the plasma concentration-time profile (Fig. 6).With the use of rat K p , the terminal gradient for pitavastatin simulations was predicted accurately as compared to the simulated profile without the use of rat K p (Fig. 6D-F).Our simulations demonstrated a short and negligible distribution phase with the majority of the drug being cleared within 10 hours.Upon including rat K p , the distribution profile for our simulations improved for pitavastatin.However, for rosuvastatin, the inclusion of rat K p did not lead to a significant change in the simulated plasma concentration-time profile (Fig. 6A-C).

Discussion
In this study, we demonstrate the utility of bottom-up PBK modeling as an alternative to animal testing to predict human biokinetic data.Bottom-up IVIVE is considered vital for prediction of chemical exposure, but the systematic suboptimal estimation quent reconversion of lactone to parent drug would increase systemic pitavastatin exposure.Further studies to investigate this reconversion are needed to improve IVIVE for repeated dosing of pitavastatin.Secondly, a pitavastatin tablet was administered in the original clinical trial.However, we did not include any dissolution data and chose to use a solution formulation instead because appropriate dissolution data for pitavastatin are unavailable.Finally, in a clinical study that reported repeated dosing of 2 mg pitavastatin over 6 days, the observed V ss and half-life (t 1/2 ) increased two-fold compared to a single dose of pitavastatin (Luo et al., 2015).This was speculated by the authors to be a result of saturable distribution (in this case of the main peripheral tissue, the liver) during multiple dosing.This phenomenon was not included in our model and could also contribute to the overprediction of CL.
While we have attempted to be as mechanistic as possible, we were unable to include active intestinal transporters in our models for rosuvastatin and fluvastatin due to a lack of appropriate in vitro data.Nonetheless, we were able to recapitulate the kinetics of rosuvastatin without the inclusion of any active intestinal transport.This could indicate that active transport is not significant in the oral absorption of rosuvastatin.However, for fluvastatin, the T max after a single oral dosing was overestimated and above two-fold of the observed T max .To fully recover the absorption kinetics of fluvastatin, further improvements of our model could be made by including the kinetics of intestinal transporters.Moreover, in our use of the mechanistic kidney model, we did not consider the possible differences in expression of transporters between isolated human tubular epithelial cells and kidney tissue (Results S1 and S5 1 ).A phenomenon similar to the downregulation observed in isolated hepatocytes leading to an underestimation of the renal clearance may have occurred.This could explain the slight underprediction of CL for rosuvastatin and pitavastatin after IV infusion.
Our simulations reveal that there is a need for better predictions of V ss using in vitro data and physicochemical properties.When our rosuvastatin and pitavastatin bottom-up models incorporated rat K p values to optimize distribution, this resulted in a clear improvement in the predicted V ss values.In contrast, the effects of rat K p on the simulated plasma concentration-time profile were mixed.For rosuvastatin, the profiles for all 3 dosing regimens did not differ much upon the inclusion of rat K p , and were relatively well predicted even without the use of rat K p .Notably, even with the use of rat K p , our simulations were unable to fully recapitulate the triphasic decline of rosuvastatin observed after an IV infusion (Martin et al., 2003).This example demonstrates the challenges that current distribution models face when predicting more complex distribution scenarios.Interestingly, in the case of pitavastatin, although incorporating rat K p allowed us to fully recover the terminal gradient of the plasma concentration-time profile for all 3 dosing regimens, this did not materially alter our prediction of biokinetic parameters.In other words, where animal distribution data are available, distribution profiles can be further improved, but the absence of such data does not appear to hinder modelling reliability.
Finally, even though active transport is involved in the tissue distribution of a compound (Grover and Benet, 2009), the Rod-in vivo uptake clearance when measured protein levels from different laboratories were used for modeling (Wegler et al., 2017).Hence, there is a need for a comprehensive meta-analysis of all available transporter abundance data in the liver and isolated hepatocytes to obtain appropriate values for scaling.This would mitigate variations caused by different proteomic techniques or different donor characteristics when calculating SF2.
In our calculation of SF1, we found a large variability of expression levels of the same transporter in the same in vitro systems reported in different references.This exemplifies the need for every laboratory to carry out its own quantitative measurements of the expression of transporters in its own particular in vitro system.The observed variability could be due to inter-laboratory differences in gene transfection protocols when using HEK-293 cells (Ooi et al., 2016) or the quantitative proteomics method used to measure the expression of transporters (Wegler et al., 2017).Similarly, when using Caco-2 cells, cell density per unit surface area (mg/cm 2 ) would be heavily influenced by cell culture techniques, area of insert used to culture the cells, and growth conditions (Sambuy et al., 2005).Without such information, inaccurate estimates of intestinal transporter kinetics would be made when scaling using Equation 9. Our analyses indicate that it is not appropriate to cross-apply scaling factors across different laboratories and in vitro systems.
Mechanistic modeling of metabolism and transporter kinetics was performed to describe the biokinetics of the three statins.For bottom-up IVIVE, the quality and appropriateness of the data are critical.In the case of transporter kinetics, it is crucial to consider: (1) type of in vitro system, (2) unit in which the data is presented, (3) availability of quantitative measurements of transporter abundances, and (4) surface area of and cell density on the transwell filter.There is often a significant overexpression of transporters in membrane vesicles, which would result in an overestimation of the contributions of the transporter without a correction based on quantitative measurements of the transporter expression.For example, expression of breast cancer resistant protein (BCRP) was found to be 3000 times greater in membrane vesicles compared to isolated HHEP, leading to an SF1 of 0.000327 (Vildhede et al., 2016).Furthermore, the unit of J max is typically expressed as pmol/min/mg protein when using membrane vesicles, which has to be converted to pmol/min/10 6 cells using transporter expression data in HHEP for subsequent scaling.Such considerations are also pertinent for other in vitro systems (HEK-293 cells, Caco-2 cells, etc.).These observations underscore the need for careful scrutiny of the experimental units and derivation of scaling factors using quantitative measurements (Prasad and Unadkat, 2014;Vildhede et al., 2018).
Interestingly, although our pitavastatin model recapitulated biokinetic profiles of both IV infusion and single oral dosing, we were unable to fully recapitulate the biokinetics of repeated pitavastatin dosing.Analysis of the model suggested that this was due to an overestimation of CL, which could be explained by several reasons: Firstly, non-enzymatic reconversion of pitavastatin lactone (a metabolite) back to pitavastatin has been reported (Fujino et al., 2004b(Fujino et al., , 2005)).Given that pitavastatin lactone is the major metabolite of pitavastatin with a half-life (t 1/2 ) of 11-16 hours (FDA, 2012), it is possible that an accumulation and subse-gers and Rowland equation used to predict K p values does not account for any active transporter process that contributes to the transport of xenobiotics into and out of tissue (Rodgers and Rowland, 2007).Also, it has been previously noted that predictions of K p and V ss using the same method are less reliable when a compound is involved in active transport processes (Rodgers and Rowland, 2007).Furthermore, if EHC occurs, it may increase the mean residence time of the compound, leading to a larger V ss (Chan et al., 2018).Since the dispositions of rosuvastatin and pitavastatin are greatly influenced by transport processes and have been noted to undergo EHC (Kimata et al., 1998;Martin et al., 2003;FDA, 2012), this could explain our sub-optimal prediction of distribution profiles for rosuvastatin and pitavastatin but not for fluvastatin.Therefore, further research is needed to improve the predictions of K p and V ss , especially for xenobiotics that rely on active transport for tissue distribution.

Conclusion
In conclusion, a bottom-up PBK model of rosuvastatin, fluvastatin, and pitavastatin was developed using mainly in vitro data, QSAR models, and animal distribution data.The IVIVE of transporter-mediated disposition was greatly improved upon the application of additional scaling factors that account for differential expression of transporter proteins between in vitro systems and in vivo organs.Our novel findings suggest that bottom-up PBK modeling is a viable alternative to animal testing for predicting human biokinetics of xenobiotics.Further research must be performed to validate this bottom-up PBK-IVIVE approach and develop new equations to better predict the distribution of xenobiotics in the body.In summary, PBK modeling has the potential to replace animal testing in predicting biokinetics of xenobiotics in humans.

Fig. 1 :
Fig. 1: Scaling of transporter expression using SF1 and SF2Jmax in the liver is estimated by (1) scaling from the in vitro expression system (membrane vesicles or HEK-293 cells) to isolated HHEP systems (SCHH or suspension hepatocytes) using SF1 and (2) scaling from isolated HHEP systems to hepatocytes in liver tissue using SF2.

Fig. 2 :
Fig. 2: Plasma concentration-time profile of simulations performed without scaling factors (dashed line) and with scaling factors (continuous line) for IV infusions of (A) 8 mg rosuvastatin and (B) 2 mg pitavastatin The blue boxes represent the observed rosuvastatin and pitavastatin concentration with error bars (SD) for observed data where available.

Fig. 3 :
Fig. 3: Bottom-up simulated rosuvastatin model and observed plasma concentration-time profile of rosuvastatin after (A) 8 mg IV infusion, (B) 20 mg single oral dose, and (C) 10 mg once daily oral dosing for 14 days The continuous line represents the predicted mean concentration for the simulated population.The dotted and dashed lines represent the 5 th and 95 th percentile of the predicted mean concentration.The blue boxes represent the observed rosuvastatin concentration with error bars (SD) if available.

Fig. 4 :
Fig. 4: Bottom-up simulated fluvastatin model and observed plasma concentration-time profile of fluvastatin after (A) 2 mg IV infusion, (B) 10 mg single oral dose, and (C) 40 mg once daily oral dosing for 6 days The continuous line represents the predicted mean concentration for the simulated population.The dotted and dashed lines represent the 5 th and 95 th percentile of the predicted mean concentration.The blue boxes represent the observed fluvastatin concentration with error bars (SD) if available.

Fig. 5 :
Fig. 5: Bottom-up simulated pitavastatin model and observed plasma concentration-time profile of pitavastatin after (A) 2 mg IV infusion, (B) 2 mg single oral dose, and (C) 4 mg once daily oral dosing for 5 days The continuous line represents the predicted mean concentration for the simulated population.The dotted and dashed lines represent the 5 th and 95 th percentile of the predicted mean concentration.The blue boxes represent the observed pitavastatin concentration.

Fig
Fig. 6: Middle-out simulated and observed plasma concentration-time profile for rosuvastatin after (A) 8 mg IV infusion, (B) 20 mg single oral dose, (C) 10 mg once daily oral dosing for 14 days, and for pitavastatin after (D) 2 mg IV infusion, (E) 2 mg single oral dose, and (F) 4 mg once daily oral dosing for 5 days The continuous and dashed lines represent the predicted mean concentration for the simulated population using rat K p and without using rat K p , respectively.The blue boxes represent the observed rosuvastatin or pitavastatin concentration with error bars (SD) if available.