Delineating the Role of Various Factors in Renal Disposition of Digoxin through Application of Physiologically Based Kidney Model to Renal Impairment Populations

Development of submodels of organs within physiologically-based pharmacokinetic (PBPK) principles and beyond simple perfusion limitations may be challenging because of underdeveloped in vitro-in vivo extrapolation approaches or lack of suitable clinical data for model refinement. However, advantage of such models in predicting clinical observations in divergent patient groups is now commonly acknowledged. Mechanistic understanding of altered renal secretion in renal impairment is one area that may benefit from such models, despite knowledge gaps in renal pathophysiology. In the current study, a PBPK kidney model was developed for digoxin, accounting for the roles of organic anion transporting peptide 4C1 (OATP4C1) and P-glycoprotein (P-gp) in its tubular secretion, with the aim to investigate the impact of age and renal impairment (moderate to severe) on renal drug disposition. Initial PBPK simulations based on changes in glomerular filtration rate (GFR) underestimated the observed reduction in digoxin renal excretion clearance (CLR) in subjects with moderately impaired renal function relative to healthy. Reduction in either proximal tubule cell number or the OATP4C1 abundance in the mechanistic kidney model successfully predicted 59% decrease in digoxin CLR, in particular when these changes were proportional to reduction in GFR. In contrast, predicted proximal tubule concentration of digoxin was only sensitive to changes in the transporter expression/ million proximal tubule cells. Based on the mechanistic modeling, reduced proximal tubule cellularity and OATP4C1 abundance, and inhibition of OATP4C1-mediated transport, are proposed as possible causes of reduced digoxin renal secretion in renally impaired patients.


Introduction
Dosage adjustment is often required in patients with impaired renal function because of the impact this condition may have in altered drug clearance through either decreased renal excretion alone or in combination with reduced metabolism. Ideally, such decisions on dosage adjustment will be supported by recommendations on drug labels that are currently based on statistical analyses of clinical data, often from dedicated but small sized studies (US Food Drug Admin, 2010;Matzke et al., 2011;European Medicines Agency, 2014). The pharmaco-statistical nature of the analysis limits the ability to extrapolate findings beyond the boundaries of the patient groups originally studied. A direct consequence of this, in conjunction with paucity of data in severe renal impairment, has been the large proportion of drug labels void of any recommendations in the most vulnerable renal impairment patients according to recent survey of labels for drugs approved by the Food and Drug Administration in 2013 and 2014 (Jadhav et al., 2015).
A possible solution is the use of physiologically based pharmacokinetic (PBPK) models with embedded organ models that incorporate mechanisms of local disposition. However, these models require careful separation of the drug, system and trial design Rostami-Hodjegan, 2012). In addition, PBPK modeling approaches require case examples to build certainty in their performance. Relatively strong confidence is now placed on the application of PBPK modeling when cytochrome P450 (CYP)-mediated metabolism is the dominant route of elimination by many investigators, including, but not restricted to, the regulatory agencies (Wagner et al., 2015). Conversely, despite recent efforts, evidence for the applicability of PBPK for the prediction of transporter-mediated disposition remains less comprehensive, particularly for the kidney (Zamek-Gliszczynski et al., 2013;Jones et al., 2015;Posada et al., 2015;Varma et al., 2015). The study of Hsu et al., (2014) investigated, through simulation, reduced proximal tubule cellularity as one possible mechanism that could cause changes in renal drug secretion in renal impairment. However, the additional plausible underlying mechanisms, namely reduced transporter expression or inhibition of transporters by uremic solutes, have not been explored. Such simulations are required to assess the implications of specific model assumptions on systemic and local drug concentrations predicted by the PBPK models. In addition, although simulations were performed to investigate changes to proximal tubule intracellular drug concentrations in a transporter mediated drug-drug interaction (DDI) scenario, such changes were not investigated in the renal impairment scenario (Hsu et al., 2014).
This study aimed to develop a PBPK kidney model for digoxin and apply the model to investigate potential effects of age and renal impairment on digoxin renal drug disposition. Rationale for selecting digoxin in this study was the large availability of clinical data to support the estimation and verification of mechanistic kidney model parameters. The initial in vitro-in vivo extrapolation (IVIVE) of digoxin renal excretion clearance (CL R ) using reported in vitro OATP4C1 transporter intrinsic clearance (CL int,T ) resulted in under prediction of its CL R . Subsequently, this parameter was estimated from collated digoxin clinical data accounting for the uncertainty in the contribution of glomerular filtration to the observed digoxin CL R . The developed model was subsequently used to simulate changes in digoxin CL R due to aging and moderate and severe renal impairment. In the case of renal impairment, the impact of reduction in OATP4C1, P-gp abundance per million proximal tubule cells, or proximal tubule cellularity on digoxin CL R and its proximal tubule cell concentrations were investigated.

Methods
The overall strategy for the development and application of the PBPK kidney model for digoxin is presented as a workflow diagram in Fig. 1.
Clinical data sources. Digoxin mean plasma concentration-time profiles and pharmacokinetic parameters were collated from the scientific literature. Pharmacokinetic parameters of interest were the area under the curve (AUC) for the plasma concentration-time profile, the intravenous clearance (CL) and oral clearance (CL/F), volume of distribution PBPK Kidney Model of Digoxin and Renal Impairment at steady state, and CL R . Where necessary, data were digitized using GetData Digitizer (version 2, www.getdata-graph-digitizer.com).
Verification of basal PBPK model of digoxin. All simulations presented herein were performed using the SimCYP population-based PBPK simulator software, version 14, release 1 (SimCYP, Sheffield, UK) (Jamei et al., 2009;Jamei et al., 2013). Initial simulations of digoxin plasma concentration-time profiles and CL R were performed using the default "healthy volunteers" population file provided with the SimCYP simulator. Optimized parameters for the full PBPK model of digoxin, including mechanistic models of liver and intestine, were recently published in version 12.2 of the simulator (Neuhoff et al., 2013b). The default compound file provided with version 14.1 of the SimCYP simulator (Supplemental Table S1) was verified against clinical data to ensure consistency between different versions of the software. This verification was performed using several clinical studies, with 10 trials for each set of simulations (Supplemental Table  S2) using the default full-PBPK model (Neuhoff et al., 2013b). Upon visual inspection, the simulated concentration-time profiles as well as key pharmacokinetic parameters were not inconsistent with observed data (see Supplemental Fig. S1).
IVIVE of renal clearance and optimization of tubular secretion in PBPK model. Simulation of renal digoxin disposition was performed using the mechanistic kidney model (MechKiM) module in the SimCYP simulator (Neuhoff et al., 2013a). This model comprises eight segments representing the glomerulus, three subregions of proximal tubule, the loop of Henle, the distal tubule, and the cortical and medullary collecting ducts. Separate compartments represent the blood and tubular filtrate of each segment, as well as the tubular cell mass for the seven tubular segments. The model can account for glomerular filtration, passive permeability within all of the tubular segments, and active transport processes within the proximal tubule segments. Readers may refer to previous publications for model equations and assumptions (Hsu et al., 2014;Burt et al., 2016). Digoxin has low passive permeability across cell monolayers, even in the presence of P-gp inhibitors (literature Caco-2 apparent permeability range from 1.15 to 8.03 Â 10 26 cm/s in various assay formats used (Neuhoff et al., 2003;Zhang and Morris, 2003;Djuv and Nilsen, 2008;Fossati et al., 2008). There are mixed reports from clinical studies suggesting possible urine flow-dependent CL R of digoxin (Steiness, 1974;Halkin et al., 1975;Steiness et al., 1982) and potential role of passive tubular reabsorption in vivo (although minor). Although prediction of tubular reabsorption clearance from physicochemical properties was recently reported, the validated quantitative structurepharmacokinetic property relationship model does not allow for prediction of passive diffusion clearance (CL PD ) values in different regions of the nephron (Dave and Morris, 2015). Recently, a static model for prediction of tubular reabsorption was reported that considered regional differences in physiologic parameters for prediction of fraction reabsorbed; this model has not been validated in the MechKiM using appropriate range of model compounds (Scotcher et al., 2016c). To account for tubular reabsorption in the current model, the same value for CL PD parameter (0.01 ml/min/million tubule cells) was assigned to all tubular compartments (i.e., proximal tubule, loop of Henle, distal tubule, and collecting duct). This assigned value for CL PD resulted in simulated fraction reabsorbed for digoxin of 0.12, in agreement with the range predicted using the static model of tubular reabsorption using the Caco-2 apparent permeability literature permeability data (predicted range using static model 5 0.064-0.34; see Supplemental Fig. S2) (Scotcher et al., 2016c). The assumption of equal CL PD for each tubular region may underestimate the CL PD in proximal tubule, because its surface area is greater than for the other tubular regions (Scotcher et al., 2016c). Such underestimation of CL PD in proximal tubule could subsequently lead to underestimation of OATP4C1 CL int,T . Fraction unbound in kidney cells (f u,kidney,cell ) of 0.51 was predicted using the Rodgers and Rowland method (Rodgers et al., 2005;Rodgers and Rowland, 2006).
Digoxin secretion in kidney is considered to be mediated predominantly by the OATP4C1 and P-gp transporters, expressed on the basolateral and apical proximal tubule membranes, respectively (Tanigawara et al., 1992;Mikkaichi et al., 2004;He et al., 2014;Lee et al., 2014). A wide range of in vitro K m and V max values were available for P-gp in the literature, with two studies reporting in vitro data for OATP4C1 (Supplemental Table S3). For consistency, the P-gp K m (mM) and V max (pmol/min/million cells) parameter values in the mechanistic kidney model were the same as those implemented in the pre-existing mechanistic models for liver and intestine (Troutman and Thakker, 2003;Neuhoff et al., 2013b). For the transporters, sensitivity analyses were performed by modifying the OATP4C1 CL int,T and/or the P-gp relative expression factor (REF) parameters. P-gp mRNA expression data in kidney and Caco-2 cells were used to inform the REF parameter for this transporter. This inherently assumes that mRNA expression at the tissue level is representative of expression at the cellular level. Lack of quantitative information on differences in P-gp protein expression and activity between different tissues at the cellular level precluded a more mechanistic approach. Because of the lack of in vitro V max values reported in the literature, OATP4C1 CL int,T was estimated from reported uptake rate (dX/dt) and initial substrate concentration (C) data (Mikkaichi et al., 2004;Chu et al., 2007), using eq. 1, which assumes that C , , K m .
As the transporter expression/abundance data needed to inform the REF scaling factor were lacking, in vitro OATP4C1 CL int,T was normalized to proximal tubule cellularity (Supplemental Table S3), as previously attempted for OCT2 uptake activity (Li et al., 2014). Normalization was performed by assuming equal expression/activity of OATP4C1 in transfected cells and proximal tubule cells and a cellularity of 13 million Madin-Darby canine kidney cells per milligram protein (Richardson et al., 1981). Although large uncertainty is associated with the number of proximal tubule cells per gram of kidney (PTCPGK) in human (Scotcher et al., 2016b), a value of 60 million cells/g kidney was used for this model parameter, in accordance with previous simulation studies (Neuhoff et al., 2013a;Hsu et al., 2014;Li et al., 2014). As the IVIVE approach outlined above did not successfully predict clinically observed digoxin CL R , OATP4C1 CL int,T was estimated using collated digoxin clinical data. First, parameter estimation was performed by fitting the model to the observed plasma concentration profiles of 9 clinical studies after intravenous administration (Johnson and Bye, 1975;Koup et al., 1975;Ochs et al., 1978;Kramer et al., 1979;Ding et al., 2004), using the automated parameter estimation module of the SimCYP simulator. As simultaneous fitting of all data were not possible in the SimCYP software, fitting was performed using each dataset separately, and the overall weighted (by subject number) mean estimate of OATP4C1 CL int,T was subsequently calculated. In addition, the observed overall weighted (by subject number) mean CL R collated from 19 clinical studies in total (Supplemental Table S4) was used to optimize OATP4C1 CL int,T using a detailed sensitivity analysis. In this approach, digoxin pharmacokinetics was simulated using various values of OATP4C1 CL int,T and serum creatinine in a population representative. The OATP4C1 CL int,T that resulted in the simulated CL R in closest agreement with the observed data when serum creatinine was fixed to 80 mmol/l (corresponding to glomerular filtration rate of 120 ml/min in the population representative) was taken forward as the optimal value. This optimization procedure was performed using two separate clinical study designs that had different doses and some differences in the age range. Significant dose-related differences were not expected, because OATP4C1 CL int,T is assumed to be within the linear range of kinetics (Kramer et al., 1979;Rengelshausen et al., 2003). The final model was verified against plasma and urinary digoxin concentration data from additional clinical studies that were not included in the model development (Johnson and Bye, 1975;Ochs et al., 1978;Lindenbaum et al., 1981). In addition, simulation of digoxin plasma concentration data using the final mechanistic model was compared with digoxin PBPK model where mechanistic kidney model was not included and CL R was defined as a single input parameter.
Simulation of digoxin CL R in elderly and renal impairment populations. Digoxin pharmacokinetics was simulated in 100 virtual subjects after intravenous infusion (30 minutes) of 0.75 mg in different virtual populations, with other trial design parameters such as age range and proportion of women of the virtual subjects determined by the general values of the relevant population files. The virtual populations used for simulations included the "Geriatric NEC", "RenalGFR_30-60", and "RenalGFR_less_30" populations, which were supplied with the SimCYP simulator. Geriatric NEC population accounts for changes in the age-sex distribution, weights, and heights and kidney size of subjects, whereas the RenalGFR_ populations account for changes in glomerular filtration rate (GFR), protein binding by albumin, hematocrit, kidney weight, and renal blood flow in different stages of renal impairment. There are conflicting reports that digoxin volume of distribution may be altered in patients with end-stage renal disease (Jusko and Weintraub, 1974;Cheng et al., 1997); because of the lack of concordant results on the magnitude of such changes and the fact that these patients were not simulated here, this aspect was not captured in the current study.
In addition the "RenalGFR_30-60" and "RenalGFR_less_30" populations were modified by reducing the value of the PTCPGK parameter. Reducing PTCPGK in the model was investigated as a way to mechanistically represent changes in secretion resulting from tubular cell damage, which may occur during renal impairment and injury (and possibly aging) (Nangaku, 2006;Andreucci et al., 2014;Wang et al., 2014). Alternatively, change in expression of drug transporters per million proximal tubule cells was explored as a mechanism contributing to reduced renal secretion, in line with limited biological data (Naud et al., 2011;Wang and Sweet, 2013). Changes to abundance of the OATP4C1 (basolateral membrane) or P-gp transporters (apical membrane) or PTCPGK were each simulated separately. Although changing the PTCPGK parameter or transporter abundance per million proximal tubule cells for a specific transporter will have the same net effect on overall transporter abundance for that specific transporter per individual (after scaling for kidney weight of the individual), some important differences should be noted. Changes in proximal tubule cell number (i.e., PTCPGK) will affect the overall abundance per kidney of all transporters, as well as scaled CL PD per kidney. The latter will not occur when the transporter abundance per million proximal tubule cells is changed. In addition, the MechKiM model uses the PTCPGK parameter to define the cellularity of all tubular regions, including loop of Henle, distal tubule, and collecting duct compartments. Thus changing the PTCPGK can affect the scaled CL PD in the distal tubular regions and therefore the tubular reabsorption, a factor that is not directly affected by changes to transporter abundance parameters. Reduced abundance of kidney drug transporters per million proximal tubule cells was represented in MechKiM by assigning relative abundances for the OATP4C1 and P-gp transporters in kidney in the "poor transporter" (PT) phenotype as a proportion the SimCYP default "extensive transporter" phenotype value of 1 and setting the frequency of PT in the modified population to 1. Separate changes to PTCPGK or transporter abundance parameters applied equally to each of the three subregion of the proximal tubule. Relative abundance of P-gp in liver and gut remained the same. AUC ratio (AUCR) was calculated using eq. 2: where AUC RI and AUC Control are the mean digoxin AUC in renal impairment subjects (GFR , 60 ml/min/1.73 2 ) and mean digoxin AUC in healthy volunteers or patients without renal impairment (GFR . 60 ml/min/1.73 2 ). CL R ratio and maximum concentration of digoxin in the cells of the first of the three proximal tubule segments (C max,PT-1 ) ratio were calculated in an analogous manner. Separately, simulations were performed in the population representative mode after changes in systems parameters in the kidney model, using the "healthy volunteers" population file as a template. In these simulations, digoxin CL R was simulated after changes in GFR either alone or in combination with proportional changes in the OATP4C1 abundance per million proximal tubule cells or PTCPGK parameter (see Table 1 for details).

Results
Optimization of digoxin kidney transporter kinetic parameters. The relevance of individual renal mechanisms for digoxin CL R was assessed using the MechKiM module in a stepwise manner. Consideration of either glomerular filtration in isolation or both glomerular filtration and passive tubular reabsorption resulted in simulated digoxin CL R of 98.8 or 87.3 ml/min, respectively, which were both lower than the overall weighted (by subject number) mean observed value of 136.1 ml/min (Supplemental Table S4). Digoxin plasma concentration time profiles for these scenarios are presented in Fig. 2. Simulated AUC 0-' was 29% and 39% higher compared with setting when MechKiM was not activated (i.e., CL R defined using single input parameter) for the "filtration only" and "filtration and reabsorption" scenarios, respectively.
The sensitivity of simulated digoxin CL R and AUC to model input parameters was assessed to determine their relative Reduction in filtration and secretion was performed to represent changes in renal impairment. Simulated population representative of the "healthy volunteers" population had an age, weight, and body surface area of 20 years, 81 kg, and 1.98 m 2 , respectively. Serum creatinine (input parameter of model) was calculated for each scenario using the Cockcroft-Gault equation (Cockcroft and Gault, 1976), based on the target GFR and the age, weight, and body surface area of the population representative. PBPK Kidney Model of Digoxin and Renal Impairment importance. Simulated digoxin CL R and systemic exposure were highly sensitive to changes in OATP4C1 CL int,T , in contrast to the marginal effect of changes in P-gp REF (Fig. 3A). Simulated C max, PT-1 was sensitive to changes in both OATP4C1 CL int,T and P-gp REF (Fig. 3B). The simulated digoxin CL R and AUC 0-' were insensitive to changes in f u,kidney,cell , with minor changes noted at f u,kidney,cell , 0.2 (Supplemental Fig. S3, A and  B). C max, PT-1 was sensitive to changes in f u,kidney,cell at values below approx. 0.4 (Supplemental Fig. S3C).
Three studies reported P-gp relative mRNA expression between kidney and intestine (Supplemental Table S5). Multiplying these values by the intestine: Caco-2 REF of 2.04 used in the SimCYP gut module resulted in kidney: Caco-2 REFs ranging from 0.78 to 5.34. However, the cellular and tissue expression data were reported in different studies. Therefore, the P-gp REF for Caco-2 cells: kidney used in the final digoxin kidney model (1.51) was based on P-gp mRNA expression data reported in both systems in the same study (Hilgendorf et al., 2007). In vitro OATP4C1-mediated uptake clearance of digoxin was calculated from two uptake rate values reported in the literature (Table 2). In vitro uptake clearance values varied by ∼3 orders of magnitude, and the use of these values to inform the OATP4C1 CL int,T parameter (i.e., IVIVE) resulted in over eightfold difference in simulated CL R (Table 2), representing 67% and 591% of the observed value (Supplemental Table S4).
Alternatively, OATP4C1 CL int,T was obtained by parameter estimation by fitting the model to the observed plasma concentration-time profiles after intravenous administration. Fitting the model separately to the plasma concentration-time profiles from nine clinical studies resulted in weighted mean OATP4C1 CL int,T of 1.85 ml/min/million PTC (Table 2). By using this weighted mean OATP4C1 CL int,T value, simulated digoxin CL R was 80% of the observed value. A sensitivity analysis based approach was next used for estimation of the OATP4C1 CL int,T parameter. The overall weighted mean observed CL R of digoxin obtained from extensive literature search (136.1 ml/min; n 5 214 healthy subjects) was used as the optimal value (Supplemental Table S4). By using a population representative, digoxin pharmacokinetics were simulated using different OATP4C1 CL int,T and serum creatinine input parameter values. The estimated OATP4C1 CL int,T value, based on the sensitivity analysis approach using a fixed serum creatinine value of 80 mmol/l, was 4.14 ml/min/million PTC (Fig. 4). A serum creatinine value of 80 mmol/l was an assumed average value, because 12 of 19 of the clinical studies used in the literature analysis did not report serum creatinine, creatinine clearance (CL CR ), or other measurements/estimates of GFR of subjects enrolled. The comparison of simulated (colored mesh) and observed (gray plane) digoxin CL R in Fig. 4 indicates a range of possible values for optimized OATP4C1 CL int,T (i.e., various intersections between simulated and observed CL R ), depending on the assumed serum creatinine value. From the mean GFR or CL CR data in the 7 digoxin clinical studies that reported these values, the minimum and maximum serum creatinine concentrations (54 and 108 mM) were calculated. Estimation of the OATP4C1 CL int,T parameter assuming these values resulted in twofold differences in the optimized value, as illustrated by Fig. 4. Based on this sensitivity analysis, the estimated OATP4C1 CL int,T value would reliably support conclusions drawn from extrapolation within the context of the current study. All digoxin MechKiM parameters in the developed model are listed in Table 3.
After optimization of the OATP4C1 CL int,T parameter, digoxin pharmacokinetics were then simulated after clinical trial designs reported in studies (Johnson and Bye, 1975;Ochs et al., 1978;Lindenbaum et al., 1981) that were not used in the model development/optimization. Simulated plasma concentration and urinary excretion rate profiles were not inconsistent with the observed data (Fig. 5). In addition, simulated digoxin plasma concentrations were comparable before and after activation of MechKiM to define the CL R of digoxin (Fig. 2).
Simulation of digoxin pharmacokinetics in special populations: effects of age and renal impairment. Mean digoxin CL R simulated in elderly virtual subjects ("Sim-Geriatric NEC" population file) was 90.7 ml/min, 31% lower than simulated CL R in for healthy volunteers (Table 4). This change in digoxin CL R was comparable to the relative change observed in clinical study [36% lower CL R in elderly subjects relative to young subjects (Ewy et al., 1969)], but lower than the relative change in simulated GFR [44% lower in elderly, compared with 54% observed (Ewy et al., 1969)].
Mean predicted digoxin CL R in moderate (GFR_30-60) and severe (GFR_less_30) renal impairment virtual populations were 50% and 67% lower than in the simulated healthy volunteers (Table 4). The mean GFR of the virtual subjects with moderate and severe renal impairment were 64% and 82% lower than in the simulated healthy volunteers. Comparison of the observed clinical data with predicted CL R and GFR in the healthy, moderate renal impairment and severe renal impairment virtual populations are shown in Fig. 6A. Although some agreement between predicted and observed data were noted, there was a general trend of overestimation of CL R in the renal impairment populations. The extent of overestimation of CL R was directly correlated with high expression of OATP4C1 in virtual subjects. Average predicted AUCR in moderate renal impairment (1.4) was in agreement with clinical data (1.3). In contrast, model predicted changes in digoxin systemic exposure in severe renal impairment (predicted AUCR of 1.5) were significantly underestimated in comparison with clinical data (AUCR 5 3.3) (Okada et al., 1978). GFR and OATP4C1 relative abundance per million proximal tubule cells had similar coefficient of determination (R 2 ) of the line of best fit with simulated CL R in healthy virtual subjects (Fig. 6B). In contrast, a stronger correlation between either OATP4C1 relative abundance or PTCPGK and CL R was noted in virtual subjects with renal impairment compared with healthy (data not shown).
Reduction in OATP4C1 relative abundance per million proximal tubule cells in the renal impairment virtual populations (decrease in REF from 1 to 0.125) or scenario with comparable changes in proximal tubule cell number (7.5-60 million proximal tubule cells/g kidney) both resulted in similar predicted impact on digoxin CL R and AUCR (Fig. 7). For example, maximal reduction in PTCPGK or OATP4C1 abundance per million proximal tubule cells simulated in the severe renal impairment population resulted in predicted CL R ratios of 0.164 or 0.152, respectively (Fig. 7). In contrast to systemic exposure, model assumptions used showed differential effect in the predicted digoxin concentration in the proximal tubule. Reduced PTCPGK had minimal impact on simulated digoxin concentrations in proximal tubule cells, whereas reduced OATP4C1 or P-gp abundance per million proximal tubule cells decreased or increased digoxin intra-cellular concentrations, respectively (Fig. 7).
Equal CL PD was assigned for each tubular region, which may underestimate the CL PD in proximal tubule because of its larger surface area. However, the impact was marginal in the case of digoxin, considering low CL PD in proximal tubule compartment relative to transporter kinetic parameters. Simulations performed to assess the potential impact of underestimation of CL PD in proximal tubule demonstrated that increase in this parameter up to fivefold would have minor overall impact on digoxin CL R and C max,PT-1 ratios (Supplemental Fig. S4).
Further simulations were performed in population representative mode, using the "healthy volunteers" population as a template, whereby GFR was altered either alone or alongside proportionally altered OATP4C1 abundance or PTCPGK (Fig.  8). Accounting for changes in tubular secretion in renal impairment, assuming that either OATP4C1 abundance per million proximal tubule cells or PTCPGK are affected proportionally to changes in GFR, resulted in improved agreement between simulated and observed digoxin CL R (Fig. 8).

Discussion
Use of PBPK kidney models is challenged by a lack of physiological data to inform system parameters and scaling  (Greiner et al., 1999). Insets show the graphs presented on logarithmic scales.

PBPK Kidney Model of Digoxin and Renal Impairment
factors, the need for various in vitro data, and availability of suitable clinical data for the drug and population(s) of interest (Scotcher et al., 2016a,b). In an attempt to overcome such challenges in the current study, digoxin was selected as a model drug because of availability of in vitro data and ample clinical data (both plasma and urine) measured in healthy subjects and special populations (e.g., elderly and renal impairment).
a CL R data presented in Supplemental Table S4. Fig. 4. Estimation of OATP4C1 CL int,T parameter using a sensitivity analysis approach by simulating digoxin CL R in population representatives with different serum creatinine values. Colored meshes and gray horizontal plane indicate the simulated CL R and the overall weighted (by subject number) mean CL R obtained from the literature analysis (136.1 ml/min; n = 214 healthy subjects), respectively. Values of OATP4C1 CL int,T and serum creatinine parameters were varied using the automated sensitivity analysis tool in the SimCYP simulator. Optimal OATP4C1 value was taken at the intersection (yellow star) of the simulated digoxin CL R with the observed CL R at a serum creatinine value of 80 mmol/l (which corresponds to simulated GFR ∼120 ml/min), as indicated by the blue arrows. Sensitivity analysis was performed twice using clinical trial designs reported previously (Kramer et al., 1979;Rengelshausen et al., 2003).
Development of mechanistic kidney model for digoxin. A PBPK model for digoxin was previously published and incorporated permeability-limited organ models for the gut and liver (Neuhoff et al., 2013b,c). In the current study the mechanistic kidney model was developed, accounting for the contribution of OATP4C1 and P-gp to digoxin renal secretion (Mikkaichi et al., 2004;Lee et al., 2014). Despite a large amount of available clinical plasma and urine concentration data for digoxin, the P-gp transporter kinetic parameter was practically nonidentifiable because of the lack of measured intracellular digoxin concentrations or sufficiently high temporal resolution of the urinary excretion rate data. As such, accuracy of the IVIVE approach for the P-gp REF parameter could not be assessed, consistent with previous efforts to model renal efflux transport of other drugs (Hsu et al., 2014;Posada et al., 2015). The developed model was therefore unsuitable for  (Hilgendorf et al., 2007) CL PD (ml/min/million PTC) 0.01 Estimated using sensitivity analysis and comparing simulated F reab with that predicted using static tubular reabsorption model (Scotcher et al., 2016c) and published Caco-2 data (Neuhoff et al., 2003;Zhang and Morris, 2003;Djuv and Nilsen, 2008;Fossati et al., 2008). Same value for apical and basolateral membranes in all segments of nephron f u,kidney , fraction unbound in kidney; f u,urine , fraction unbound in urine; F reab , fraction reabsorbed; CL int,T , transporter intrinsic clearance; CL PD , permeability diffusion clearance; K m , Michaelis constant; MechKiM, mechanistic kidney model in SimCYP simulator; PTC, proximal tubule cells; RAF, relative activity factor; REF, relative activity factor; V max , maximal velocity. PBPK Kidney Model of Digoxin and Renal Impairment simulation of P-gp-mediated digoxin DDIs, although sensitivity analyses indicated that inhibition of renal P-gp would unlikely substantially affect digoxin plasma concentrations (Fig. 3), as proposed by Neuhoff et al. (2013b).
There was also insensitivity of digoxin AUC and CL R to changes in f u,kidney,cell (Supplemental Fig. S3) because of very low CL PD relative to the active transport parameters. Changes to simulated unbound digoxin concentrations in the proximal tubule cells would affect the rate of digoxin diffusion back to the systemic circulation. However, because of low passive permeability, this change was marginal compared with uptake via active transport, and therefore AUC and CL R were mostly unaffected.
Reliable estimate of the OATP4C1 CL int,T could not be obtained based on reported in vitro transporter data and IVIVE scaling factors. Use of in vitro uptake transporter kinetic data from different literature sources resulted in approximately eight-fold difference in predicted CL R (Table 2). This finding further indicates the need for quality in vitro kinetic data and transporter abundance for in vitro systems and human kidney.
The weighted mean OATP4C1 CL int,T obtained by fitting to clinical plasma concentration data were approximately half of the value obtained by using CL R data with a sensitivity analysis ( Table 2). The resultant simulated digoxin CL R values were each within 1.5-fold of the observed value. The use of CL R data for optimization was viewed as more reliable because this approach focuses on the parameter of interest with less noise than when using plasma concentration data. Furthermore, separate fitting to mean plasma digoxin concentration data from individual studies was performed rather than global fitting, which can lead to bias in parameter estimation.
Nevertheless, sources of uncertainty of the optimized OATP4C1 CL int,T parameter should be noted. First, the "optimal" CL R value used was obtained from a meta-analysis of clinical data from 19 studies (weighted mean CL R 5 136.1 ml/min). This value differed slightly from those obtained in previous literature analyses [151.1 ml/min (Scotcher et al., 2016c) and 160.8 ml/min (Neuhoff et al., 2013b)], although in the current study some data were assigned for model evaluation. Second, the estimated OATP4C1 CL int,T was dependent on the serum creatinine parameter (Fig. 4). Ideally the serum creatinine value would be informed by data obtained from subjects participating in the specific clinical studies, although in the current study, suitable data were missing in the majority of clinical study reports used in the CL R data analysis.
Simulation of digoxin renal drug disposition in renal impairment: implications for drug toxicity. For drugs eliminated predominantly by renal excretion, dosage adjustment (e.g., for elderly patients or with impaired renal function) is informed by the ratio of the estimated GFR (eGFR) or CL CR in patients relative to subjects with normal renal function (Elinder et al., 2014). In the current study the difference in predicted digoxin CL R between elderly and young virtual subjects was smaller (31%) than corresponding differences in simulated GFR [44%; calculated using the Cockcroft-Gault equation (Cockcroft and Gault, 1976)], in agreement with the findings of a clinical study (Ewy et al., 1969). This supports the proposal that despite physiologic changes in kidney during aging (Darmady et al., 1973), proximal tubule secretion is largely retained in elderly subjects without kidney disease (Musso and Oreopoulos, 2011).
Dose adjustment in the clinic tends to use a reduction in renal plasma clearance by a ratio equivalent to the ratio of eGFR (or estimated CL CR ) in renally impaired patient compared with  (Bloom et al., 1966;Okada et al., 1978). Solid black line represents linear line of best fit using total least squares regression (which recognizes experimental error in both variables) for the observed clinical data. (B) Simulated CL R and GFR in healthy (yellow solid squares) and moderate (green open circles) and severe (turquoise solid diamond) renal impairment virtual subjects. Solid lines represent linear lines of best fit using ordinary least squares regression for data from each simulation, with relevant equations and R 2 shown in boxes.

492
patient with normal eGFR (Bloom et al., 1966;Okada et al., 1978). This assumes that all processes responsible for renal handling, including tubular secretion, decline in parallel with GFR (Bricker et al., 1960;Naud et al., 2011;Schnaper, 2014). As such, when active secretion occurs, mechanistic models that account only for changes in GFR cannot accurately simulate the decline in drug CL R (Grillo et al., 2012;Hsu et al., 2014;Li et al., 2014). A caveat is that simulated GFR values are sometimes necessarily and inappropriately compared with observed CL CR data (Bauer et al., 1982;Lin et al., 2013). Various underlying physiological changes have been proposed to cause reduced tubular secretion in renal impairment, including transporter inhibition by uremic solutes, loss of proximal tubule cells, and decreases in transporter expression levels (Naud et al., 2011;Hsueh et al., 2016). Therefore, to mimic renal impairment, changes to transporter abundance (OATP4C1 and P-gp) per million proximal tubule cells and proximal tubule cellularity parameters of the model were considered in the current study. The effect of uremic toxins was not investigated due to limited availability of inhibition data. Equivalent changes to the OATP4C1 abundance or PTCPGK parameters had comparable impact on the predicted AUCR (e.g., reduction of either parameter by 50% in severe renal impairment resulted in 5% increase in AUCR) and CL R ratio (e.g., reducing OATP4C1 abundance per million proximal tubule cells or PTCPGK by 50% in severe renal impairment resulted in 31% or 27% decrease in CL R ratio, respectively). The minor differences between changing OATP4C1 abundance and PTCPGK occurred because the PTCPGK parameter affects also tubular reabsorption, which is not affected by the OATP4C1 abundance parameter (see Methods). In contrast, changes in renal P-gp abundance per million proximal tubule cells had negligible effect on these pharmacokinetic parameters for digoxin (Fig. 7).
Conversely, simulated digoxin C max,PT-1 was insensitive to changes in proximal tubule cell number but was affected by changes in the transporter abundance parameters (Fig. 7). Overall, the results show that the high degree of correlation between OATP4C1 abundance per million proximal tubule cells and PTCPGK with respect to effect on systemic exposure is not apparent when the dynamic situation within the proximal tubule cell is considered. Therefore, improved understanding of underlying mechanisms behind changes in tubular secretion in renally impaired patients is crucial to determine the increased risk of proximal tubule drug-related toxicity reported in such patients (Naughton, 2008) and for projecting the combined impact of multiple factors (e.g., transporter mediated DDIs in renally impaired subjects). Fig. 8. Simulation of digoxin CL R in population representative mode with changes in different systems parameters performed to represent changes in the case of renal impairment. Glomerular filtration rate (range 20-140 ml/min/1.73 m 2 ) was changed by altering the serum creatinine parameter (74.5-695.6 mmol/l); OATP4C1 abundance and PTCPGK parameters were altered by a factor proportional to the relative change in GFR from the population representative of the default "healthy volunteers" population (GFR = 136.4 ml/min/1.73 m 2 ; serum creatinine = 76.5 mmol/l). Lines represent simulations performed with changes in GFR alone (orange dashed line), both GFR and OATP4C1 abundance (blue dashed line), or both GFR and PTCPGK (green dashed lines). Reported clinical data (black X) are overlaid (Bloom et al., 1966;Okada et al., 1978). Fig. 7. Impact of reduced renal secretion on simulated digoxin AUC ratio (A), CL R (B) and C max,PT-1 ratio (C) in renal impairment populations. Renal secretion was reduced either by changing the kidney OATP4C1 or P-gp relative abundance parameters or by reducing the PTCPGK parameter by a proportional amount. Lines represent changes in PTCPGK in moderate renal impairment (light green solid line) and severe renal impairment (green dashed line), OATP4C1 abundance in moderate renal impairment (light blue solid line) and severe renal impairment (blue dashed line), and P-gp abundance in moderate renal impairment (light purple solid line) and severe renal impairment (purple dashed line). Each scenario was simulated in 100 virtual subjects. Solid horizontal black line (ratio = 1) represents the healthy volunteer population; estimated CL R ratios for the average moderate (GFR = 46.5 ml/min/1.73 m 2 ; orange solid line) and severe (GFR = 23.5 ml/ min/1.73 m 2 ; orange dashed line) renal impairment were calculated based on the correlation of GFR and CL R in the observed data (Bloom et al., 1966;Okada et al., 1978) (A). Relative change of PTCPGK or transporter abundance of 1 indicates that the default moderate or severe renal impairment population in the SimCYP simulator was used PBPK Kidney Model of Digoxin and Renal Impairment Conclusion A mechanistic kidney model for digoxin was developed accounting for the roles of OATP4C1 and P-gp in its tubular secretion and subsequently applied for the prediction of digoxin pharmacokinetics in special populations. Consideration of reduced GFR in renal impairment in isolation was insufficient to capture changes in digoxin CL R in this patient group. Different mechanisms associated with reduced active tubular secretion in renal impairment were explored in the kidney model, namely reduced transporter abundance per million proximal tubule cells and decrease in proximal tubule cellularity. Although quantitative transporter abundance data in normal human kidney samples are emerging and are likely to become available for diseased tissue, data on proximal tubule cellularity in normal and diseased kidneys are lacking. Reduction in OATP4C1 expression or PTCPGK each caused comparable changes on the predicted digoxin systemic exposure and CL R . In contrast, predicted proximal tubule concentration of digoxin was only sensitive to changes in the transporter expression parameters. These results suggest that depending on the output parameter of interest, accurate model specification of pathophysiology may or may not be important. However, the implications of potential misspecification could be more severe if the model developed is applied to extrapolate to more complex scenarios (e.g., transporter-mediated DDIs in renally impaired patients).