Several values related to physicochemical and ADME properties of BIBR0951 and BIBR1087, including plasma stability, plasma protein binding, blood-to-plasma partitioning, enzyme kinetic, and permeability/transport) were determined by in vitro experiments, as described in the section of Supplementary Methods.

Simcyp population based ADME simulator software version 20 (Certara, Sheffield, United Kingdom) was used for the PBPK analysis in this study. All simulations were performed using a virtual “Sim-Healthy Volunteer” population aged 20–50 years for 100 subjects (10 subjects per trial, 10 trials per run). Other parameters for the virtual trial design were set to mimic the clinical designs of observed data as much as possible.

Supplementary Table S1 summarizes the observed datasets for model development and qualification. The clinical PK profiles and parameters of free DAB (unconjugated DAB) were gathered and included in this study. All observed plasma concentration-time profiles data were digitized from the published literature using DigitizeIt software version 2.5.9 (I. Bormann, Braunschweig, Germany). The full plasma concentration-time profiles and associated PK parameters of DABE and two intermediate metabolites (BIBR0951 and BIBR1087) were not available in the literature.

Figure 1 illustrates the stepwise processes for PBPK model development, qualification, and application.

Step 1

Development of the jointed PBPK model of DABE with its metabolites

A mechanistically jointed PBPK model of DABE composed of four PBPK sub-models of DABE and its metabolites (BIBR0951, BIBR1087, and DAB) connected in a sequential fashion, similar to what occurs in humans (Figure 2). The DAB model was adopted from a default model (SV-Dabigatran) in the software without any modification. The PBPK models of DABE, BIBR0951, and BIBR1087 were newly developed via a middle-out approach, based on both in vitro and observed clinical data.

For the BIBR0951 model, the intestinal disposition of BIBR0951 was described using an Advanced Dissolution, Absorption, and Metabolism (ADAM) model (Ezuruik et al., 2021). An effective permeability coefficient in humans (P_eff,man), which relates to the transport rate/permeability across the enterocytes, of BIBR0951 was predicted based on the polar surface area (PSA) and hydrogen bond donor (HBD) at 137 and 3, respectively, (PubChem, 2023), and used in the ADAM model to predict drug absorption. Based on the observed profiles, a basolateral (BL) global scalar governing the passing of the gut formed BIBR0951 into the hepatic portal vein was optimized to 0.1 unit, (Delavenne et al., 2013; Prueksaritanont et al., 2017). A whole-body distribution PBPK model was applied for prediction of steady-state volume of distribution (V_d,ss) (Rodgers and Rowland, 2007). The metabolic conversion of BIBR0951 to DAB was governed by the kinetics of two hydrolytic pathways: 1) CES1-mediated hydrolysis using the kinetic parameters of BIBR0951 hydrolysis in HLM (Udomnilobol et al., 2023); and 2) plasma esterase-mediated hydrolysis using an in vitro half-life of BIBR0951 disappearance in the human plasma. According to previous in vitro findings in which gut and hepatic CYP3A significantly metabolized BIBR0951 to several oxidative metabolites (Udomnilobol et al., 2023), the enzyme kinetic parameters for CYP3A4/5-mediated BIBR0951 metabolism were incorporated into the model, as the parallel pathways competing with the DAB formation pathway. Intersystem extrapolation factors (ISEF) for CYP3A4/5 enzymes were calculated based on the in vitro metabolism data in HLM and recombinant enzymes (Udomnilobol et al., 2023), and used in the extrapolation of in vitro to in vivo CL parameters (Supplementary Methods).

For the BIBR1087 model, the V_d,ss value was predicted from the whole-body distribution model (Rodgers and Rowland, 2007). The hydrolysis of BIBR1087 to form DAB was mediated by the kinetics of CES1 and CES2 enzymes.

For the DABE model, all physicochemical properties were obtained from the published literature (FDA, 2010b; Doki et al., 2019). An intestinal absorption of DABE was described by the ADAM model using an in vitro passive apparent permeability coefficient in Caco-2 cells (P_app,caco-2) as an input parameter (Ishiguro et al., 2014). Distribution of DABE was described by the whole-body distribution model (Rodgers and Rowland, 2007). Enzyme kinetic parameters of CES1- and CES2-mediated DABE hydrolysis were obtained from the previous literature (Udomnilobol et al., 2023), and incorporated into the model for driving the formation of 2 primary metabolites (BIBR0951 and BIBR1087). In addition, the in vitro half-life of plasma esterase-mediated hydrolysis was also considered for the conversion of DABE to BIBR1087. As shown previously that intestinal CYP3A significantly metabolized DABE to several oxidative products (Udomnilobol et al., 2023), the kinetic parameters of CYP3A4/5-mediated DABE metabolism was added into the model as a parallel pathway competing with CES1/2 hydrolysis. The ISEF values for CYP3A4/5 enzymes were calculated based on the in vitro metabolism data in HIM and recombinant systems (Udomnilobol et al., 2023). To model transporter-mediated DABE efflux, the P-gp K_m was fixed at 2.6 µM according to previously reported values (Yamazaki et al., 2019; Lang et al., 2021). The maximum transport rate (J_max) of P-gp was then top-down estimated from the observed data (Delavenne et al., 2013; Prueksaritanont et al., 2017).

Step 2

Top-down parameter optimization to capture DAB profiles

To capture the plasma DAB profiles following microdose DABE administration, a global sensitivity analysis (GSA) using Morris’s method was utilized for screening of additional input parameters highly influencing the plasma profiles of DAB. The intrinsic clearance for biliary excretion of BIBR0951 was top-down estimated by fitting DAB PK profile following microdose administration of DABE (Prueksaritanont et al., 2017), using Simcyp-plugin parameter estimation (PE) function to obtain the best model-fitted intrinsic clearance value.

To capture the plasma DAB profiles following DABE administration at therapeutic dose, the input parameters related to solubility, precipitation, and formulation were added into the DABE model. Either solution with precipitation or solid immediate release (IR) dosage form was selected as the formulations of DABE, depending on the clinical study design (Supplementary Table S1). A diffusion layer model (DLM) was used for modeling dissolution and precipitation of DABE in the gut lumen (Farhan et al., 2021). An aqueous phase solubility of DABE at pH 7.4 was obtained from a published value of 0.003 mg/mL (FDA, 2010a). Parameters controlling the precipitation of DABE, including the critical supersaturation ratio (CSR) and precipitation rate constant (PRC), were top-down estimated from the observed plasma DAB profile following oral administration of 200 mg DABE solution (Blech et al., 2008).

Step 3

Establishment of the final PBPK models

The final PBPK models of DABE and its metabolites were re-evaluated using the DABE-CTC interaction dataset following microdose and therapeutic dose of DABE (Delavenne et al., 2013; Prueksaritanont et al., 2017). Supplementary Tables S2–S4 list the final input parameters for all drug models.

Step 4

PBPK model qualification

The developed models were qualified by simulating several clinical PK scenarios including the single-dose (SD) PK, multiple dose (MD) PK, and DDI with CYP3A/P-gp inhibitors (CTC; itraconazole, ITZ; rifampin, RF and verapamil, VP). The simulated results were then compared to the observed qualification datasets (Supplementary Table S1). The satisfactory qualification was justified based on the pre-specified criteria, described in the “Model evaluation” section.

For inhibitor PBPK models, the ITZ and its metabolite hydroxy-itraconazole (OH-ITZ) models were modified from Chen et al. (2019) by addition of inhibitory constants (K_i) of ITZ and OH-ITZ for intestinal P-gp of 2 and 5 μM, respectively (Lang et al., 2021). For other perpetrators, the default models in the software were adopted without modification. The final input parameters for all inhibitor PBPK models are listed in Supplementary Table S5. Based on the final parameters and the dosing regimens (SD or MD), all the inhibitors employed were considered dual inhibitors of CYP3A (reversible and/or mechanism-based) and P-gp, but not inhibitors of CES1/2.

A top-down estimation was carried out via a default Simcyp-plugin parameter estimation (PE) module. Weighted least square and Nelder-Mead methods were chosen as the objective function and minimization method, respectively.

A sensitivity analysis of multiple input parameters was performed by a global sensitivity analysis (GSA) using Morris’s method. The numbers of levels, trajectories, and replications were set as default values in the software. The higher µ* value indicates greater importance of that parameter to overall simulation outputs.

Evaluation of PBPK models at all steps was justified based on the following three pre-specified criteria: 1) visual predictive check (VPC), 2) acceptance criteria for PK parameters, and 3) acceptance criteria for DDI prediction. For VPC, all observed plasma concentrations of DAB at various time points should be within the 5%–95% confidence interval of the simulated results. For acceptance criteria of plasma PK parameters (C_max and AUC_0-inf), the alternative success criteria at a 99.998% confidence interval for PK parameters were calculated from the observed mean and percentage coefficient of variation (% CV) (Abduljalil et al., 2014). The prediction of PK parameters was considered successful when the observed PK parameters were within the lower and upper boundaries. For the datasets without standard deviation (s.d.) or % CV, the simulated PK parameters should be within twofolds of the observed data. For DDI prediction, the satisfactory prediction for DDI magnitudes was evaluated based on Guest’s DDI prediction criteria (Guest et al., 2011). The predicted C_max and AUC_0-inf ratios should be within the lower and upper limits of Guest’s acceptance range. The precision and bias for overall prediction were evaluated by calculating the geometric mean fold error (GMFE) using a formula: GMFE = 10   ∑ log s imulated observed number   of   scenarios . A model with good precision and no bias should have a GMFE within the range of 0.80–1.25.

To determine the input parameters for hydrolysis of DABE and its intermediate metabolites (BIBR0951 and BIBR1087) in human plasma, the compounds were separately incubated in the plasma for up to 2 h. Following incubation at 1 and 100 µM of DABE and BIBR0951, there was no obvious difference in parent disappearance or formation of hydrolytic products in the plasma between the 2 drug concentrations tested (Supplementary Figures S1A, B). The mean t_1/2 values of DABE and BIBR0951 in plasma were at 364 and 55 min, respectively. These results showed that DABE and BIBR0951 were hydrolyzed in human plasma in a concentration-independent manner. In contrast, BIBR1087 was stable up to 2 h in human plasma (Supplementary Figure S1C).

For modeling BIBR1087 elimination, the kinetic parameters of BIBR1087 hydrolysis were determined based on the metabolite formation in HIM and HLM. The CES-mediated BIBR1087 hydrolysis was described by the Michaelis-Menten kinetics, with the V_max and K_m values shown in Supplementary Figure S2.

As shown in Supplementary Table S3, the fractions unbound in plasma (f_u,p) of BIBR0951 and BIBR1087 were at 0.227 and 0.018, respectively. The blood/plasma partition coefficient (B:P) of both compounds was approximately at 0.6 unit.

In Step 1 of model development, the developed PBPK model of DABE linking DABE with all 3 metabolites was first evaluated by simulating the DDI between DABE and CTC. As shown in Supplementary Table S6, the model well described the magnitudes of changes in DAB plasma exposure in the DABE-CTC interaction studies following both microdose and therapeutic dose of DABE. However, in both scenarios, the plasma concentration-time profiles and PK parameters (C_max and AUC_0-inf) of DAB were overpredicted (Supplementary Table S6; Supplementary Figure S3), suggesting the need for parameter optimization to capture the DAB profiles while maintaining the magnitudes of the DABE-CTC interaction.

In Step 2 of model development (Figure 1), the GSA using Morris’s method was applied to explore an additional mechanism to resemble the plasma DAB profile following microdose DABE administration (training dataset 1). The results suggested that the parameters related to BIBR0951 clearance (e.g., intrinsic clearance for biliary excretion, CL_{int, bile}) had a remarkable influence on the DAB exposure (Supplementary Figure S4). After incorporating the estimated CL_{int, bile} of BIBR0951 at 273 μL/min/million cells, the models adequately captured the PK profile of DAB following microdose DABE (Supplementary Figure S5A, left and middle panels), but still overpredicted the DAB profile following DABE at the therapeutic dose (Supplementary Figure S5B, left and middle panels). To recover the DAB profile following therapeutic dose DABE, the input parameters related to formulation, aqueous solubility, and precipitation were included into the DABE model, as described in the “Material and methods” section. Based on model training set 3 (Blech et al., 2008), the CSR and PRC values were estimated at 17.9 unit and 2.88 h⁻¹, respectively. After the incorporation of solubility/precipitation model and formulation of DABE, the models could reasonably predict the plasma DAB profiles following DABE administration at both micro- and therapeutic doses (Supplementary Figures S5A, B, right panel).

In Step 3 of model development (Figure 1), the final PBPK model of DABE linking with its metabolites was then used for simulation of DDI between DABE and CTC again. As shown in Figure 3, the model well captured the observed plasma concentration-time profiles of DAB following microdose and therapeutic dose either in the presence or absence of CTC. Furthermore, the plasma PK parameters and DDI magnitudes were reasonably predicted within the acceptance ranges (Supplementary Table S7).

The developed PBPK model of DABE linked with its metabolites was then qualified against several independent datasets, including SD PK, MD PK, and DDI with CYP3A/P-gp inhibitors (Supplementary Table S1).

As shown in Supplementary Figure S6, the model well simulated the observed plasma concentration-time profiles of DAB following single and multiple dosing of DABE at both the microdose and therapeutic dose. In addition, the model could adequately predict the plasma PK parameters of DAB in various scenarios, except for qualification dataset 9 (Table 1). Noteworthy that the observed values in this dataset 9 are higher than the corresponding values in dataset 2 with the same 300 mg DABE dose, but in a comparable range to those reported in dataset 10–11 following 150 mg DABE dose. Overall, the GMFE values for predicting DAB C_max and AUC_0-inf were at 1.07 and 1.09, respectively, indicating successful predictions with negligible bias.

For DDI prediction, the developed model was able to reasonably simulate the plasma concentration-time profiles of DAB following coadministration of DABE and various CYP3A/P-gp inhibitors (Figures 4A–F). Notably, there have been no reports on the plasma DAB profiles of unconjugated DAB for the data set 9–11 (Figures 4D–F); thus, we could not evaluate the simulated profiles using the VPC method. In addition, the DDI magnitudes (C_max and AUC_0-inf ratio) following DABE administration at the microdose and therapeutic dose were all well predicted within the acceptance ranges (Table 2). The GMFE values for predicting C_max and AUC_0-inf ratio were at 0.88 and 0.87, respectively, suggesting that the model can accurately predict DDI magnitudes with minimal bias.

Taken together, the results suggested that the comprehensively developed PBPK model of DABE with its metabolites adequately predicted the plasma PK profiles and parameters of DAB after oral administration of DABE with and without CYP3A/P-gp inhibitors. Therefore, the model was sufficiently qualified for further applications.

The qualified PBPK model of DABE with a link to its metabolites was subsequently applied to gain a mechanistic understanding of the DDI magnitude disparity of the DABE-CTC interaction observed following administration of DABE at the microdose (375 µg) versus therapeutic dose (300 mg). The simulated results of the DABE-CTC interaction (Figure 3) were dissected for the inhibitory effects of CTC on the intestinal P-gp and gut/hepatic CYP3A. When comparing the DABE-CTC interaction using the microdose and therapeutic dose of DABE, there was a remarkable difference in the F_g’/F_g ratio (∼3.3 vs. 1.5) of the gut metabolite BIBR0951 (Figure 5A). On the contrary, either small or no difference was observed for other presystemic PK parameters (F_a’/F_a and F_g’/F_g of DABE, and F_h’/F_h of BIBR0951) (Figure 5A). These results indicated that the disparity in the magnitudes of the DABE-CTC interaction observed between two doses of DABE was driven mainly by the intestinal CYP3A-mediated oxidative metabolism of BIBR0951, with a minimal contribution from the gut P-gp-mediated efflux of DABE or any other presystemic events. The higher magnitude of F_g’/F_g ratio at the micro-vs therapeutic dose is consistent with the in vitro report that the CYP3A-mediated oxidation was saturable (Udomnilobol et al., 2023). Interestingly, following the therapeutic dose of DABE with CTC, the intestinal CYP3A-mediated oxidation of BIBR0951 also played a significant, albeit relatively small role, with the F_g’/F_g ratio for BIBR0951 being higher than any other ratios shown in Figure 5A.

The qualified model was also used to dissect the inhibitory effects of VP on the intestinal P-gp and gut/hepatic CYP3A for the clinical DABE-VP interaction results (Table 2 dataset 10–11) following administration of DABE at the therapeutic dose (150 mg). Based on the inhibitor model input parameters, VP is much more potent than CTC as a P-gp inhibitor (Supplementary Table S5), and thus could serve as a better tool for assessing the potential role of the gut P-gp in the intestinal absorption of DABE. Consistent with this, VP appeared to have a slightly higher impact than CTC on the gut P-gp, with F_a’/F_a ratio of 1.2 versus 1.05 for CTC (Figure 5B). However, the relative impact of the gut P-gp-mediated efflux of DABE was still rather limited, when compared with the intestinal CYP3A-mediated oxidation of BIBR0951 (F_g’/F_g ratio of 1.5) (Figure 5B).

We first simulated the plasma DAB profiles over a wide dose range of 0.1–400 mg following single oral administration of DABE. When compared with the observed data, the model could reasonably predict the systemic exposure of DAB throughout the dose range employed. Based on the simulation outputs, the AUC_0-inf values of DAB were approximately dose-proportional from 75 to 400 mg DABE. When the dose of DABE was <75 mg (lower than therapeutic doses), the plasma DAB exposure was slightly over 2-fold lower than that extrapolated from the therapeutic dose range of DABE (Figure 6A).

To gain an insight into the underlying mechanisms of nonlinearity, the presystemic PK parameters were computed from the simulation outputs. By increasing the doses of DABE from 0.3 to 300 mg, there was a dramatic decrease in the F_a of parent DABE from approximately 1.0 to 0.4 (Figure 6B), consistent with its poor solubility. The decreases in F_a with increasing doses were not supportive of the gut P-gp-mediated DABE efflux being a major contributing factor (versus the solubility factor) in the overall intestinal absorption of DABE across dose levels. On the other hand, the F_g of BIBR0951 was increased following the ascending doses of DABE; this was in line with the in vitro result of the saturable CYP3A-mediated metabolism of BIBR0951 being the sole metabolic pathway mediating the gut metabolism of BIBR0951 (Udomnilubol et al., 2023). Noteworthy that the minimal changes in the F_g of DABE or F_h of BIBR0951 across different doses of DABE also are not inconsistent with the earlier in vitro findings that CYP3A contributed only partly to the metabolism at the gut and liver of DABE and BIBR0951, respectively (Udomnilubol et al., 2023). Taken together, the simulated results showed that both the solubility-limited DABE absorption and the saturation of intestinal CYP3A-mediated BIBR0951 oxidation were the two major contributing factors, opposing to each other, to the somewhat modest nonlinearity observed in the plasma DAB exposure.