Physiologically Based Pharmacokinetic Models of Probenecid and Furosemide to Predict Transporter Mediated Drug-Drug Interactions

Purpose To provide whole-body physiologically based pharmacokinetic (PBPK) models of the potent clinical organic anion transporter (OAT) inhibitor probenecid and the clinical OAT victim drug furosemide for their application in transporter-based drug-drug interaction (DDI) modeling. Methods PBPK models of probenecid and furosemide were developed in PK-Sim®. Drug-dependent parameters and plasma concentration-time profiles following intravenous and oral probenecid and furosemide administration were gathered from literature and used for model development. For model evaluation, plasma concentration-time profiles, areas under the plasma concentration–time curve (AUC) and peak plasma concentrations (Cmax) were predicted and compared to observed data. In addition, the models were applied to predict the outcome of clinical DDI studies. Results The developed models accurately describe the reported plasma concentrations of 27 clinical probenecid studies and of 42 studies using furosemide. Furthermore, application of these models to predict the probenecid-furosemide and probenecid-rifampicin DDIs demonstrates their good performance, with 6/7 of the predicted DDI AUC ratios and 4/5 of the predicted DDI Cmax ratios within 1.25-fold of the observed values, and all predicted DDI AUC and Cmax ratios within 2.0-fold. Conclusions Whole-body PBPK models of probenecid and furosemide were built and evaluated, providing useful tools to support the investigation of transporter mediated DDIs. Supplementary Information The online version contains supplementary material available at 10.1007/s11095-020-02964-z.


INTRODUCTION
Many important drug transporters are members of the solute carrier (SLC) family, which is widely expressed throughout the human body and mediates influx or efflux of endogenous and exogenous substrates (1). For the approval of new drugs, the U.S. Food and Drug Administration (FDA) and the European Medicines Agency (EMA) require various in vitro, and in many cases, clinical studies to characterize the transporter mediated drug-drug interaction (DDI) potential of investigational drugs. Based on the outcome of these investigations, recommendations for dose adjustments are given in the accompanying label of a new drug (2,3). Organic anion transporter (OAT) 1 and OAT3 are members of the SLC transporter family (SLC22A6, SLC22A8) and recognized as important drug transporters from the perspective of their potential to be involved in clinically relevant DDIs. OAT1 and OAT3 are predominantly expressed in the kidney at the basolateral membrane of proximal tubule cells (4), where they facilitate the uptake of endogenous (e.g. paminohippurate, estrone sulfate) and exogenous (e.g. diuretic drugs) organic anions from the blood into the proximal tubule cells (5), from where they can be secreted into the nephron lumen for excretion with the urine. Several polymorphisms have been identified and investigated; however, variants of OAT1 or OAT3 have not been shown to significantly impact the renal clearance of OAT substrates in clinical studies (6,7). To characterize the OAT mediated DDI potential, the FDA recommends furosemide as clinical OAT1/OAT3 substrate and probenecid as clinical OAT1/OAT3 inhibitor (8). In addition, probenecid can also be used to investigate organic anion transporting polypeptide (OATP) 1B1 mediated DDIs (9). OATP1B1, another clinically relevant member of the SLC transporter family (SLCO1B1), is exclusively expressed at the sinusoidal membrane of hepatocytes, where it is responsible for the uptake of endogenous (e.g. bile acids) and exogenous (e.g. statins, rifampicin) organic anions from the blood into the hepatocytes (10,11). As probenecid also inhibits OATP1B1, the prediction of the probenecid-rifampicin DDI was included into the presented study, to further evaluate the performance of the probenecid model as a DDI perpetrator.
Physiologically based pharmacokinetic (PBPK) modeling is encouraged and supported by the FDA and EMA as a valuable tool to quantitatively describe and predict the pharmacokinetics of drugs, to evaluate DDI potential and to support clinical study design, dose selection and labeling during drug development (2,3,(12)(13)(14). The objectives of this study were to provide whole-body PBPK models of probenecid and furosemide, incorporating the transporters and enzymes involved in the pharmacokinetics of these drugs. In vitro measurements were used to parametrize the respective incorporated processes. The models were built and evaluated to adequately predict the plasma concentration-time profiles and the fractions excreted unchanged in urine. Furthermore, the models were used to predict probenecid DDIs, with probenecid as potent clinical OAT1/OAT3 inhibitor (8,15), moderate OATP1B1 inhibitor (9), weak inhibitor of multidrug resistance-associated protein (MRP) 4 (16) and weak inhibitor of uridine 5′-diphospho-glucuronosyltransferase (UGT) 1A9 (in-house measurement), furosemide as clinical OAT1/OAT3 substrate and rifampicin as OATP1B1 substrate. The comprehensive Electronic Supplementary Material (ESM) to this manuscript provides detailed information on the developed PBPK models, including all model parameters and a complete documentation of the extensive model evaluation. The model files will be shared in the Open Systems Pharmacology PBPK model library (www.open-systems-pharmacology.org).

Software
PBPK modeling was performed with the open source PK-Sim® and MoBi® modeling software (version 8.0, part of the Open Systems Pharmacology Suite (17,18), www.opensystems-pharmacology.org). Published plasma concentrationtime profiles were digitized using GetData Graph Digitizer (version 2.26.0.20, S. Fedorov) (19). Parameter optimizations were accomplished with the Monte Carlo algorithm as well as the Levenberg-Marquardt algorithm using the "multiple random starting values" function implemented in PK-Sim®. The final optimizations were run using the Levenberg-Marquardt algorithm. PK parameter analysis and calculation of model performance measures was performed with R (version 3.6.2, The R Foundation for Statistical Computing) and graphics were compiled with R and RStudio (version 1.2.5033, RStudio, Inc., Boston, MA, USA). Sensitivity analysis was performed using the implemented Sensitivity Analysis tool in PK-Sim® (20).

PBPK Model Building
PBPK model building was started with an extensive literature search to collect physicochemical parameters, information on absorption, distribution, metabolism and excretion (ADME) processes and clinical studies of intravenous and oral administration of probenecid and furosemide in single-and multiple dose regimens.
To build the datasets for PBPK model development, the reported observed plasma concentration-time profiles were digitized and divided into a training dataset for model building and a test dataset for model evaluation. Model input parameters that could not be informed from experimental reports were optimized by fitting the model simultaneously to the observed data of all studies assigned to the training dataset. To limit the parameters to be optimized during model building, the minimal number of processes necessary to mechanistically describe the pharmacokinetics and DDIs of the modeled drugs were implemented into the models. If two transporters show very similar expression patterns and affinity for the same compound, optimizing the transport rate constants of both transporters would lead to identifiability issues. Therefore, only the transporter with the higher affinity for the respective substrate was implemented, to describe a transport that probably is accomplished by both transporters in vivo.
Model evaluation was carried out based on the clinical data of the test dataset. Descriptive (training dataset) and predictive (test dataset) performance of the model for all analyzed clinical studies is transparently documented in the ESM.

Virtual Individuals
The PBPK models were built based on data from healthy individuals, using the reported sex, ethnicity and mean values for age, weight and height from each study protocol. If no demographic information was provided, the following default values were substituted: male, European, 30 years of age, 73 kg body weight and 176 cm body height (characteristics from the PK-Sim® population database (21,22)). ADME transporters and enzymes were implemented in accordance with literature, using the PK-Sim® expression database to define their relative expression in the different organs of the body (23). Table S7.0.1 summarizes all system-dependent parameters on the implemented transporters and enzymes.

Virtual Population Characteristics
To predict the variability of the simulated plasma concentration-time profiles, virtual populations of 100 individuals were generated, consisting of either European or Asian individuals. The percentage of female individuals and the ranges of age, weight and height were set according to the reported demographics. If not specified, virtual populations containing 100 male subjects 20-50 years of age were used, with body weight and height restrictions from the PK-Sim® population database (22). For details on the study populations see Tables S2.2 In the generated virtual populations, organ volumes, tissue compositions, blood flow rates, etc. were varied by an implemented algorithm within the limits of the International Commission on Radiological Protection (ICRP) (21,22) or Tanaka (24) databases. In addition, the reference concentrations of the implemented transporters and enzymes were lognormally distributed around their mean values, using reported variabilities for their expression from the PK-Sim® database (25) or from literature. Table S7.0.1 summarizes the modeled transporters and enzymes with their reference concentrations and variabilities.
As the clinical plasma concentration data from literature is mostly reported as arithmetic means ± standard deviation, population prediction arithmetic means and 68% prediction intervals were plotted, that correspond to a range of ±1 standard deviation around the mean assuming normal distribution.

PBPK Model Evaluation
Model performance was evaluated using various methods. The population predicted plasma concentration-time profiles were compared to the data observed in the clinical studies. Furthermore, predicted plasma concentrations of all studies were compared to the observed plasma concentrations in goodness-of-fit plots. In addition, the model performance was evaluated by comparison of predicted to observed areas under the plasma concentration-time curve (AUC) from the time of drug administration to the time of the last concentration measurement (AUC last ) and peak plasma concentration (C max ) values. As quantitative measures of the model performance, the mean relative deviation (MRD) of all predicted plasma concentrations (Eq. 1) and the geometric mean fold error (GMFE) of all predicted AUC last and C max values (Eq. 2) were calculated. MRD and GMFE values ≤2 characterize an adequate model performance.
where c predicted, i = predicted plasma concentration, c observed, i = corresponding observed plasma concentration and k = number of observed values.
where predicted PK parmeter i = predicted AUC last or C max , observed PK parameter i = corresponding observed AUC last or C max and m = number of studies. Furthermore, sensitivity analysis results were assessed. A detailed description of the sensitivity calculation is provided in the ESM.

PBPK DDI Modeling
As an additional means of model evaluation, the DDI performance of the developed models was assessed. To model the probenecid-furosemide DDI, inhibition of OAT3, UGT1A9 and MRP4 by probenecid was implemented. To predict the probenecid-rifampicin DDI, inhibition of OATP1B1 by probenecid was incorporated (Fig. 1). The rifampicin model applied was developed by Hanke et al. (26) and is freely available in the Open Systems Pharmacology repository on GitHub (27). The model includes rifampicin transport via OATP1B1 and P-glycoprotein (Pgp), metabolism via the arylacetamide deacetylase (AADAC), as well as auto-induction of OATP1B1, Pgp and AADAC (26). The good DDI performance of the model was demonstrated in many different applications (26,(28)(29)(30)(31). Mathematical implementation of the DDI processes is specified in Section 1 of the ESM.
Inhibition constants characterizing the inhibition of OAT3, UGT1A9 (in-house measurement) and OATP1B1 by probenecid were taken from in vitro experimental reports (9,15). Short descriptions of the respective in vitro assay conditions are provided in Sections 5.1 and 6.1 of the ESM. To describe the competitive inhibition of MRP4 by probenecid, the corresponding inhibition constant was optimized during the furosemide PBPK model parameter identification, using the clinical data of one of the probenecid-furosemide interaction studies (32) (see Table S5.2.1). The DDI parameter values are listed in the probenecid drug-dependent parameter table (Table S2.3.1).

PBPK DDI Performance Evaluation
All DDI predictions were evaluated by comparison of predicted versus observed victim drug plasma concentration-time profiles alone and during co-administration, DDI AUC last ratios and DDI C max ratios (Eq. 3). DDI PK parameter ratio ¼ PK parameter victim drug during co−administration PK parameter victim drug alone ð3Þ where PK parameter = AUC last or C max . As quantitative measure of the DDI prediction accuracy, GMFEs of the predicted DDI AUC last ratios and DDI C max ratios were calculated according to Eq. 2.

PBPK Model Building and Evaluation
The probenecid PBPK model was developed using 27 different clinical studies, including intravenous (single dose) and oral (single-and multiple dose) administration. A complete list of the clinical studies used in the presented analysis is provided in Table S2.2.1. In addition, five studies reported fraction excreted unchanged in urine profiles following oral administration. In the intravenous studies, probenecid was administered in doses of 464-1860 mg. In the oral studies, probenecid was administered in doses of 250-2000 mg. The training dataset included 11 probenecid plasma concentration-time profiles and one fraction excreted unchanged in urine profile. The final probenecid model applies uptake into kidney cells via OAT3, glucuronidation mainly in the renal cells by UGT1A9, glomerular filtration and tubular reabsorption,   Table S2.2.1. iv: intravenous, po: oral, sd: single dose, tab: tablet. a b which was modeled as a reduction of the glomerular filtration rate (GFR fraction <1). The drug-dependent parameters are summarized in Table S2.3.1.
Population predicted compared to observed plasma concentration-time profiles of selected studies as well as the probenecid goodness-of-fit plot with the predicted plasma concentrations of all studies are presented in Fig. 2. Semilogarithmic and linear plots of the plasma profiles of all 27 clinical studies included in the analysis are shown in the ESM. Population predicted compared to observed fraction excreted unchanged in urine profiles are also presented in the ESM. Table S2   The correlation of predicted and observed furosemide AUC last and C max values is presented in Fig. 5, further demonstrating the good model performance with 41/42 predicted AUC last and 24/25 predicted C max values within 2-fold of the observed data. The individual values, mean GMFE values and ranges are listed in Table S3.5.2.
The sensitivity analysis results of a simulation of 80 mg furosemide once daily as a tablet are illustrated in Fig.  S3.5.3. Applying a threshold of 0.5, the furosemide model is sensitive to the values of furosemide fraction unbound in plasma (literature value) and the OAT3 catalytic rate constant (optimized).

PBPK DDI Modeling and Evaluation
The developed PBPK models were applied to model the probenecid-furosemide and probenecid-rifampicin DDIs and the DDI performance was evaluated using the clinical data of six studies investigating the probenecid-furosemide DDI and one study of the probenecid-rifampicin DDI. For all studies, plasma concentration-time profiles of the victim drugs, administered alone and during probenecid co-administration, were predicted and compared to observed data. In addition, four studies of the probenecid-furosemide DDI reported fraction excreted unchanged in urine profiles, allowing the comparison of predicted and observed urinary excretion under control and DDI conditions. Administration protocols, study population details and references of the clinical DDI studies are listed in Tables S5.2.1 and S6.2.1.
To predict the probenecid-furosemide DDI, competitive inhibition of OAT3 (K i = 5.41 μmol/l) (15) and noncompetitive inhibition of UGT1A9 (K i = 242.0 μmol/l) (inhouse measurement) by probenecid were implemented using interaction parameter values measured in vitro. As no information regarding the inhibition of MRP4 could be obtained, the K i to describe the competitive inhibition of MRP4 (K i = 87.4 μmol/l) by probenecid was optimized during the parameter identification of the furosemide model. To predict the probenecid-rifampicin DDI, competitive inhibition of OATP1B1 (K i = 39.8 μmol/l) (9) was implemented using an interaction parameter value measured in vitro.
The coupled models adequately describe and predict all furosemide and rifampicin plasma concentration-time profiles, as well as fraction excreted unchanged in urine profiles of furosemide, under control conditions and during probenecid co-administration, over the full range of reported DDI administration protocols. Population predicted compared to observed plasma profiles of furosemide and rifampicin, administered alone and during the DDIs, are presented in The correlation of predicted and observed DDI AUC last ratios and DDI C max ratios of all analyzed clinical DDI studies is shown in Fig. 7, further demonstrating the good DDI performance with all predicted DDI ratios within 2-fold of the observed data. The individual ratios, mean GMFE values and ranges for both DDI combinations are listed in Table I.

DISCUSSION
The newly developed whole-body PBPK models for probenecid and furosemide accurately describe the observed plasma concentration-time profiles and fraction excreted unchanged in urine data over the full range of reported doses and administration protocols. Furthermore, these models adequately describe the available clinical data from probenecid-furosemide and probenecid-rifampicin clinical DDI studies.
Development of the probenecid model was particularly challenging. The number of published clinical studies is low and the quality of the available data varies considerably, as probenecid was approved in the late 1940s. Hence, careful consideration of the study protocols and the presented data was required, and studies in patients or elderly volunteers, studies with co-medication and studies using others than the marketed formulation were excluded for probenecid model development. In addition, information on the ADME processes that govern the pharmacokinetics of probenecid is very limited.
Probenecid shows a low solubility and permeability, indicating an important role of transporters in its absorption and distribution. However, neither in vitro nor in vivo studies describing transporters involved in probenecid absorption, organ uptake, secretion or reabsorption are available in literature. Therefore, probenecid absorption was described by optimization of the passive transcellular intestinal permeability (optimized value: 3.97 · 10 −4 cm/min, lipophilicity based calculated value: 3.12 · 10 −6 cm/min). The uptake of probenecid into the kidney, which is its site of action and metabolism (39)(40)(41), was assumed to be mediated by OAT3. The parameters (K M , k cat ) to describe this transport were optimized during parameter identification (see Table S2.3.1). Furthermore, a low renal clearance (42) and fraction excreted unchanged in urine of only 0.3% to 5% (43,44) are reported, indicating tubular reabsorption. Due to our current lack of knowledge regarding transporters that may contribute to probenecid reabsorption, the GFR fraction was optimized to 0.03. This reduced GFR fraction substitutes for the implementation of active reabsorption processes of probenecid (20) and correctly captures the low probenecid fraction excreted unchanged in urine.
In the clinical studies conducted by Vree et al. (43,44), probenecid tablets were broken in half prior to oral administration. The corresponding plasma concentration-time profiles display an earlier time to peak plasma concentration (T max ) of 1.6 h compared to the other clinical studies with a T max of 3.3 h. Given the low solubility of probenecid it is possible that the broken tablets show a different dissolution behavior, resulting in faster release and absorption. Therefore, a different dissolution profile was used to describe the studies by Vree et al. (43,44). The parameters to model the two different dissolution profiles are listed in Table S2.3.1.
Similar to probenecid, furosemide also demonstrates low solubility and permeability and is classified as a BCS class IV drug (45). Therefore, transporters play an essential role in furosemide absorption, distribution and elimination. Furosemide bioavailability is highly variable (37%-83%) (46), and influenced by dosage form and fasted/fed state of the patient (46). For the absorption of furosemide, Flanagan et al. postulate a saturable active transport process and passive diffusion with paracellular contribution in Caco-2 cells (47). As no further information on transporters that may contribute to the absorption of furosemide is available, absorption was modeled as passive transcellular (5.06 · 10 −7 cm/min) and paracellular (2.32 · 10 −6 cm/min) intestinal permeability. These processes together allow a rapid absorption in the small intestine to describe the early furosemide T max of 1.0-1.5 h (46), while limiting the furosemide absorption in the large intestine. Other crucial transport processes take place in the kidney, which is the main organ for furosemide metabolism and excretion (fraction excreted unchanged in urine: 20%-80% (46)).
Uptake of furosemide into renal cells in vivo is probably facilitated via OAT1 and OAT3. Both transporters are predominantly expressed in the kidney (48) and show a similar affinity to furosemide (OAT1 furosemide K M = 38.9 μmol/l, OAT3 furosemide K M = 21.5 μmol/l (49)). Without additional information to distinguish the furosemide transport via these two transporters, the furosemide OAT1 and OAT3 transport rate constants would be highly correlated in a model parameter optimization. To avoid indentifiability issues, renal uptake of furosemide was incorporated via OAT3 only, as a substitute for transport by both, OAT1 and OAT3. The probenecid inhibition potency towards these transporters is also similar (probenecid OAT1 K i = 11.4 μmol/l, probenecid OAT3 K i = 5.41 μmol/l (15) and the DDI was predicted via inhibition of OAT3. Furosemide metabolism was modeled using UGT1A9 (50). The pronounced urinary excretion is accomplished by glomerular filtration and active tubular secretion. Renal secretion was incorporated via MRP4 (51,52), based on in vitro furosemide transport measurements.
In the goodness-of-fit plot, there seems to be an underprediction of the lower furosemide plasma concentrations for two of the intravenous studies (Rupp 1974 andWaller 1982). These are the only two studies that published furosemide plasma concentrations later than 6 h after intravenous dosing, and they did not report their lower limits of quantification. The model is therefore not qualified for the prediction of plasma concentrations later than 6 h after intravenous administration of furosemide. However, this does not affect the prediction of higher plasma concentrations and AUC last values, and the developed model shows a good performance, with 41/42 predicted AUC last and 24/25 predicted C max values within 2-fold of the observed data.
To model the probenecid-furosemide DDI, inhibition of OAT3 alone, using published in vitro inhibition parameters, was not sufficient to describe the clinically observed data. Therefore, inhibition of UGT1A9 and MRP4 were added to adequately capture the impact of probenecid on the furosemide pharmacokinetics. Inhibition parameter values for OAT3 (15) and UGT1A9 (in-house measurement) were available from in vitro studies. Van Aubel et al. reported MRP4 inhibition by probenecid (16), but so far, no inhibition parameter values have been published. Therefore, competitive MRP4 inhibition was assumed and the K i value was optimized during parameter identification of the furosemide model. Applying the inhibition of OAT3, UGT1A9 and MRP4, all reported plasma concentration time-profiles and fraction excreted unchanged in urine profiles of furosemide during probenecid co-administration are well described by the presented models. The probenecid-rifampicin DDI was predicted using a probenecid OATP1B1 inhibition value measured with 2′,7′dichlorofluorescein as the substrate (9). For OATP1B1, it has been demonstrated that its inhibition can strongly depend on the employed substrate (53). However, since there are no in vitro reports of probenecid OATP1B1 inhibition using rifampicin, the DDI was predicted applying the probenecid OATP1B1 K i = 39.8 μmol/l measured with 2′,7′-dichlorofluorescein (9). Rifampicin is typically not used as OATP1B1 victim drug, but rather as OATP1B1/OATP1B3 inhibitor. However, as the clinical data of this probenecid-rifampicin trial were published, we wanted to utilize them to test our probenecid model. In addition to OATP1B1/OATP1B3 inhibition, rifampicin is also inhibiting and inducing further enzymes and transporters (13). Therefore, an impact of rifampicin on the perpetrator drug probenecid may have influenced the results of this particular DDI administration protocol. Taking into account that the only available clinical study has been published in 1973 and that the reported rifampicin plasma concentrations show considerable standard deviations, this DDI is also well described (see Fig. 6 and Table I).

CONCLUSIONS
The presented whole-body PBPK models of probenecid and furosemide have been carefully built and evaluated for their ability to predict the pharmacokinetics of these drugs, using a multitude of clinical studies. In addition, the models adequately describe the available clinical data of the probenecid-furosemide and probenecid-rifampicin DDIs and will be shared in the Open Systems Pharmacology PBPK model library (www.open-systemspharmacology.org) as tools to support the investigation of the DDI potential of new compounds during drug development. The ESM to this paper has been compiled to serve as transparent and comprehensive documentation of the probenecid and furosemide model development and evaluation.

ACKNOWLEDGMENTS AND DISCLOSURES
We would like to thank Solvo Biotechnology (Budapest, Hungary) for providing MRP4 vesicles used in the transporter studies. In vitro work on MRP4 mediated furosemide transport kinetics was conducted by Revathi Chapa and  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.