Physiologically-based pharmacokinetic model-based translation of OATP1B-mediated drug-drug interactions from coproporphyrin I to probe drugs.

The accurate prediction of OATP1B-mediated drug-drug interactions (DDIs) is challenging for drug development. Here, we report a physiologically-based pharmacokinetic (PBPK) model analysis for clinical DDI data generated in heathy subjects who received oral doses of cyclosporin A (CysA; 20 and 75 mg) as an OATP1B inhibitor, and the probe drugs (pitavastatin, rosuvastatin, and valsartan). PBPK models of CysA and probe compounds were combined assuming inhibition of hepatic uptake of endogenous coproporphyrin I (CP-I) by CysA. In vivo Ki of unbound CysA for OATP1B (Ki,OATP1B ), and the overall intrinsic hepatic clearance per body weight of CP-I (CLint,all,unit ) were optimized to account for the CP-I data (Ki,OATP1B , 0.536 ± 0.041 nM; CLint,all,unit , 41.9 ± 4.3 L/h/kg). DDI simulation using Ki,OATP1B reproduced the dose-dependent effect of CysA (20 and 75 mg) and the dosing interval (1 and 3 h) on the time profiles of blood concentrations of pitavastatin and rosuvastatin, but DDI simulation using in vitro Ki,OATP1B failed. The Cluster Gauss-Newton method was used to conduct parameter optimization using 1000 initial parameter sets for the seven pharmacokinetic parameters of CP-I (β, CLint, all , Fa Fg , Rdif , fbile , fsyn , and vsyn ), and Ki,OATP1B and Ki,MRP2 of CysA. Based on the accepted 546 parameter sets, the range of CLint, all and Ki,OATP1B was narrowed, with coefficients of variation of 12.4% and 11.5%, respectively, indicating that these parameters were practically identifiable. These results suggest that PBPK model analysis of CP-I is a promising translational approach to predict OATP1B-mediated DDIs in drug development.

area under the concentration time curve

area under the concentration time curve ratio (rifampicin/control)

breast cancer resistance protein

Cluster Gauss–Newton method

maximum concentration

coefficient of variation

cyclosporin A

drug–drug interaction

inhibition constant

multidrug resistance protein 2

organic anion transporting polypeptide 1B1

organic anion transporting polypeptide 1B3

physiologically‐based pharmacokinetic

time to maximum concentration

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Physiologically‐based pharmacokinetic (PBPK) models are used to predict transporter‐mediated drug–drug interactions (DDIs). Endogenous OATP1B biomarkers, such as coproporphyrin I (CP‐I), are strongly predicted to improve DDI prediction in drug development.

WHAT QUESTION DID THIS STUDY ADDRESS?

Can PBPK model analysis of the clinical CP‐I data successfully predict OATP1B‐mediated DDIs using probe drugs?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

The value of the most important DDI parameter, Ki,OATP1B, estimated by PBPK model‐based analysis of clinical CP‐I data, was able to overcome the discrepancy between the in vitro and in vivo effects of CysA on OATP1B, and could be applied to predict the change in the blood concentration time profiles of OATP1B probe drugs.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

The collection of endogenous OATP1B biomarker data is a feasible strategy to capture DDI potential. PBPK models aids in the prediction of its clinical impact more precisely, help to reduce risk in drug development, and impact the regulatory decision tree for DDI risk assessment.

INTRODUCTION

OATP1B1 and OATP1B3 (OATP1B) are multispecific transporters mediating the hepatic uptake of various drugs from blood circulation into hepatocytes. The hepatic clearance of OATP1B substrates follows the extended clearance concept, where overall intrinsic hepatic clearance is described by intrinsic parameters for sinusoidal uptake mediated by OATP1B, sinusoidal efflux and metabolism, and canalicular efflux. An OATP1B‐mediated process governs the first step of the overall hepatic elimination for its substrates, and, hence, changes in OATP1B activity may alter the systemic exposure of the substrate compounds when the liver is the major clearance organ in the body.

The physiologically‐based pharmacokinetic (PBPK) model considers both physiological parameters (tissue volumes and blood flow rates), and anatomic characteristics (organ layout in the body). The PBPK model is essential in drug development for analyzing and predicting the impact of drug–drug interactions (DDIs) of new investigational drugs as victims and perpetrator. , We have developed PBPK models to account for the plasma and tissue concentration time profiles of an OATP1B substrate, pravastatin, for the first time, and applied them to the quantitative analysis of the OATP1B‐mediated DDIs caused by rifampicin and cyclosporin A (CysA). The PBPK model for OATP1B substrates is characterized by the liver model, which consists of capillary and parenchymal compartments, and is represented by five tandemly connected units along the blood flow, so that the investigator can use the same PBPK model for the analysis of even high clearance compounds. , After optimization of the parameters, PBPK models of the victim and perpetrator compounds, which were combined assuming a competitive inhibition of OATP1B by the perpetrator compounds, could reasonably explain the alteration in the time profiles of the plasma concentrations of the OATP1B substrate drugs and endogenous substrates. , A discrepancy was observed in the values for the inhibition constant (Ki) of unbound CysA for OATP1B (Ki,OATP1B) estimated using the PBPK models and the values determined in vitro (1.6 nM vs. 26 nM, respectively), and OATP1B‐mediated DDI caused by CysA was underestimated using in vitro Ki,OATP1B.

The endogenous substrates of OATP1B, such as coproporphyrin I (CP‐I), is considered as functional biomarkers to detect changes of OATP1B activity by OATP1B inhibitor drugs, when the administration of probe drugs is not practical. , CP‐I has also been used to investigate genetic variations in OATP1B1*15 , , and the effect of disease conditions on OATP1B activities. , Significant effect of OATP1B1 genotypes on the plasma CP‐I concentrations supports major contribution of OATP1B1 to the hepatic uptake of CP‐I. Likewise, after administration of an inhibitor of OATP1B1/1B3, the plasma concentration of CP‐I displayed a transient increase, and the magnitude of change in the maximum plasma concentrations (Cmax) and the area under the plasma concentration time curve (AUC) depend on the magnitude of OATP1B inhibition during the study. , , The PBPK model‐based analysis of CP‐I together with the concentration time profiles of rifampicin could yield a converged Ki,OATP1B, , which, in turn, is anticipated to help DDI simulation with OATP1B substrate drugs and reduce the uncertainty involved in extrapolating in vitro Ki,OATP1B to in vivo Ki,OATP1B.

Recently, we conducted a clinical DDI study using CysA as an OATP1B inhibitor where CysA was given orally at 20 and 75 mg with a 1‐h interval before probe drug administration, and at 75 mg with a 3‐h interval before probe drug administration (Mochizuki et al., submitted). PBPK‐model based analysis of the blood concentrations of CP‐I and probe drugs (pitavastatin and rosuvastatin) determined in the same subjects was conducted to support the feasibility of the translation of CP‐I data to the DDI prediction with actual probe drugs. The PBPK model‐based analysis of clinical data essentially includes optimization of the selected parameters, such as Ki,OATP1B, and the initial parameter estimation is often laborious and has some uncertainty. To address this issue, Aoki et al. developed Cluster Gauss–Newton method (CGNM), a comprehensive parameter finding the algorithm to yield multiple parameter sets that satisfy the observed values. Our first application of this newly developed method was to analyze the nonlinear pharmacokinetics of bosentan, an endothelin receptor antagonist. Using the PBPK model, saturation of hepatic uptake and target binding were considered. In the present study, the CGNM was also applied to the PBPK‐model analysis of the OATP1B‐mediated DDIs caused by CysA to elucidate parameter sets that satisfy the CP‐I data.

METHODS

Clinical study design

The data were cited from our clinical DDI study registered as the specified clinical trial in Japan Registry of Clinical Trials (https://jrct.niph.go.jp/en‐latest‐detail/jRCTs031200012). The subjects in the study were 10 healthy men recruited in Japan: aged 22–39 years, height 164–182 cm, and weight 53–76 kg. Neoral was given orally at 9 or 7 a.m. for a 1‐h or 3‐h dosing interval study, respectively. Pitavastatin (0.2 mg), rosuvastatin (1 mg), and valsartan (2 mg) were administered in a cocktail as test OATP1B probe. These drugs were administered under fasting conditions.

PBPK model analysis and parameters

The structures of the PBPK models for CysA and pitavastatin, CP‐I, and rosuvastatin are reported as documented without any modification, whereas some parameters in the models were optimized for the data obtained in this study. The fixed and optimized values of each parameter for CysA and OATP1B probes are summarized in Tables 1 and 2. Most parameters are cited from previous reports , unless otherwise noted. Body weight (kg), the blood flow rates (L/h/kg), and tissue volumes (L/kg) are as follows; 64.5 kg, Qliver of 1.24, Qmuscle of 0.642, Qskin of 0.257, and Qadipose of 0.223, and Vcentral of 0.0743, VHC of 0.0174, VHE of 0.0067, Vmuscle of 0.429, Vskin of 0.111, and Vadipose of 0.143, respectively.

TABLE 1

Pharmacokinetic parameters for PBPK model, and Ki of CysA

Parameter	Value	Units	Comments
FaFg	0.337 ± 0.017 (20 mg), 0.572 (75 mg)		Optimized/Previous report 5
Kp,liver	25.3 ± 2.3		Optimized
fhCLint	1.08 ± 0.08	L/h	Optimized
fp	0.0242		Geometric mean of previous reports c
Rb	2.15		Geometric mean of previous reports d
fb	0.0113		fp/Rb
Tlag	0.254	h	Previous report 5
ka	0.999	/h	Previous report 5
Xps	2.20 ± 0.27		Optimized
PSm	517		Previous report 5
PSs	78.9		Previous report 5
PSa	21.5		Previous report 5
Yft	0.0387 ± 0.0038		Previous report 5
fm	0.116		Previous report 5
fs	0.0255		Previous report 5
fa	0.0201		Previous report 5
Ki,OATP1B a	0.536 ± 0.041	nM	Optimized using CP‐I data
Ki,MRP2 b	13.3	μM	Geometric mean of previous reports e

Abbreviations: FaFg, product of fraction absorbed and availability in the intestine; fb, unbound fraction in the blood; fhCLint, product of unbound fraction in the liver and intrinsic clearance; fp, unbound fraction in the plasma; ka, absorption rate constant; Kp,liver, liver‐to‐plasma concentration ratio; PBPK, physiologically‐based pharmacokinetic; Rb, blood‐to‐plasma ratio; Tlag, lagtime for intestinal absorption; Xps, scaling factor for the clearance of membrane transport in the peripheral organs (common value across the peripheral tissues); PSm, PSs, PSa, membrane transport in the muscle (m), skin (s) and adipose (a), respectively (included as Xps × PSm, Xps × PSs, and Xps × PSa in the model); Yft, scaling factor for the unbound fraction in the peripheral organs; fm, fs, fa, unbound fraction in the muscle (m), skin (s) and adipose (a), respectively (included as Yft × fm, Yft × fs, and Yft × fa in the model).

Defined with regard to the unbound CysA concentration in the blood.

Defined with regard to the unbound CysA concentration in the liver.

References 30, 31.

References 5, 30, 32.

References 33, 34, 35.

TABLE 2

Pharmacokinetic parameters for PBPK modeling of CP‐I, pitavastatin and rosuvastatin

Parameter	Values			Units	Comments
	CP‐I	Pitavastatin	Rosuvastatin
CLint,all,unit	41.9 ± 4.3	80.8 ± 6.9	6.00 ± 0.55	L/h/kg	Optimized
Kp,muscle	0.103	0.113	0.144		Previous report a
ka	3	0.775	0.0546 ± 0.0074	/h	Previous report a (CP‐I) Optimized by previous data 16 (pitavastatin) Optimized b (rosuvastatin)
ktransit	5.15	0.403	2.73	/h	Previous report a (CP‐I and rosuvastatin) Optimized by previous data 16
fh	0.0697	0.035	0.179		Previous report a
fbile	0.844	0.406	0.943		Previous report a (CP‐I and rosuvastatin) Optimized by previous data 16
Kp,adipose	0.079	0.086	0.087		Previous report a
fb	0.0105	0.009	0.174		Previous report a
Rb	0.628	0.578	0.69		Previous report a
FaFg	0.309	1	0.429		Previous report a
Kp, skin	0.442	0.481	0.439		Previous report a
Rdif	0.0352	0.0345	0.00502		Previous report a
γ	0.02	0.244	0.25		Previous report a
Baseline value	0.247	–	–	μg/L	Mean of observed blood concentration in control
β	0.5	0.5	0.5		Fixed
CLr,unit	0.0421	0	0.201	L/h/kg	Geometric mean of observed CLr, unit in control phase
Molecular weight	655	421	482
fsyn	1	–	–		Previous report a
Tlag	–	0	0	h	Previous report a
kstomach	–	–	0.873	/h	Previous report a

Abbreviations: CLint,all,unit, overall intrinsic hepatic clearance per unit body weight; Kp,muscle, Kp,adipose, Kp, skin, tissue‐to‐plasma concentration ratio; ka, rate constant for absorption; ktransit, rate constant for transit from the bile compartment to the intestine; fh, unbound fraction in the liver; fbile, fraction of biliary excretion, (CLint,bile/[CLint,bile + CLint,met]); fb, unbound fraction in the blood; Rb, blood‐to‐plasma ratio; FaFg, product of fraction absorbed and availability in the intestine; Rdif, PSdif,inf/PSact,inf (passive clearance/active transport); γ, PSdif,inf/PSeff (clearance ratio for passive influx to the sinusoidal efflux); β, (CLint,met + CLint,bile)/(PSeff + CLint,met + CLint,bile); CLr,unit, renal clearance per unit body weight; fsyn, fraction of synthesis in the liver in the body to the synthesis in the whole body; Tlag, lag time for drug absorption of pitavastatin and rosuvastatin; kstomach, rate constant for the transit from the stomach to the intestine.

Parameters were cited from previous report; CP‐I ref. 6, pitavastatin ref. 5, rosuvastatin ref. 21.

ka was optimized to account for the plasma concentration time profiles in the control phase.

Previous report (CP‐I)

Optimized by previous data (pitavastatin)

Optimized (rosuvastatin)

Previous report (CP‐I and rosuvastatin)

Optimized by previous data

Previous report (CP‐I and rosuvastatin)

Optimized by previous data

To reproduce the observed blood concentrations of CysA, tissue‐to‐blood concentration ratio (Kp,liver) and the product of unbound fraction in the liver and intrinsic hepatic clearance (fhCLint) were optimized using all the observed data (20 mg, 75 mg [−1 h] and 75 mg [−3 h]). FaFg was also set as a floating parameter only in the condition of CysA 20 mg (−1 h). Then, the PBPK model of CysA was combined with that of CP‐I, assuming the competitive inhibition of OATP1B as described previously. The overall intrinsic hepatic clearance per body weight (CLint,all,unit) was optimized for CP‐I using all of the observed clinical data (baseline, CysA 20 mg, 75 mg [−1 h] and 75 mg [−3 h]). The basal blood concentration of CP‐I was defined as the average of the CP‐I blood concentration in the control phase. The renal clearance of CP‐I represents the geometric mean of the observed renal clearance per body weight (CLr, unit) in the control phase. The value of β, ratio of the overall intrinsic hepatic clearance to the hepatic uptake clearance, was fixed as 0.5, because this parameter could not be identified uniquely in the previous analysis, and its absolute value did not affect the estimation of free Ki,OATP1B.

To simulate the effect of CysA on DDIs, the PBPK models of CysA and pitavastatin or rosuvastatin were combined by assuming the competitive inhibition of OATP1B as described. The CLint,all,unit of pitavastatin and rosuvastatin, and absorption rate constant (ka) of rosuvastatin were optimized by the nonlinear regression method using the control data. The Ki,OATP1B optimized for CP‐I data was used in the calculation, or optimized to account for the effect of CysA on the blood concentration time profiles of pitavastatin and rosuvastatin.

Estimation of the distribution of parameters using the Cluster Gauss–Newton method

In the PBPK model of CP‐I, nine parameters were estimated using the CGNM which samples multiple solutions of nonlinear least‐squares problems. The CGNM simultaneously estimates parameters from multiple initial estimates within the wide initial range (between the base values multiplied by 10−2 and 102). This method can estimate the distributions of the parameter values in the accepted parameter sets. The pharmacokinetic parameters of CysA were previously reported ones or optimized ones in this study by conventional parameter estimation as summarized in Table 1, and only those of CP‐I, and Ki,OATP1B and Ki,MRP2 were estimated by the CGNM.

The CGNM finds multiple sets of the best‐fit parameters in terms of the sum of squares residual (SSR). In theory, it is assumed that multiple sets of parameters exist with exactly the minimum SSR; however, due to the artifact of numerical computation, multiple sets of parameters with exactly the same SSR are not identified. Therefore, the parameter sets with SSRs that are considered very similar are determined in the following manner.

First, all parameter sets with SSRs that are significantly larger than the minimum SSR found by the CGNM are rejected. The statistical significance is determined by assuming the chi‐square distribution of SSR and a chosen alpha value of 0.05. Then, among the SSRs that are not rejected, the SSRs were plotted against the rank of the SSRs, and the rank at which the SSR increased sharply was identified using the knee of the rank–SSR curve (this method is also known as the elbow method; Figure S2A). The rank where the “elbow” is located was identified as follows; let rank Z be the largest rank that is accepted by the cut off defined by assuming chi‐square distribution of SSR. We find rank A that minimizes the following and defines it to be the location of the “elbow” and accept all points that are below or equal to rank A: (SSR at rank 1 + SSR at rank A)*(A − 1) + (SSR at rank A + SSR at Rank Z)*(Z‐A). This is equivalent to minimizing the area under the line segments (can be thought of as two trapezoids) connecting SSR at rank 1 and SSR at rank A, SSR at rank A, and SSR at rank Z in the plot presented in Figure S2A (the SSR vs. Rank plot). The elbow method is a commonly used heuristics in mathematical optimization when balancing the quantity and quality (in our case, quantified by SSR) of the accepted points.

The coefficient of variation (CV%) of the parameters was estimated by two methods using accepted parameter sets (observed data), and a bootstrap method, which was conducted as follows. First, bootstrap data sets were created by resampling the residuals with replacement. As we will use the accepted parameter sets as the initial estimates when re‐estimating for the bootstrap data sets, we create the same number of bootstrap data sets as the number of the accepted parameter sets. Second, the parameter sets are re‐estimated for the bootstrap datasets using the accepted parameter sets as the initial values for the parameter estimation. Third, CV%s are calculated by computing the means and SDs of the re‐estimated parameters. Through this approach, we calculate the CV% values including both residual variability as well as the uncertainty resulting from the choice of initial parameter sets for the optimization.

Software

The nonlinear least squares fitting software Napp version 2.31 on MacOS (version 10.14.6), was used in the optimization and simulation processes applied in this study. The weights for the nonlinear least‐squares calculations were the square root of the original values. The CGNM calculation was run on R implementation of CGNM software, version a.6.4.3 (http://www.bluetree.me/CGNmethod_for_PBPKmodels/, https://cran.r‐project.org/web/packages/CGNM/index.html).

RESULTS

Determination of the Ki, of CysA based on CP‐I kinetics

Kp,liver and fhCLint, were optimized to reproduce the blood concentration time profiles of CysA using the measured blood levels of CysA given orally at 20 mg and 75 mg (Figure 1a). The parameters are summarized in Table 1. Kp,liver and fhCLint were 1.52‐ and 1.85‐times greater than the previously reported values (16.7 and 0.587, respectively). Apart from FaFg (0.337 ± 0.017 at 20 mg, and 0.572 at 75 mg), the same parameter values (Table 1) were used to reproduce the blood concentration time profiles of CysA after oral administration.

FIGURE 1

Blood concentration time profiles of CysA and blood concentration time profiles of CP‐I under base line conditions, or after CysA administration. (a) CysA was given to participants 3 h or 1 h before probe drug administration in the clinical study. Time zero was set when baseline CP‐I was determined (−3.5 h before probe drug administration). Then, CysA was administered at 0.5 h (−3 h) for 75 mg (blue), and 20 mg or 75 mg CysA (green and red, respectively) was given at 2 h on the same time scale. The solid lines represent the lines calculated using the fitted parameters (summarized in Tables 1 and 2). (b) Pharmacokinetic parameters (CLint,all,unit) of CP‐I, and Ki,OATP1B were optimized by simultaneous optimization using the all data sets (20 and 75 mg [−1 h], and 75 mg [−3 h]). The solid lines represent the lines calculated using the fitted parameters. The optimized parameters are summarized in Tables 1 and 2

The CLint,all,unit of CP‐I and the Ki,OATP1B of CysA were optimized by fitting the PBPK model to account for the observed data under all dosing conditions (control, 20 mg, 75 mg [−1 h] and 75 mg [−3 h]) simultaneously (Figure 1b; Table 2). The CLint,all,unit of CP‐I was 1.09‐fold greater than the corresponding value in our previous analysis.

Prediction of the OATP1B‐mediated DDI using the optimized Ki, value for CP‐I

The Ki,OATP1B optimized for CP‐I in this study was used to simulate the blood concentration of the probe drugs (pitavastatin and rosuvastatin) in a DDI. The CLint,all,unit were optimized to account for the blood concentration time profiles of the probe drugs. For rosuvastatin, ka was also optimized to account for the blood concentrations in the control phase. The dose‐dependent effect of CysA could be sufficiently reproduced for both drugs (Figure 2a and b). The Cmax, AUC, Cmax ratio (CmaxR), and AUC ratio (AUCR; CysA administered phase/control phase) of pitavastatin and rosuvastatin were selected as metrics to indicate the precision of the predicted values. AUC and AUCR of pitavastatin were close to the observed values (Figure 2c). The Cmax and CmaxR of pitavastatin displayed a relatively large difference between the observed and estimated values (5.7–32% and 29–36%, respectively) using Ki,OATP1B optimized for CP‐I (Figure 2c; Table S1). In the case of rosuvastatin, the simulation underestimated all the metrics: Cmax (8.7–49%) CmaxR (28–44%), AUC (2.5–31%), and AUCR (22–30%). The optimized Ki,OATP1B values based on the pitavastatin and rosuvastatin data were 0.631 ± 0.106 nM and 0.448 ± 0.191 nM, respectively (see fitted lines in Figure S1A). This optimization of Ki,OATP1B did not improve the fit for pitavastatin and rosuvastatin data (Figure 2c, d; Table S1). The simulation was also conducted using the in vitro Ki,OATP1B1 of CysA (19 nM) obtained with a 1‐h pre‐incubation, which failed to reproduce the observed DDI impact (Figure S1B).

FIGURE 2

Simulation of dose‐dependent and dosing interval effects of CysA on the blood concentration time profiles of pitavastatin (a) and rosuvastatin (b). (a, b) Using Ki,OATP1B optimized for the CP‐I data, the dose‐dependent and dosing interval dependent effects of CysA were simulated using the PBPK models of pitavastatin a and rosuvastatin b given orally (0.2 and 1 mg, respectively) combined with the PBPK model of CysA. Pitavastatin or rosuvastatin was given 3 h or 1 h after CysA administration. Solid lines represented fitted lines, and broken lines represent the simulated lines. (c, d) Cmax, CmaxR, AUC, and AUCR of pitavastatin c and rosuvastatin d were calculated by connecting their PBPK models to the CysA PBPK model using the three different Ki,OATP1B; Ki,OATP1B (CP‐I), optimized for CP‐I data; Ki,OATP1B (PTV or RSV), optimized for PTV or RSV data; Ki,OATP1B (CP‐I, corrected), optimized for CP‐I corrected by the in vitro Ki,OATP1B1 ratio. Blood concentration time profiles calculated using the Ki,OATP1B (PTV or RSV) are shown in Figure S1A, and those calculated using the Ki,OATP1B (CP‐I, corrected) are shown in Figure S1C. Values of Cmax, CmaxR, AUC, and AUCR are summarized in Table S1 AUC, area under the plasma concentration time curve; AUCR, area under the plasma concentration time curve ratio; Cmax, maximum plasma concentration; CmaxR, maximum plasma concentration ratio; CP‐I, coproporphyrin I; PBPK, physiologically‐based pharmacokinetic; PTV, pitavastatin; RSV, rosuvastatin

We previously proposed a correction for substrate dependence in Ki,OATP1B using the in vitro Ki,OATP1B1 values, and in vivo Ki,OATP1B optimized for CP‐I in analyzing the effect of a single dose of rifampicin on the blood concentration profiles of CP‐I and probe OATP1B drugs. According to this approach, Ki,OATP1B,PTV (0.367 nM) was calculated from the ratio of the in vitro Ki,OATP1B1 ratio (0.684; Mochizuki et al., submitted) of pitavastatin to CP‐I, and this value was used in the DDI simulation (Figure S1C). However, the specific uptake of rosuvastatin by OATP1B1 in our in vitro assay was not high enough to determine the Ki,OATP1B1 for use in this analysis.

The time profiles of OATP1B inhibition were simulated at different CysA doses (beyond those administered to healthy subjects; Figure 3a). When CysA is administered at high clinical doses (600 mg) for tissue transplantation, 41 h are required to recover OATP1B function to 80%. We also simulated the impact of CysA dose and dosing interval on the AUCR and CmaxR for pitavastatin and rosuvastatin caused by OATP1B inhibition (Figure 3b). The effect of dosing interval was drug dependent. For example, the simulated DDI impact was highest for pitavastatin when CysA was administered simultaneously or 1 h before pitavastatin administration, whereas that for rosuvastatin was highest when CysA was given 1 or 2 h after rosuvastatin administration.

FIGURE 3

Time‐ and dose‐ dependency of OATP1B inhibition by CysA. (a) Time profiles of OATP1B inhibition by CysA were simulated at doses of 20, 75, 300, and 600 mg. The blood concentration time profiles of CysA at the doses of 20 mg and 75 mg were validated (Figure 1). Linearity in the pharmacokinetic parameters of CysA was assumed to calculate the blood concentration time profiles at 300 and 600 mg. (b) Simulation of CysA impact on the AUCR and CmaxR of pitavastatin and rosuvastatin. CysA is assumed to be administered at doses of 20, 75, 300, and 600 mg before 2 h to after 5 h of probe administration. AUCR, area under the plasma concentration time curve ratio; CmaxR, maximum plasma concentration ratio

Parameter search of CP‐I model using the Cluster Gauss–Newton method

In the previous analysis, most of the pharmacokinetic parameters for CP‐I were fixed, hence, the estimated value of Ki,OATP1B may depend on these fixed values. To confirm the independence of the optimized Ki,OATP1B from these fixed values, we used the CGNM. For the analysis using CGNM, pharmacokinetic parameters of CysA were cited from Table 1; however, the pharmacokinetic parameters for CP‐I were estimated from the wide range of initial iterates. The method successfully determined the distribution of the solutions of nine parameters (pharmacokinetic parameters for CP‐I, Ki,OATP1B, and Ki,MRP2; Figure 4) simultaneously, over a relatively wide range (102 – 104). Of 1000 initial parameter sets, 546 were selected according to the criteria using the chi‐squared test and elbow method (Figure S2A). These 546 parameter sets could accurately explain CP‐I blood concentrations under all dosing conditions (Figure 4a).

FIGURE 4

CGNM analysis of the CP‐I plasma concentration with or without CysA administration, and distribution of parameter values of the initial and corresponding optimized values. (a) Summary of the plasma concentration time profiles of CP‐I using the optimized parameters generated in 1000 cases. The green lines represent the profiles using the parameter sets accepted based on the chi‐square distribution of SSR and elbow method (Figure S2A). (b) Violin plots of the initial and optimized parameters in the selected 546 parameter sets. In each plot, a gray area indicates the distribution of the parameter values, a black dot in the center indicates the median, a vertical bar indicates interquartile range, solid lines stretched from the bar indicate the 25 percentile and 75 percentile values, and broken lines indicate the lower and upper adjacent values. The units of values are shown in Table S2. CGNM, Cluster Gauss–Newton method; CP‐I, coproporphyrin I; SSR, sum of squares residual

The comparison of initial and final distributions of the parameters were described in a violin plot (Figure 4b) and the median, maximum, minimum, and 25% and 75% percentiles were determined (Table S2). The distributions of the solutions of the parameters, including β, FaFg, the intrinsic biliary clearance per the sum of the intrinsic biliary and metabolic clearance (fbile), Ki of unbound CysA for MRP2 (Ki,MRP2), the ratio of CP‐I synthesis in the liver (fsyn), and the passive hepatic uptake clearance divided by the active hepatic uptake clearance (Rdif), did not change much from the initial distribution (Figure 4b). By contrast, distributions of the solutions of the overall intrinsic clearance (CLint, all) and synthesis rate (v syn) of CP‐I, and Ki,OATP1B were narrowed down considerably compared to the initial distributions. CLint, all is a hybrid parameter representing the rate of overall hepatic clearance, which is composed of the clearances for the sinusoidal influx and efflux, and canalicular efflux, according to the extended clearance concept. , Using β and Rdif, we calculated the values of these parameters for every accepted parameter set (Figure S2B). The ranges narrowed slightly but not as much as that observed for CLint,all (Figure 4b). Based on the CV% of each parameter estimated by a bootstrap method (Table S2), the CLint,all and Ki,OATP1B were considered identifiable, because their CV was 12.4% and 11.5%, respectively, even considering residual variability. The distribution of parameter values after bootstrapping is shown in Figure S2C. In addition, the medians of CLint, all and Ki,OATP1B (Table S2) were similar to the corresponding parameters calculated using the Napp program (Tables 1 and 2). This confirms that estimated Ki,OATP1B is independent of the values that were fixed during the analysis using Napp program. The median value of v syn (12.1 nmol/h = 0.188 nmol/h/kg) was similar to that (0.212–0.442 nmol/h/kg) reported for rifampicin DDI study.

DISCUSSION

The present study conducted a PBPK model‐based analysis using a conventional method and CGNM to support the feasibility of the translation of clinical CP‐I data to DDI predictions using OATP1B probe drugs, pitavastatin and rosuvastatin.

The PBPK model analysis of CP‐I and CysA data (20 mg and 75 mg with a 1‐h interval) yielded an estimate of Ki,OATP1B of CysA, which could reasonably describe the dose‐dependent effect of CysA on the blood concentration time profiles of CP‐I (Figure 1). Notably, using the Ki,OATP1B obtained using CP‐I data, simulation of the DDI with OATP1B probes could reproduce the observed dose‐ and dosing interval dependent effect of CysA on the CmaxR and AUCR of both pitavastatin, although the simulation somewhat underestimated the DDI impact (percentage of the difference from the observed values, 29–32% and 1.6–14%, respectively) and rosuvastatin (28–44% and 22–30%, respectively; Figure 2; Table S1). Ki,OATP1B optimization to account for the pitavastatin and rosuvastatin data yielded the values similar to the Ki,OATP1B optimized for the CP‐I data (18% and 16% difference for pitavastatin and rosuvastatin, respectively), but it could improve the prediction of both CmaxR and AUCR for rosuvastatin (Figure 2c and d). Previously, we proposed correction of the substrate dependence in the DDI prediction by taking the ratio of the in vitro Ki,OATP1B1 of rifampicin for CP‐I and pitavastatin. In the case of CysA, no obvious substrate‐dependence of the in vitro Ki values to be corrected for were observed between pitavastatin and CP‐I (CP‐I and pitavastatin; 19 and 13 nM, respectively), and thus, the correction produced only a negligible impact on the DDI simulation (Figure S1C). On the other hand, the specific uptake of rosuvastatin in our OATP1B1 expression system was not high enough to determine Ki,OATP1B1 for CysA to support the substrate dependence in Ki,OATP1B between rosuvastatin and CP‐I. We note that, even though CysA has been reported to inhibit BCRP in vitro, in the current analysis, we did not consider inhibition of BCRP‐mediated transport of rosuvastatin by CysA in the PBPK model. This limitation may have resulted in underestimation of the k a and/or FaFg of rosuvastatin when co‐administered with CysA. This remains a topic for future work to describe the gastrointestinal kinetics of CysA precisely and its effects on intestinal BCRP in the PBPK model.

Although the optimized Ki,OATP1B in this study was 3‐to‐4 fold smaller than the previously reported values for pitavastatin and fluvastatin (2.16 and 1.76 nM, respectively), considering the difference in fp used in the analysis (0.11 in the previous study vs. 0.024 in this study), it should be noted that pharmacokinetic data for CysA and probe drugs (pitavastatin and rosuvastatin) in the previous PBPK model‐based analysis were derived from separate studies, and thus, there is a possibility that inter‐experimental differences modified Ki,OATP1B. The successful translation of CP‐I data to the probe drugs strongly supports the feasibility of examining the effect of new investigational drugs on the CP‐I concentration as a promising approach to predict the DDI impact on the AUC of these probe drugs, by overcoming the discrepancy between in vivo Ki,OATP1B and in vitro Ki,OATP1B1.

The parameter optimization was conducted additionally using the in vitro Ki,OATP1B as an initial value to elucidate the impact of the limited data points when estimating Ki,OATP1B. First, it was confirmed that the optimized Ki,OATP1B values of CysA for CP‐I using 20 mg or 75 mg data alone were consistent with each other, and that these values were consistent with those obtained using all the data (Figure S3). Hence, the collection of data on dose‐dependent effects is not essential to estimate reliable Ki,OATP1B values in clinical studies. Second, the impact of removing data points on parameter estimates was explored. Ki,OATP1B could be estimated reasonably, but more data points were required to optimize CLint,all,unit than Ki,OATP1B because the estimated values had large deviations or reached the limit of the range that was set to avoid divergence under most conditions in which we removed data points (Figure S3). These analyses supported the idea that CP‐I is a robust endogenous OATP1B biomarker for obtaining an estimate of Ki,OATP1B.

Parameter optimization by conventional analysis requires expert knowledge. Previously, we used the CGNM to investigate the PBPK‐model based analysis of P‐450‐mediated DDI. This approach circumvents investigators from the laborious initial parameter estimation, and selection of identifiable parameters to avoid divergence. Recently, CGNM was developed for comprehensive parameter set survey to yield multiple approximate solutions. The CGNM can easily elucidate the characteristics of the parameters, uniquely identifiable or not, and any cross‐correlation between the parameters. Conventional algorithms require some parameter values to be fixed. Hence, sensitivity analysis is routinely conducted as a post hoc analysis to examine the effect of fixed values on the optimized parameter values. An analysis of the CGNM conducted in this study to explore such effects revealed that based on CV% (Table S2), the CLint, all of CP‐I and Ki,OATP1B of CysA were uniquely identifiable using the observed data set. CLint, all is a complex parameter composed of the clearances for the sinusoidal influx and efflux, and canalicular efflux; however, the observed data sets (blood concentration time profiles alone) were insufficient to identify every intrinsic clearance (Figure S2B). Sampling of a wide set of initial values in the analysis was supported by the fact that the accepted parameter sets included extreme cases where both β and fsyn were always close to one (Figure S2D). In this extreme case, the input rate of CP‐I was determined solely by biliary excretion followed by absorption from the intestine to the blood circulation; eventually, the v syn had to be greater (Figure 4) because of the low FaFg of CP‐I to account for the observed blood concentrations. Nonetheless, the Ki,OATP1B values in the parameter sets, even in such an extreme case, were within the deviation limits of this parameter.

Among the endogenous OATP1B biomarkers, CP‐I is frequently used to demonstrate the OATP1B‐mediated DDI potential of approved drugs, , , , , and new investigational drugs. , , A highlight of this study was the successful translation of CP‐I data using the PBPK model to overcome the discrepancy between in vivo Ki,OATP1B and in vitro Ki,OATP1B1, and to predict more precisely the OATP1B‐mediated DDI with two OATP1B probe substrates, pitavastatin and rosuvastatin. Based on this successful translation, a workflow for DDI risk assessment using the endogenous biomarkers was developed (Figure 5). CmaxR is considered a better metric than AUCR because the AUCR value depends on the time interval to calculate AUC. Because CmaxR can practically approximate 1 + (CysA concentration at time to Cmax [Tmax] of CP‐I)/Ki,OATP1B, assuming fOATP1B = 1 (Supplemental Text), we can obtain an estimate of Ki,OATP1B easily. Indeed, Ki,OATP1B estimates yielded by this approximation were reasonable considering the optimized value (0.536 nM in this study; Table S3; Figure S4). Similarly, the Ki,OATP1B of rifampicin estimated from CmaxR in our previous studies were 0.119 and 0.126 μM (Table S3A; Figure S4), which were also similar to the optimized values (0.0824–0.106 μM ). It should be noted that AUCR of CP‐I, and mean concentrations of CysA or rifampicin (AUC/calculation interval) could yield Ki,OATP1B values similar to the corresponding estimates using CmaxR (Table S3B; Figure S4).

FIGURE 5

Scheme of the workflow for predicting OATP1B‐mediated DDIs using an endogenous OATP1B biomarker in drug development. DDIs, drug–drug interaction

Ki,OATP1B optimized for CP‐I could yield reasonable prediction values at different doses (20 and 75 mg), and dosing intervals (1 and 3 h; Figure 2; Table S1). The accumulation of such successful translation data could be expected to lead to a waiver for the clinical DDI studies using OATP1B probe drugs in the future. A reliable estimate of in vivo Ki,OATP1B enables the extrapolation of CysA doses to those used for tissue transplantation in DDI predictions, as confirmation using OATP1B probe drugs is not practical in healthy subjects because of the severity of the adverse reactions (Figure 3a), and also the extrapolation of CysA dosing intervals. Generally, the OATP1B perpetrators were given simultaneously or before administration of the probe drug. Our simulation suggests the maximum DDI impact will be obtained when CysA is given after rosuvastatin administration, while simultaneously or before pitavastatin administration as far as we considered OATP1B inhibition (Figure 3b). PBPK model analysis has the potential to suggest better design to assess the magnitude of DDIs depending on the probe drug in drug development, and to avoid clinically relevant DDIs. Furthermore, PBPK model analysis is anticipated to discriminate OATP1B‐mediated DDI from overall DDIs caused by multiple mechanisms. For instance, the DDI of rosuvastatin includes breast cancer resistance protein modulation in addition to OATP1B inhibition.

In conclusion, PBPK model‐based analysis using conventional parameter estimation methods and CGNM yielded a Ki,OATP1B of CysA using CP‐I data, which could aid in the translation of CP‐I data to DDI predictions with OATP1B probe drugs.

CONFLICT OF INTEREST

Tadayuki Takashima is an employee of Asahi Kasei Pharma; Kenta Yoshida and Jialin Mao are employees of Genentech; Yurong Lai is an employee of Gilead Sciences; Kunal Taskar and Maciej J. Zamek‐Gliszczynski are employees of GlaxoSmithKline; Kevin Rockich is an employee of Incyte Research Institute; Xiaoyan Chu is an employee of Merck & Co., Inc.; Yoshiyuki Yamaura is an employee of Ono Pharmaceuticals; and Hideki Hirabayashi is an employee of Takeda. All other authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

T.M., K.Y., Y.L., H.H., Y.Y., K.R., K.T., T.T., X.C., Y.S., and H.K. wrote the manuscript. T.M., K.Y., Y.L., H.H., Y.Y., K.R., K.T., T.T., X.C., M.J.Z., J.M., K.F., Y.S., and H.K. designed the research. T.M., T.Y., Y.S., and H.K. performed the research. T.M., Y.A., K.M., Y.S., and H.K. analyzed the data. Y.A., T.Y., and Y.S. contributed new reagents/analytical tools.

Supinfo

Click here for additional data file.