Models analysis

The population parameters with interindividual variability in the PK model were Cl, Ka, MTT, rchr and tchr (variable parameters). The population PK parameters without interindividual variability were Vc, Vp and Cld (nonvariable parameters). This set of parameters agrees with those selected by Andreu et al.¹⁶, except for new circadian parameters. Another nonvariable parameter was the coefficient of the final accepted covariate model, which is presented below.

All population parameters in the PBPK model have interindividual variability and included Clih, Bmax, and Kd. The variables rchr and tchr were added for circadian modulation as in the PK model.

The interindividual error model for variable population parameters was exponential, assuming log-normal distribution for error (η), with a mean value θ, as shown in Eq. (⁷) for the Cl parameter.

We then defined the normal distribution for the η vector (dimension equal to the number of variable population parameters) using the zero mean and covariance matrix Ω. We employed the log-transformed whole blood TAC concentration (ln C_c in PK and ln C_b in PBPK) as the model output to be fitted. We tested the proportional and additive residual error to establish the model that best fits the clinical data.

Interoccasion variability (IOV) reflects random IPV that may be transferred to the residual variability term if it is not considered³⁸. However, the term IOV is associated with occasions defined by dosing periods, number of samples, or administration instants, which must be anticipated to be properly implemented. This is not the case of a model that may be applied as a customized knowledge engine in a real-time MIPD for different clinical environments and changing drug administration scheduling. Therefore, we decided not to implement IOV in the models structure.

The variables BW, body mass index, body surface area (BSA, calculated from BW and height using the Du Bois formula), sex, age, height and haemoglobin concentration were tested as potential covariates with an influence on PK model parameters, using a forward stepwise approach³⁹. The PBPK model included available covariates by means of equations describing mechanistic models, as described in the previous section.

We analysed the correlations during the initial exploration, and the correlation models were linear or exponential for continuous covariates. The linear model correlation is as follows:

whereas the exponential model is as follows:

in which θˆi is the mean population variable without covariate influence, θc is the correlation coefficient, x is the covariate and xr is the reference value for this covariate. The last value was taken as the median value. The value of θc was estimated during the fitting population process. Categorical covariates such as sex were analysed assigning discrete values to each categorical level of the covariate.

The best correlation models were selected to be evaluated during the forward stepwise approach. We accepted a covariate-population variable relationship if the reduction in the minimum value of the objective function (MVOF) was somewhat less than 3.84 (significance level of 5% according to χ2 distribution for 1 degree of freedom in the nested model). In addition, an approximate reduction of 10% in intraindividual variability (coefficient of variation [CV]) was required in the associated population parameter for minimum clinical significance.

We performed the population model fitting with the blood TAC concentrations measured during the first monitoring day in the hospital (Prograf administration), executed in 2 steps. During the first step, we also performed an initial exploration of the models to ensure their robustness and reliability.

We estimated the population parameters as a second step by maximising the likelihood function associated with the hierarchical models, using first-order conditional estimation methods (FOCE and FOCE-I)^40,41. The absolute simulation time at the first hospital admission was set at 48 h. Additional details can be found in the Model analysis – Population fitting section in the ^{First Supplementary Document}.

The initial evaluation of the models’ behaviour was based on individual residual errors, conditional weighted residual, ε-shrinkage, η-shrinkage⁴² and relative standard errors (RSEs) of estimated parameters.

We also performed a noncompartmental analysis (NCA) based on TAC blood concentration measurements to compute the NCA AUC24 (log-trapezoidal method), differences model AUC24 minus the NCA AUC24 (Δ_AUC) and the root mean square error (RMSE)⁴³ of the blood TAC concentrations minus the model TAC concentrations (log transformed) RMSEc.

We performed goodness of model predictions by plotting the predictive simulated TAC concentrations against the measured blood TAC concentrations, showing the accuracy and characteristic dynamics of the computed TAC concentrations. We performed the prediction-corrected visual predictive check (pcVPC)⁴⁴ with 500 simulated replicas of the patients’ data set to confirm the proper adjustment of the data to the population distribution. Observed values and their 5%, 50% and 95% percentiles were presented against the 95% confidence/prediction interval corresponding to the same model-predicted percentiles. The pcVPC graphics were obtained with R scripting from Ron Keizer⁴⁵, using the replicas file computed using PhysPK software. Stratification with respect to the normalised dose was applied in the case of the pcVPC plot. Finally, the conditional weighted residual versus time after dose plot was presented as a diagnostic instrument to detect potential model misspecification⁴⁶.

We performed 600 Monte Carlo simulations with estimated population parameters for each TAC model and 3 TAC dose levels (0.025, 0.042, and 0.05 mg/kg, which refer to 1.5, 2.5 and 3.0 mg normalised for 60 kg of body weight) to calculate the mean AUC24 associated with each dose. Simulations were performed 36 h and the final 24 h were used, with the aim of ensuring a customised near steady-state for the patients. Covariates and input variables, such as BW, body mass index, sex, age and height, varied according to our paediatric population’s measurements. The bioavailability F (mean ± SD) and the time that the TAC concentration exceeded a gross value of 20 ng/ml as maximum tolerated concentration⁴⁷ were also added to the calculated exposures.

We adapted the models during the second phase (Advagraf data) to correct the prediction errors after the initial fitting with Prograf (day 1). Deviations in the model TAC concentrations could be induced by conversion to Advagraf (day 2) and intrapatient variability during the week. New TAC measurements were taken on day 8 (second monitoring day). The model was adapted using the following procedures (additional details can be seen in the Adapting Evaluation section of the ^{First Supplementary Document}):

Procedure 1: The model was first adapted (stage 1) to correct the prediction error of the blood TAC concentrations (C0) measured just before the morning administration of Advagraf on day 8. We employed the WLS adaptive technique using one model parameter (LibT for PBPK and MTT for PK). The parameter adjustment occurred at the instant of TAC formulation conversion (beginning of day 2). We performed a second parametric adaptation (stage 2) to maximise the prediction accuracy of the other TAC measurements during day 8. The parametric change occurred at the beginning of day 8. We applied adaptive Bayesian and WLS techniques in this second stage.

Procedure 2. We adapted the models to correct the prediction errors in the blood TAC concentrations, using all day 8 measurements simultaneously. The parametric adjustment occurred at the beginning of day 2. We applied adaptive Bayesian and WLS techniques.

The first adaptation of procedure 1 sought to adjust the model parameter that had a greater effect on the drug release rate. After the adjustment had been successfully applied and the first TAC measurement on day 8 had been accurately predicted, the need for a second model adaptation was evaluated to analyse whether there was significant intraindividual parameter variation. The Bayesian technique cannot be applied in procedure 1 - stage 1 because LibT and MTT are not in the set of the models’ population parameters, as defined in the Model analysis subsection (additional details appear in the ^{First Supplementary Document}, Adaptive modelling techniques).

We employed a log-trapezoidal integration method to calculate the reference AUC24 (NCA) on day 8, which was then compared with the AUC24 from the models. The accuracy of the TAC model predictions was quantified by means of the RMSE_c⁴³. We compared the results of these 2 procedures, along with the procedures’ effect on the models’ suitability as predictive engines for dosage recommendations.