Revisiting Nonlinear Bosentan Pharmacokinetics by Physiologically Based Pharmacokinetic Modeling: Target Binding, Albeit Not a Major Contributor to Nonlinearity, Can Offer Prediction of Target Occupancy s

Bosentan is a high-affinity antagonist of endothelin receptors and one of the earliest examples for target-mediated drug disposition [a typeofnonlinearpharmacokinetics(PKs)causedbysaturabletarget binding]. The previous physiologically based PK (PBPK) modeling indicated that the nonlinear PKs of bosentan was explainable by considering saturable hepatic uptake. However, it remained un-examined to what extent the saturable target binding contributes to the nonlinear PKs of bosentan. Here, we developed a PBPK model incorporating saturable target binding and hepatic uptake and analyzed the clinical bosentan PK data using the cluster Gauss-Newtonmethod(CGNM).ThePBPKmodelwithouttargetbindingfell shortincapturingthebosentanconcentrationsbelow100nM,based on the PK profiles and the goodness-of-fit plot. Both global and local identifiability analyses (using the CGNM and Fisher information respectively) informed target binding lowest bosentan contributor to nonlinear bosentan PKs, target occupancy PK achieving for efficacy dose cover


Introduction
Bosentan (Tracleer), an endothelin receptor antagonist, is used for the treatment of pulmonary arterial hypertension. The pharmacological target of bosentan is abundantly expressed in vascular endothelium and smooth muscle (Kuc and Davenport, 2004;Maguire and Davenport, 2015). Bosentan displays nonlinear pharmacokinetics (PKs), and in fact it was one of the earliest examples that initiated the concept of target-mediated drug disposition (TMDD, a type of nonlinear PK caused by high-affinity binding of a drug to its pharmacological target) (Mager and Jusko, 2001). When 10-750 mg bosentan doses were administered intravenously to healthy volunteers, its systemic clearance and volume of distribution at steady state decreased with escalating doses (Weber et al., 1996). When bosentan was administered orally to healthy volunteers, the dosenormalized C max values of bosentan were comparable for doses up to 600 mg, but decreased at doses above 600 mg, suggesting that the dissolution rate limited absorption from an aqueous suspension of bosentan (Weber et al., 1996).
The mechanisms of disposition of bosentan have been extensively investigated. After intravenous administration of radiolabeled bosentan, 92.9% and 5.2% of the radioactivity were collected from feces and urine, respectively, and unchanged bosentan in feces accounted for 3.7% of the dose, indicating that most of bosentan undergoes metabolism and subsequent excretion to feces (Tracleer package insert). Bosentan is metabolized mainly by cytochrome P450s 3A4 and 2C9 and also s This article has supplemental material available at dmd.aspetjournals.org.
ABBREVIATIONS: CGNM, cluster Gauss-Newton method; CL, clearance; detFIM, determinant of the Fisher information matrix; FIM, Fisher information matrix; f u,b , fraction unbound in the blood; K d , dissociation equilibrium constant; K m,act,inf , Michaelis-Menten constant for active influx; k off , dissociation rate constant; k on , association rate constant; PBPK, physiologically based pharmacokinetic; PK, pharmacokinetic; SSR, sum of squared residuals; TMDD, target-mediated drug disposition; V d,ss , volume of distribution at steady state; X TotalR , total amount of the target receptor.
handled by the organic anion transporting polypeptide 1Bs (Dingemanse and Van Giersbergen, 2004;Treiber et al., 2007). With the intent of probing the mechanisms involved in the nonlinear PKs of bosentan after intravenous administration, a previous report compared the performances of the physiologically based PK (PBPK) models that incorporated the saturable components for hepatic uptake, hepatic metabolism, or both (Sato et al., 2018). The results indicated that the incorporation of saturable hepatic uptake (but not saturable hepatic metabolism) can describe the nonlinearity observed in the clinical PK data of bosentan (Sato et al., 2018). Other researchers examined the contribution of saturable target binding to the nonlinear PKs of bosentan. Volz et al. (2017) reported a compartment model incorporating TMDD, which included the formation and internalization of the drug-receptor complex as a clearance mechanism for bosentan. But the authors did not consider any other saturable processes involved in bosentan disposition. Another PBPK model reported by Li et al. (2018) considered nonlinear components of hepatic uptake, hepatic metabolism, target binding (but with no internalization of the drug-receptor complex), and enteric metabolism. Their comprehensive PBPK modeling successfully described various clinical scenarios of bosentan use, but missing in their analysis was the assessment of the relative contribution of saturable target binding to the nonlinear PKs of bosentan.
This study aimed to assess the extent to which the saturable target binding contributes to the nonlinear PKs of bosentan (relative to the saturable hepatic uptake) by developing a PBPK-TMDD model that incorporates saturable target binding and hepatic uptake. Our PBPK-TMDD model incorporating target binding yielded practically identifiable target binding parameters and estimated the target occupancy-time profiles by bosentan, based on the analysis of the blood PK profiles alone. This type of analysis may predict target occupancy by a drug in vivo and potentially offer guidance on achieving optimal efficacy outcomes, under the condition when the high-affinity drug target is responsible for the efficacy of interest and when the dose ranges cover varying degrees of target binding.

Materials and Methods
Structure of the Bosentan PBPK Model. As shown in Fig. 1, our PBPK model of bosentan incorporated nonlinear components for hepatic uptake and target binding and linear components for hepatic/intestinal metabolism. The structural base of our PBPK model was similar to that reported previously (Sato et al., 2018), with the modifications of adding target binding and absorption components. The central compartment was connected to three tissues (adipose, muscle, and skin) that can contribute to the distribution volume of bosentan. Rapid equilibrium between tissues and blood was assumed in our model, and the tissue/ blood drug concentration ratios were calculated as described previously (Rodgers and Rowland, 2006;Sato et al., 2018).
The current study analyzed the reported clinical PK data after the intravenous doses (10, 50, 250, 500, and 750 mg) (Weber et al., 1996) and the oral doses (62.5, 125, and 500 mg) (Weber et al., 1999;Dingemanse et al., 2002;van Giersbergen et al., 2002) as listed in Supplemental Table 1. For oral dosing, the intestinal absorption model assumed a linear process when considering the undissolved drug, dissolved drug, and the enterocyte compartment and used the same rate constant for dissolution and absorption, similar to the previous report (Li et al., 2018). Based on the previous findings that the nonlinear bosentan PKs were explainable by incorporating saturable hepatic uptake (Sato et al., 2018), the current model employed saturable hepatic uptake and nonsaturable hepatic metabolism. Target binding was modeled to be connected directly to the central compartment, described using the dissociation rate constant (k off ) and the dissociation equilibrium constant [K d ; the association rate constant (k on ) was defined as k off /K d ]. Receptor occupancy was defined as the amount of receptordrug complex divided by the total amount of receptor (X TotalR ). In total, 30 parameters were included in our PBPK model. The model equations are listed in the Supplemental File. As we had no access to individual-level data, we did not include interindividual variability in our model.
Parameters for the Bosentan PBPK Model. Ten of the total 30 parameters in the PBPK model were estimated, whereas the remaining 20 parameters were fixed based on information available in the literature (Supplemental Table 2). Physiologic parameters were calculated using a body weight of 78 kg, which was Fig. 1. Structure of the bosentan PBPK model. The PBPK model incorporates saturable components for bosentan hepatic uptake and binding to its pharmacological target and linear components for the intestinal absorption (limited by the dissolution rate) and hepatic metabolism. Target binding is modeled to be directly connected to the central compartment and described using k off and K d (k on = k off /K d ). Detailed information including the ordinary differential equations and the nomenclature is available in Supplemental Material. CL met , hepatic metabolic clearance; CL r , renal clearance; D i.v. , the dose administered intravenously; D p.o. , the dose administered orally; EH i (i = 1-5), the ith extrahepatic compartment; F a , the fraction of the dose absorbed; f u,ent , fraction of the unbound drug in the enterocyte compartment; f u,h , the fraction of the unbound drug in the hepatocellular compartment; HC i (i = 1-5), the ith hepatocellular compartment; k a , the absorption rate constant; K p,a , adipose-to-blood drug concentration ratio; K p,m , muscle-to-blood drug concentration ratio; K p,s , skin-to-blood drug concentration ratio; V max,act,inf , the maximum velocity for the active influx through the basolateral membrane of hepatocytes; PS ent,b,eff , the efflux intrinsic clearance via passive diffusion through the basolateral membrane of enterocytes; PS dif,eff , the efflux intrinsic clearance via passive diffusion through the sinusoidal membrane of hepatocytes; PS dif,inf , the influx intrinsic clearance via passive diffusion through the sinusoidal membrane of hepatocytes; Q h , Q m , Q s , and Q a , the blood flow rate to the liver, the muscle, the skin, and to the adipose, respectively; R met , the ratio of intestinal/hepatic metabolic clearance; V a , V m , and V s , the volume of the adipose compartment, the muscle compartment, and the skin compartment, respectively; V central , V ent , V eh , and V hc , the volume of the central compartment, the enterocyte compartment, the extrahepatic compartment, and the hepatocellular compartment, respectively.
Target Occupancy of Bosentan Estimated by PBPK Modeling 299 at ASPET Journals on July 16, 2021 dmd.aspetjournals.org the mean value for the participants in the clinical study of intravenous administration (Weber et al., 1996), and tissue volume was converted from tissue weight (Davies and Morris, 1993) by assuming a tissue density of 1 g/ml. The fraction unbound in the blood (f u,b = 0.033) was calculated from that in plasma (0.02) and the blood/plasma concentration ratio of 0.6 reported in the Tracleer package insert. The fraction unbound in hepatocytes was from the previously reported in vitro study on ice (Sato et al., 2018). The renal clearance of bosentan was calculated using the fraction of the drug excreted unchanged in the urine of 0.008 in humans (Weber et al., 1996). The ratio of intestinal/hepatic metabolic clearance was reported previously (Li et al., 2018). The fraction absorbed was assumed to be unity. The fraction unbound in enterocytes was also assumed to be unity, based on the results from the previous report (an optimized median value of 0.964; 95% confidence interval of 0.50 and 1.0) (Li et al., 2018). Efflux clearance of bosentan via passive diffusion in the liver was calculated by dividing influx clearance via passive diffusion by 0.243 (g value) as described previously (Yoshikado et al., 2016;Sato et al., 2018).
Parameter Estimation and Global Practical Identifiability Analyses Using the Cluster Gauss-Newton Method. The ten unknown parameters were estimated using the cluster Gauss-Newton method (CGNM), which samples multiple solutions of nonlinear least-squares problems (Aoki et al., 2020); i.e., CGNM finds multiple parameter sets that minimize the sum of squared residuals (SSR). The SSR as an objective function was defined by the following equation: SSR ¼ + n i¼1 ðlog 10 y obs;i 2 log 10 y pred;i Þ 2 where y obs,i and y pred,i are the ith observed value and model prediction, respectively. The key characteristic of CGNM is that it simultaneously estimates parameters starting from multiple initial estimates. (We shall refer to this collection of initial estimates as the "initial cluster.") In generating an initial cluster for CGNM, we set the base values for each unknown parameter based on the reported values in the literature, and the initial lower and upper boundaries were set to the base values multiplied by 10 22 and 10 2 , respectively (Table 1). Within the initial range, we generated 1000 parameter sets as an initial cluster, and the parameter sets were improved by 100 iterations. The CGNM calculation was run on CGNM software version a.4.4.0 (http://www.bluetree.me/CGNmethod_ for_PBPKmodels/). The CGNM is a computationally efficient algorithm to estimate multiple parameters by starting a conventional parameter estimation method with multiple initial estimates. Unlike the conventional parameter estimation methods, the CGNM assumes that the parameter cannot be uniquely determined from the observation (i.e., there is no unique best fit parameter).
Conversely, if the multiple parameter sets found by the CGNM converge to one parameter set, the parameter can be deemed to be globally practically identifiable (the parameter can be uniquely determined from the given experimental data) (Aoki et al., 2020).
Local Parameter Identifiability Analyses Using Fisher Information. For local practical parameter identifiability analysis, we computed Jacobian matrices of the simulation with respect to the estimated parameters at the best fit parameter sets found by the CGNM using finite difference scheme. The observed Fisher information matrix (FIM) was calculated using the Jacobian, and we used the determinant of the FIM (detFIM) to quantify the information content for the three target binding-related parameters (K d , k off , and X TotalR ) (Fedorov, 1972;Aoki et al., 2016). The larger detFIM, the more information content to estimate the parameters, and when detFIM = 0, then one or more parameters are not locally practically identifiable.

PBPK Modeling of Bosentan and Simulation of Target
Occupancy. When sorted in ascending order, the SSR values of the optimization results (1000 parameter sets) displayed minimal changes from ranks 1 to 401 (the SSR values ranging from 2.21 to 2.30), followed by an abrupt increase from rank 402 (Fig. 2, inset). The parameter sets from ranks 1 to 401 were thus chosen for further analysis. These parameter sets (ranks 1-401) were used to draw the blood bosentan concentration-time profiles, which were in good agreement with the observed data and very close to each other and appeared overlapping at each dose (Fig. 3). For 10 of the estimated parameters, the calculated values from ranks 1 to 401 fell within narrow ranges, indicating that these parameters could be practically identified from the given clinical data ( Fig. 4; Table 1). When the procedure was repeated five times using different random seeds, all five runs yielded comparable outcomes in terms of the SSR values (Supplemental Fig. 1) and the estimated values for individual parameters (Supplemental Fig. 2).
The parameters from ranks 1 to 401 were used to simulate the time profile of receptor occupancy for each dose level. The 401 simulation results also fell within narrow ranges, but a slight degree of variation was noticeable, as shown in shaded regions (Fig. 5). The maximum values of receptor occupancy for single oral clinical doses b, a hybrid parameter to determine the rate-determining step in the hepatic elimination; CL int,all , an overall intrinsic clearance of the liver; CL met , hepatic metabolic clearance; k a , the absorption rate constant; V max,act,inf , the maximum velocity for the active influx through the basolateral membrane of hepatocytes; PS dif,inf , the influx intrinsic clearance via passive diffusion through the sinusoidal membrane of hepatocytes; PS ent,b,eff , the efflux intrinsic clearance via passive diffusion through the basolateral membrane of enterocytes; V central , the volume of the central compartment. Secondary parameters (k on , b, and CL int,all ) were calculated using the following equations: k on = k off /K d ; b = CL met /(CL met + PS dif,eff ); CL int,all = (PS dif,inf + V max,act,inf /K m,act,inf ) Â b, where PS dif,eff is the efflux intrinsic clearance via passive diffusion through the sinusoidal membrane of hepatocytes. Model Fit With or Without Target Binding. To assess the impact of including saturable target binding on capturing the clinical bosentan data, we compared the PBPK models with and without target binding. For the PBPK model with target binding, the minimum SSR value was 2.21 (hence, the Akaike information criterion of 2410). For the PBPK model without target binding, the minimum SSR value was 4.53 (hence, the Akaike information criterion of 2337) (Supplemental Fig. 3A). As shown in Fig. 6A (PK profiles) and Fig. 6B (goodness-of-fit plot), the model without target binding has significant model misspecification at low bosentan concentrations, compared with the model with target binding. For most of the optimized parameters, the estimated values were similar between the models with and without target binding (Supplemental Table 3). We calculated the volume of distribution at steady state (V d,ss ) and clearance (CL) values based on the obtained PK profiles using the PBPK model with or without target binding and the parameter set with the smallest SSR values (rank 1). As shown in Fig. 6C, the PBPK model with target binding (blue dotted lines) yielded the V d,ss and CL values that were in better agreement with the reported data than the PBPK model without target binding (red dotted lines).

Parameter Identifiability With or Without Low-Dose Data.
To clarify the importance of the low-dose data in determining parameter values, we estimated the parameters using the CGNM and the clinical bosentan PK data that excluded the data for 10 and 50 mg doses. As shown in Fig. 7, Supplemental Fig. 4, and Supplemental Table 4, the CGNM found multiple parameter sets that minimized the SSR and captured the overall profiles of the observed clinical data, but the parameter values varied to a much greater extent than the parameter sets estimated from the data including low-dose data. These results suggest that the parameters were not uniquely identifiable without the low-dose data. To quantify the information content of the data set, the FIM values were calculated ( Table 2). The determinant of FIM was markedly reduced (approximately 0 when back-transformed from natural logarithmic values) when 10 mg dose data were excluded, indicating that the target binding-related parameters were not locally practically identifiable without 10 mg dose data. The SSR values of 1000 parameter sets (100 iterations) are plotted in ascending order. Red dots indicate SSR values less than 2.3 from ranks 1 to 401, which were then used for further analysis. The best SSR (rank 1) was 2.21.  Plot showing a range of the estimated parameters from ranks 1 to 401. The parameter set of rank 1 is shown as a red line and those of ranks 2 to 401 are shown as black lines that fall within a very narrow range and appear overlapping. Initial ranges for the parameters (10 22 -10 2 -fold ranges from the base values) are depicted as a shaded area. CL met , hepatic metabolic clearance; k a , the absorption rate constant; PS dif,inf , the influx intrinsic clearance via passive diffusion through the sinusoidal membrane of hepatocytes; PS ent,b,eff , the efflux intrinsic clearance via passive diffusion through the basolateral membrane of enterocytes; V max,act,inf , the maximum velocity for the active influx through the basolateral membrane of hepatocytes.
Target Occupancy of Bosentan Estimated by PBPK Modeling 301 at ASPET Journals on July 16, 2021 dmd.aspetjournals.org

Discussion
The current study aimed to assess the extent to which saturable target binding contributes to the nonlinear PKs of bosentan. Building upon the previous PBPK model incorporating saturable hepatic uptake (Sato et al., 2018), we constructed a PBPK model of bosentan that incorporated saturable target binding and hepatic uptake and assessed the performance of the PBPK model in capturing the clinical PKs of bosentan. Our results indicated that target binding is not a major contributor to the nonlinear bosentan PKs. However, the incorporation of target binding improved the goodness of fit, especially for blood bosentan concentrations below 100 nM (Fig. 3 vs. Fig. 6). Per our results from the global and local identifiability analyses (using the CGNM and FIM, respectively), the three target binding-related parameters (K d , k off , and X TotalR ) were practically identifiable from the analysis of blood bosentan PK profiles alone ( Fig. 4; Table 1). The current analyses yielded reasonable predictions of the maximum target occupancies of 0.6-0.8 for the clinical bosentan oral doses (62.5 and 125 mg) (Fig. 5). The availability and inclusion of the low-dose PK data were important in identifying the target binding-related parameters and thereby predicting target occupancy in vivo ( Fig. 7; Table 2). It awaits further validation to support the utility of the current approach in predicting in vivo target occupancy for other drugs. Overall, our results indicate that target binding, albeit not a major contributor to the nonlinear bosentan PKs, may offer a prediction of target occupancy from blood bosentan PK profiles alone. Incorporation of saturable target binding in the PBPK modeling efforts may provide potential guidance on achieving optimal efficacy outcomes, under the condition when the high-affinity drug target is responsible for the efficacy of interest and when the dose ranges cover varying degrees of target binding.
When target binding was not considered, the model reasonably captured the blood bosentan concentrations at the intravenous doses ranging from 250 to 750 mg but underestimated the blood bosentan concentrations at the intravenous doses of 10 and 50 mg (Fig. 6A). In our PBPK modeling, the K d value (inversely related to the binding affinity of unbound bosentan to its target) was estimated to be 11 nM, falling in the ranges of the reported values of K d from in vitro experiments (3-100 nM in human tissues; Supplemental Table 5). The estimated K d value corresponds to 333 nM as the total bosentan concentration in the blood (K d /f u,b = 11/0.033 nM). For an intravenous dose of 10 mg, the bosentan blood concentrations were below 333 nM (K d /f u,b ) after 15 minutes after dosing (Fig. 3A). As such, it is not entirely surprising that the availability and inclusion of the low-dose PK data (below target saturation; e.g., 10 mg dose) are essential in estimating the target binding-related parameters (as quantified by calculating the FIM values; Table 2). When a drug interacts with a highly specific and abundant target, its target binding may represent a route by which a considerable fraction of the drug distributes to the body (the volume of distribution mediated by target binding, which can be described as X TotalR /K d , based on the receptor-ligand binding theory). If the drug displays low lipophilicity and limited nonspecific (non-target-mediated) tissue binding, the volume of distribution mediated by target binding (corresponding to X TotalR /K d ) may be ascertained by analyzing the systemic PK data from low doses (well below target saturation, for instance, microdoses) together with those from intermediate doses (near target saturation) as well as high doses (at target saturation) (Burt et al., 2020).
On the other hand, the parameter values for hepatic active influx were minimally impacted by incorporating target binding (Supplemental Table 3). The estimated values of the Michaelis-Menten constant for the active influx through the basolateral membrane of hepatocytes (K m,act,inf ) were 0.713 and 0.555 mM for the PBPK model with and without target binding, respectively (K m,act,inf /f u,b corresponding to 21.6 and 16.8 mM). When we simulated the hepatic blood concentrationtime profiles of bosentan using the estimated parameters (Supplemental Fig. 5), the bosentan concentration of the first hepatic blood compartment were above K m,act,inf /f u,b (21.6 mM) for intravenous doses of 250 mg and higher, but not for any of the oral doses. These results indicate that the nonlinearity observed at high intravenous bosentan doses can be explained mainly by the saturation of hepatic uptake.
Recently, our group and others reported the findings from their modeling efforts to analyze the nonlinear PK data of bosentan, warranting some discussion on how they differ. Our group previously reported that saturation of hepatic uptake (rather than hepatic metabolism) likely accounts for the PK nonlinearity of bosentan by comparing the performances of the PBPK models that incorporated the saturable components for hepatic uptake, hepatic metabolism, or both (Sato et al., 2018). In that report, no consideration was however given to saturable Fig. 5. Temporal profile of receptor occupancy after intravenous and oral administration of bosentan. Receptor occupancy is defined as the amount of the drug-receptor complex divided by the amount of total receptor in the body. Estimated parameters from ranks 1 to 401 were used to simulate temporal profiles of receptor occupancy (each set of estimated parameters was used for simulation, but the results fell within a very narrow range and appear overlapping). The receptor occupancy was defined as the molar ratio of receptor-drug complex to the total (free and bosentan-bound) receptor. Dark lines indicate the simulation results using the parameter sets of rank 1, whereas light, largely overlapping lines indicate those using the parameter sets of ranks 2-401.  Li et al. (2018) was comprehensive regarding saturable components incorporated (including hepatic uptake, metabolism, target binding, and other processes) and clinical scenarios analyzed. But the authors did not extend their modeling efforts to the prediction of target occupancy. In that report, the reversible, saturable binding of bosentan to its target was assumed to occur in the systemic plasma compartment and not serve as a drug clearance mechanism. The former assumption was based on the fact that endothelin receptors are mainly located in the vasculature (endothelial and smooth muscle cells, in close contact with the systemic compartment), likely with a minimal contribution of target binding in other tissue compartments. The latter assumption was based on the lack of experimental data supporting that the bosentan-target complex could eliminate bosentan, similar to the general framework laid out in the early study (Mager and Jusko, 2001). On the other hand, Volz et al. (2017) reported a compartmental PK model that incorporated saturable target binding (but no consideration given to other saturable processes) and assumed the formation of the bosentan-target complex as a clearance mechanism. In describing the receptor interactions, the authors also included the parameters for degradation and synthesis of the unoccupied receptors. As far as we are aware, the evidence is scarce to support the temporal variations in the level of the unoccupied endothelin receptors. Thus, our current analysis constructed the simplest model that can incorporate saturable target binding consistent with the currently available experimental evidence by modifying the previous PBPK model (Sato et al., 2018). Our current PBPK model of bosentan incorporated the fivecompartment liver model, similar to the previous report (Sato et al., 2018). Whereas the one-compartment liver model is likely sufficient to capture the hepatic disposition of relatively low clearance drugs, the use of a multiple-compartment liver model would be necessary to capture the hepatic disposition of high clearance drugs. With the intent of establishing a flexible and versatile model structure that can accommodate drugs with varying ranges of clearance (i.e., not requiring drug-dependent changes), our current PBPK model incorporated the five-compartment liver model. When the analyses were performed using the PBPK model with the one-compartment liver model, the estimated parameter values were mostly comparable to those obtained using the five-compartment liver model (Supplemental Table 3). Estimated parameters from ranks 1 to 790 were used to draw bosentan blood concentration-time profiles (shown as 790 lines that appear overlapped) along with the observed data (shown as circles). (B) Goodness-of-fit plot of the model with and without target binding. The residual (the difference between the log-transformed observed values and PBPK model-based values) plotted against the logtransformed observed values. The trend curves are local polynomial regression fit with 95% confidence interval. (C) Dose-dependent changes in the observed and PBPK model-based V d,ss and CL values after the intravenous administration of varying bosentan doses. The observed data are shown in closed circles (from Weber et al., 1996). Dotted lines represent the trend for the V d,ss and CL values calculated from the PK profiles drawn using the parameter set with the smallest SSR values (rank 1) (blue dotted lines, from the PBPK model with target binding; red dotted lines, from the PBPK model without target binding). In conclusion, our study indicated that the saturation of hepatic active influx is likely a major contributor to nonlinear bosentan PK profiles observed at high doses and that the saturation of target binding, albeit not a major contributor to the PK nonlinearity, was evident at low doses of intravenous bosentan administration. The bosentan PBPK model incorporating target binding could produce a reasonable estimate of target occupancy based only on the blood concentration-time profiles. In the future, this approach of estimating drug target occupancy from its blood concentration based on a PBPK model incorporating target binding might be of use in predicting the drug efficacy in the early phase of drug development.