PK data from the clinical phase I study (clinical trial approval of TDI01 was obtained from Beijing Tide Pharmaceutical Co., Ltd. Study participants received single-dose or multiple-dose oral administration of TDI01. The single-dose regimens included: 400 mg QD, 800 mg QD, and 1200 mg QD; the multiple-dose regimens included: 200 mg QD, 7 days and 400 mg QD, 7 days. Blood samples in single dose groups were collected prior to drug administration (0 h) and at 1, 2, 3, 3.5, 4, 5, 6, 8, 10, 12, 24, 36, 48, and 72 h post administration. Blood samples in multiple dose groups were collected at 0, 1, 2, 3, 3.5, 4, 5, 6, 8, 10, 12, 24, 48, 72, 96, 120, 144, 145, 146, 147, 147.5, 148, 149, 150, 152, 154, 156, 168, 180, 192, and 216 h.

Conventional compartmental PK models were tested sequentially. Model selection was based on changes in the Akaike information criterion (AIC), the precision of parameter estimates (relative standard errors, RSE) and standard goodness-of-fit (GoF) plots (Lindbom, Pihlgren and Jonsson, 2005).

One-compartment, two-compartment and hepato-enteral circulation models were investigated. After comparing the AIC and the model-fitting degree, optimal model was utilized as the structural model. The covariates known with significant impact were also included in the structural model firstly. The First-order conditional estimation with interaction (FOCEI) was used to fit the parameters in this analysis.

The random effects model included inter-individual random effects and residual random effects. An exponential model (Equation 1) was used to describe inter-individual variation, shown as follows: P i = θ ⋅ exp η i (1) where P_i is the PK parameter of each individual, θ is the PK parameter of the population and η_i represent the inter individual variation which follows a logarithmic normal distribution. A combined error model (Equation 2) was used to describe residual error: Y i j = C i j ⋅ 1 + ε 1 i j + ε 2 i j , (2) where, Y i j represents the observed values, Y i j represents the population predicted values and, ε₁ and ε₂ represent additive and proportional residual error, respectively.

The final PK model was evaluated graphically on GoF plots, including observed values versus individual prediction or population prediction, conditional weighted residuals (CWRES) versus TIME, absolute individual weighted residuals (|IWRES|) versus individual predictions, and the normality test of CWRES.

Bootstrap was performed to internally validate the final model. The original dataset was used to simulate 1,000 additional datasets, each of which was used to re-estimate the parameters with the final model. Median values and 95% confidence intervals (CI) were calculated and compared with the final model parameter estimates to assess the robustness of the final model.

Visual predictive checks (VPCs) were used to evaluate the predictive power of the final model. One thousand simulations were implemented and the 2.5th, 50th and 97.5th percentiles of the observed versus simulated data were compared.

Analyte concentrations that were below the lower limit of quantification (LLOQ) were reported as data below the quantification limit (BQL). If the proportion of BQL data was larger than 10%, the likelihood-based approach (such as the M1 method) along with Laplacian estimation was used (Ahn et al., 2008; Beal, 2001). Observations corresponded to absolute conditional weighted residuals of more than 5 (|CWRES| > 5) were regarded as outliers and omitted from the model-building process. These omitted points were re-introduced into the dataset and model evaluation at a later stage (sensitivity analysis) was performed to assess their impact on parameter estimates.

The population PK model was performed using NONMEM (version 7.5, ICON Development Solutions, Ellicott City, MD, United States). First order conditional estimation with interaction (FOCEI) method was used for parameter estimation (Keizer et al., 2013). R (version 3.5.3) was used for data preparation, graphical analysis, model diagnostics, and statistical summaries. Xpose^® (Jonsson and Karlsson, 1999) and Pearl Speaks NONMEM (PsN^®, version 4.8.0) (Lindbom, et al., 2005; Lindbom et al., 2004) were used for model diagnostics and to facilitate other model development tasks.

In this study, the proportion of BQL data in the overall data is 1.3% (10/786). A total of 776 blood drug concentrations from 39 subjects was included in the data set. As shown in Figure 1, the double peak phenomenon observed after single dosing indicated hepato-enteral circulation model should be included in the final model, and the starting and ending time of hepato-enteral circulation was set as 8–10 h.

The PK properties of TID01 were described by a one-compartment model with dosage effect and hepato-enteral circulation (Figure 2), which was described by the equations (Equations 3–5) as follows: K 20 = CL / Vc DADT 1 = FLAG * K G 1 * A 2 − K a * A (3) DADT 2 = K a * A − K 2 G * A 1 − K 20 * A 1 (4) DADT 3 = K 2 G * A 1 − FLAG * K G 1 * A 2 (5) where K20 is the clearance rate constant of central compartment; K_a represents first-order absorption rate constant; K_2G represents rate constants for TDI01 from the central compartment to the gallbladder compartment; K_G1 represents rate constants for TDI01 from the gallbladder compartment to the absorption compartment. A, A (1) and A (2) represent absorption, central and gallbladder compartments, respectively. FLAG means the time switch of hepato-enteral circulation.

The established model reasonably fit the time courses of TDI01 in the circulatory system with clearance (CL/F) of 95 L/h, volume of distribution (V _c/F) of 1400 L and the absorption rate constant (k_a/F) of 0.345 h⁻¹, respectively (Table 1). As shown in the table, all the model parameters were precisely estimated with relative standard error below 20%. Bootstrap results suggested that the estimated values of the final model parameters were close to the median and within the 95% CI from the non-parametric bootstrap. The goodness-of-fit plots for the developed models are shown in Figure 3. The observed values versus the population and individual predicted values were closely distributed around the line of identity. The conditional weighted residuals were randomly and homogenously distributed around 0. One thousand simulations were implemented and the 2.5th, 50th, and 97.5th percentiles of the observed versus simulated data were compared, as shown in the VPCs of single-ascending doses and multiple-ascending doses (Figures 4, 5). The VPC results for the Pop-PK model indicated that the established models were able to adequately describe the observed plasma concentrations, where most observed plasma concentrations fell within the 95% prediction intervals.

The simulated concentration-time profiles of TID01 under the designed dose regimens are shown in Figure 6. The regimens are listed in Table 2. The simulated PK parameters are listed in Table 3, which indicated that the C_max and AUC values at steady state of 200 mg BID (tested dose regimen) were higher than those of 400 mg QD (dose regimen used in clinical trials), but the differences were less than 25%. As shown in Figure 6, the simulated concentrations of TID01 in blood under the dosages of 200 mg BID, 7 days and 400 mg QD, 7 days are less than the C_max of 1,391 ng/mL under the highest single dose (1,200 mg) in clinical trial. The simulation results demonstrated the drug accumulation caused by hepato-enteral circulation would not affect the safety of TID01 usage under the tested multiple dosage regimen.