Predictive engines based on pharmacokinetics modelling for tacrolimus personalized dosage in paediatric renal transplant patients

The development of predictive engines based on pharmacokinetic-physiological mathematical models for personalised dosage recommendations is an immature field. Nevertheless, these models are extensively applied during the design of new drugs. This study presents new advances in this subject, through a stable population of patients who underwent kidney transplantation and were prescribed tacrolimus. We developed 2 new population pharmacokinetic models based on a compartmental approach, with one following the physiologically based pharmacokinetic approach and both including circadian modulation of absorption and clearance variables. One of the major findings was an improved predictive capability for both models thanks to the consideration of circadian rhythms, both in estimating the population and in Bayesian individual customisation. This outcome confirms a plausible mechanism suggested by other authors to explain circadian patterns of tacrolimus concentrations. We also discovered significant intrapatient variability in tacrolimus levels a week after the conversion from a fast-release (Prograf) to a sustained-release formulation (Advagraf) using adaptive optimisation techniques, despite high adherence and controlled conditions. We calculated the intrapatient variability through parametric intrapatient variations, which provides a method for quantifying the mechanisms involved. We present a first application for the analysis of bioavailability changes in formulation conversion. The 2 pharmacokinetic models have demonstrated their capability as predictive engines for personalised dosage recommendations, although the physiologically based pharmacokinetic model showed better predictive behaviour.


Models development and adaptive techniques
Pharmacokinetics model Figure 1: PhysPK diagram of the PK model for TAC with a 2-compartment structural main system (CentralP, Barrier, PeriphP) and a 3-compartment transit and absorption (ColAbsCAT, CAT3) subsystem.
The pharmacokinetic (PK) model is described by a 2-compartment structural system and a 3-compartment transit and absorption subsystem (CAT3). Figure 1 shows a diagram of the model built with PhysPK. The structural system is composed of compartments CentralP (Central) and PeriphP (Peripheral), with a Barrier element controlling the TAC flow rate between them. The CAT3 system is shown at the bottom of the diagram and is connected with the Central compartment through a collector that sums the TAC absorptions from the 3 sections of the small intestine. Lines connecting cited elements in this diagram represent TAC mass flow rates The remaining lines that connect PK model elements with circular connecters, such as p3dCentralP (CentralP), are provided by the simulation tool to facilitate access to the variables inside. For instance, a computational integrator component can access the elimination TAC mass flow rate from the Central compartment and calculate the accumulated elimination TAC or access the TAC concentration and calculate the area under the curve (AUC). The diagram is provided to clarify the model's scope and is governed by the following equation system (1): The variables are defined in the main text of the article. Index j = 1, 2, 3 refers to the small intestine section number. Additional variables used in equation (1) include TAC masses in the stomach, small intestine section j, and colon: M s , M si,j and M c , as well as the TAC transit flow rate from the stomach, φ t,0 and from small intestine section j, φ t,j .
Gastrointestinal fluid volumes (CAT volumes) are defined as V cat,j = V cat f j , in which V cat is the total CAT volume, and f j is the volume fraction of the j zone, which were calculated as V cat = 0.475 l, and f s = 0.84 (stomach), f si = 0.13 (small intestine) and f c = 1 − f s − f si (colon), according to (Schiller et al., 2005) . Absorption from the small intestine to the central compartment has first-order kinetics, with rate constant k ta = K a /V si , where K a (ml/h) is the absorption constant and V si (ml) is the small intestine fluid volume.
Linear assumptions in central clearance and gastrointestinal transit were compared with saturable (Michaelis-Menten) assumptions using the parametricfitting methodology stated in the Methods section, with no improvements in TAC predictions or interindividual parameter variability.
The absolute sensitivity of TAC concentrations with respect to changes in CAT volume fractions f j and in the small intestine M T T was much less than 1, confirming the robustness of the model with respect to those parameters. This robustness is a desired feature in the PK model, given the high variability of gastrointestinal fluid volumes and mean transit times (Kimura and Higaki, 2002).
Circadian transitional zones in f w (t) ensure a smooth change between diurnal and nocturnal phases. The durations of the diurnal and nocturnal phases were assumed to be equal, and transitions were estimated to be 20% of the phase durations. A sensitivity analysis of TAC concentrations with respect to changes in duration parameters of the circadian waveform function confirmed the robustness of the model (absolute sensitivity values < 1). Model parametric fitting was also performed with various transition durations with no improvements in TAC predictions either interindividual variability The drug administration subsystems of the PK model include the drug form delivery time F DT (h), which is the time required for TAC to be fully available for transit. This was calculated at 0.083 h (5 min) in a fast-release formulation such as Prograf. Absolute values for TAC sensitivity with respect to F DT were much less than 1, confirming the robustness of the model in terms of F DT perturbations The time-dependent bioavailability F (t) is equal to the F value calculated with the standard equation F · D = AUC · Cl when Cl is constant. On the other hand, the mean value of F (t) converges to the standard bioavailability F for sufficiently long-term values (several administration periods) (Prado-Velasco, 2018). Additional model data. The oil -water partition coefficient (log) for TAC was assumed to be 3.3 (DrugBank, 2018). The partition coefficient of other tissues was considered equal to that of muscle tissue due to the nearly identical relationship observed between muscle and other body tissues (Poulin and Theil, 2000).

Physiologically based pharmacokinetic model
Tissue volumes are proportional to body weight (BW) based on ratio values obtained from published physiological databases (ILSI, 1994), except for the fat volume ratio P f (l/kg of BW), which was calculated though a correlation for children and extracted from Table 3 in (Deurenberg et al., 2007) to give: in which age is given in years, sex is 1 (males) or 0 (females), standard body mass index (BMI) is in kg/m2 units and BW is in kg. Fat volume is then V f = P f ·BW·10 3 (ml). Paediatric scaling for this reduced physiologically based pharmacokinetic (PBPK) model was completed by considering changes in water and lipid fractions in fat tissue according to Price et al. (2003). Exploratory simulations have shown that paediatric fat scaling is required to achieve optimal TAC concentration predictions, whereas scaling of other model tissues is not needed. This conclusion agrees with the largest changes in total volume, water, and lipid fractions of adipose tissue during ontogeny (Price et al., 2003). No other ontogeny scaling was implemented, given the objective of this study was to assess a reduced PBPK model as a predictive engine for a real-time model-informed precision dosing. Variables associated with gut absorption, liver metabolism and red blood cell binding (which are related to proteins with a significant influence on the accuracy of TAC distribution predictions) have been defined as population parameters to be fitted. Table 1: Physiological values for tissues in the PBPK model extracted from (ILSI, 1994;Pelekis et al., 1995;Poulin and Theil, 2000). See main text for nomenclature and units.
Tissue  Table 1 shows physiological data related to PBPK tissues. P i (ml/kg) provides the volume of tissue i through the equation V i = P i · BW (ml), in which BW is in kg. Index i takes values l (liver), k (kidney), g (gut), f (fat) or o (others). P f does not appear since it is found in equation (2). Gut lumen volume was considered equal to gut tissue volume. Others tissue volume was V o = BW · 10 3 − i V i (ml), in which the sum is extended by liver, kidney, gut tissue, fat and blood. The value P o is equal to V o /BW. The mean body mass density was considered 1 g/L as is usual in PBPK (ILSI, 1994). Heart volume V h = 150/70 · BW ml refers only to the heart chamber. Plasma volume (V p ) and red blood cell (RBCs) volume (V rbc ) are calculated from total blood volume through haematocrit (htc) (pu) as V rbc = htc · V b and V p = V b − V rbc . Water, neutral lipids and phospholipid fractions for each tissue are W f , nL f and phL f , respectively. These fractions are not defined for blood, given this tissue was separated into plasma and RBCs, which are considered homogeneous The venous perfusion blood flow rates for various tissues were proportional to BW, according to the equation Q i = q i · BW (ml/min), for the values q i as shown in Table 1. TAC binding is assumed to be a reversible quasi-static process with a single family of binding sites (Toutain and Bousquet-Melou, 2002;Copeland, 2000), according to the following relationship between bound concentration C b and unbound concentration C u : In short, saturable reversible binding is assumed according to equation (3) for RBCs; linear binding for plasma is f u p = C up /C p , in which C p is the total plasma concentration and C up is the unbound plasma concentration. Linear binding for the non-homogenous tissue phase of the remaining organs is f u t = C ut /C t , in which C ut is the unbound tissue concentration and C t is the total tissue concentration. The value of f u t is calculated according to the Langmuir approximation given by equation 4.21 of (Copeland, 2000), as follows: in which C prt is the tissue binding protein concentration and P T P = C prt /C prp is the ratio of binding tissue / plasma protein concentrations. Given that f u p verifies the first equality of equation (4) for plasma binding protein concentration C prp , it follows that: from which the term nK −1 d C prp can be solved and used in equation (4) to obtain f u t , once the binding tissue / plasma protein ratio P T P is known.
The main binding macromolecules present in plasma and tissue (albumin, globulins, lipoproteins) maintain a tissue plasma ratio from 0.3 to 1 in mammalian tissues, with 0.5 the most frequently reported intermediate value (Poulin and Theil, 2000). Therefore, 0.5 is commonly used as the default value for P T P in PBPK modelling software tools. We used 0.5 for P T P after verifying that the absolute sensitivity of TAC blood concentration to P T P is less than 1, and TAC predictions and interindividual variability were not affected when P T P was perturbed around the value 0.5. Liver elimination is considered linear and is provided by the intrinsic hepatic clearance per liver volume Cl ih (1/h), which is equal to Cl max ih · f h , in which f h is the hepatic TAC metabolic activity given by: The variable f a is the normalised activity of the hepatic referred cytochrome P450 (CYP), f m is the fraction of hepatic clearance due to the pathway (CYP3A4/5 or others) and IR3A4 is the inhibition ratio of CYP3A4 activity due to other compounds. The usual values of fractional pathways were taken as f m 3A4 = 0.25, f m 3A5 = 0.40 and f m oth = 0.35 (Gérard et al., 2014), whereas normalised activities were defined as 1 for the kidney transplant recipients in our study. Finally, IR3A4 was set at 0 according to the full pharmacological prescription of our study patients.
Assuming that LibT depends mainly on the drug release time, its value for Prograf was calculated using in vitro (simulated gastric fluid, pH 1.2) data for immediate release TAC formulation given by Shin et al. (2018). These data indicate that around 90% (0.9 pu) of the TAC mass is released 1.5 hours after administration . Solving LibT for the zero-order liberation equation (the last equation of Gut system in the Supplementary Material) results in: Therefore LibT = 100 min. We have verified proper TAC predictions and model robustness during the exploration phase for Prograf with LibT = 100 min.
The value of K elib for LibT = 100 min was defined to provide an average F ≈ 0.2 per unit (20%), in agreement with the reported mean value (Drug-Bank, 2018). This calculation was performed based on a parametric analysis of the PBPK model, using parameters fitted in the initial exploration phase and published values for the rest of model parameters. We obtained K elib (LibT = 100) = 13.5 1/h. The hematocrit value htc (pu) was calculated based on its relationship with measured haemoglobin concentration Chb (g/dl) (Worrall et al., 2016), equation (8).
Mathematical model. The model diagram in Figure 2 complements the standard flow diagram of the main text and is built and interpreted in the same manner as the diagram in Figure 1 for the PK model. The detailed scope of each tissue and the entire model are described below by their equations. The main variables and units are defined in the main text.
Heart tissue. This is a flow-limited tissue governed by the following equation (9): Heart volume flow rate Q h is the total cardiac output (CO).
Blood system. We modelled Blood using a membrane-limited tissue approach, with plasma and RBC compartments. The approximate distribution of TAC is 95% in RBCs and 5% in plasma, with linear and strong binding in plasma (f u p ≈ 0.012) (Nagase et al., 1994). Binding in RBC is extensive and saturable with maximum binding capacity B max and dissociation constant K d (see main text) (Jusko et al., 1995). The TAC distribution between RBC and plasma can be written by means of the following equation: in which C rbc is the RBC TAC concentration. The transfer of TAC across the RBC barrier is a complex issue, and its active mechanisms are not completely understood. However, a discussion of these mechanisms exceeds the objectives of this model (Tron et al., 2019;W et al., 1993)3. The mechanisms result in lower TAC transport from RBCs to plasma, which is addressed by non-symmetric mass flow rate transport RBCplasma coefficients, which account for the TAC distribution data (Nagase et al., 1994).
K pRBC and K RBCp (ml/min) denotes the mass flow rate transport clearances from plasma to RBC and vice versa. If the net RBC -plasma TAC flow rate is zero (RBC-plasma TAC fast-equilibrium) then K pRBC C up = K RBCp C urbc . The concentrations refer to unbound values in plasma and RBCs. A gross estimate of the mass transport K pRBC and K RBCp relationship can be obtained approximating the previous equation with the following equation: in which hematocrit has been calculated as the mean value 0.45 pu (45 %). Fast-equilibrium RBC-plasma is achieved through a sufficiently high value of K RBCp . We set K RBCp = 10 4 ml/min, resulting in K pRBC = 24.4 · 10 4 ml/min. Preliminary fitting of population PBPK model variables and a sensitivity analysis confirmed the accuracy of the TAC predictions and the robustness of the model to K RBCp and K pRBC perturbations.
TAC concentration in plasma is provided by the following equations: in which i takes values f (fat), l (liver), k (kidney) and o (others).
TAC concentration in RBCs is given by the following equations: where C brbc is the RBC TAC bound concentration. Whole blood TAC concentration is calculated with the following equation: Measurements of TAC concentrations in our clinical study are referred to whole blood; therefore, C b is the model output used in the parametric fitting.
Fat tissue. This is a flow-limited tissue with a partition ratio R f and linear binding. TAC tissue concentration is provided by the following equations: Variables C uf and C uvf are the unbound fat tissue and unbound venous fat tissue concentrations. Variables C f and C vf are the total fat tissue and total venous fat tissue concentration. Equations (15) gives unbound and total concentrations in fat and venous fat phases over time.
We should note the use of the tissue-venous partition ratio R t , which relates to phase unbound concentrations, unlike the partition coefficient R tt used in other studies, which relates to total concentrations. R tt and R t are related through the following equation: in which f u p and f u t are the unbound plasma fraction and unbound tissue fraction, respectively, described and calculated in the main text.
Liver tissue. Liver tissue is a flow-limited tissue with a partition ratio R l , linear binding and hepatic clearance, using the following equations: Concentrations follow the same designation as fat (see main text), replacing subindex f with l (subindex g is used for the gut venous blood flow rate). Clearance is chronomodulated as described in the main text.
Kidney tissue. Kidney tissue is a flow-limited tissue with a partition ratio R k and linear binding. TAC tissue concentration is provided by the following equations: Variables follow the same designation used in previous tissues.
Others tissue. This tissue combines the remaining body tissues not explicitly considered. It is considered a flow-limited tissue with a partition ratio R o (see main text) and linear binding. TAC concentrations are solved using the following equations: Gut system. This system consists of the gut lumen, where TAC is administered orally and absorbed by gut tissue through a one-order absorption chronomodulated process. The system is described by the following equations: TAC mass in Gut lumen and oral form, M gl and M f , are additional variables to those defined in the main text. Monte Carlo exploratory analysis of unbound plasma fraction. As described in the main text, the value of f u p (pu) was set to the mean value 0.012 pu (Zahir et al., 2001).
We performed a Monte Carlo study to verify the robustness of the model with a total of 100 simulations: 10 patients with parameters and dispersions obtained in a preliminary population fitting of the model and 10 f u p values normally distributed (mean = 0.012, standard deviation = 0.006). Fig. 3 shows that the TAC blood concentration profiles for the indicated percentiles maintain similar temporal dynamics. The trough TAC concentration C 0 varied from 3.3 to 8.2 ng/ml for the morning dose and from 4.0 to 9.9 ng/ml for the evening dose.

Adaptive modelling techniques
This section details all the boundary equations used with the adaptive Bayesian and weighted least squares (WLS) techniques presented in the main text. Decision variables in the Bayesian technique are based on the population variables selected in the population model fitting. Additional model variables are added in the adaptive WLS technique.

Adaptive Bayesian Optimisation.
Model adaptation was focused on correcting blood TAC prediction errors during the second monitoring day.
Boundary equations for adaptive Bayesian optimisation in the PK model are as follows: The value η 0 refers to the value of η (interindividual variation) at the fitting state prior to the adaptation. Incremental bounds and weights are defined to restrict the change in Bayesian model parameters according to the new information. The circadian modulation factor r chr and temporal phase t chr are strongly limited during the adaptation process. These boundaries were evaluated during preliminary explorations of the adaptive phase, and the results confirmed that final optimised variables remained inside the boundary range, with proper convergence of the numerical method. The same methodological strategy was used with the other boundary equations Boundary equations for adaptive Bayesian optimisation in the PBPK model are as follows: Adaptive WLS optimization. This technique was used in two patient scenarios. The weight of the objective function was set to 1 (W = 1).
• First scenario. Model adaptation attempts to correct only the trough concentration at the beginning of the second monitoring day.
• Second scenario. Model adaptation attempts to correct the prediction errors during the second monitoring day.
Only one decision variable was defined for the first scenario. This variable was M T T for the PK model and LibT for the PBPK model, with boundary equations as follows: Additional model parameters are added in the second scenario. In the PK model, decision variables and bounds are similar to the incremental Bayesian method, except for the addition of F DT to facilitate the drug formulation conversion. The variable t chr was removed here due to its negligible influence.
The low incremental limit for F DT ensures that this variable will be greater than 0.0833-0.08 hours (0.2 min), in which 0.0833 h (5 min) was the value set at Prograf fitting.
Decision variables and bounds for PBPK are defined similarly. That is, clearance and absorption parameters, together with drug liberation time, are allowed to change due to the new drug formulation. The negative incremental limit of K abs guarantees that the final value is greater than zero. Population fitting. Fig. 4 shows the sequence of stages followed by a general population parameter estimation study in PhysPK. All of them are optional, and their execution depends on the value of the associated flag (OTB_TwoStage, OTB_POP, OTB_postanalysis). Results of each stage have been saved in a XML file that can be read by the following step (automatic mode)

Model analysis
The first stage performs a preliminary analysis before executing the final population estimation stage, which is not based on the likelihood maximisation technique but on least squares minimisation. This old population parameter estimation method gives gross estimates for central parameters θ and population variances. The central parameter is related to the individual parameter φ through the interindividual dispersion model, which is given by the following equation: in which η is the interindividual dispersion. Taken advantage of the low CPU demand of this computational step, we used it for three main tasks: • Computing initial guesses for central parameters, theta, covariance matrix, Ω, and residual variances.
• Analysing potential relationships among individual φ and covariate candidates.
• Modifying and evaluating various population variable sets, their population variability and the model robustness for those selections.
We followed this technique during the initial exploration of our models. The first simulations were executed without covariate regression models. The first individual parameter fitting based on minimisation of the least-square residual error was performed initially for randomly selected patients. Plots with selected magnitudes gave us an initial idea of the model's general behaviour for the chosen parameter values.
Subsequent executions of the two-stage (TS) method gave us initial guesses for θ and variances, as well as suggestions for nonvariable population parameters.
Covariate candidates were analysed, seeking correlations between potential covariates and the values of each φ estimated for each patient. Selected covariate relationships were subsequently added to the preliminary TAC model. We then performed parametric analyses modifying covariate values and model correlation coefficients for the operational point defined by the previous estimated θ values. Correlation models that gave robust numerical values closest to the known range of AUC 24 were selected. We also evaluated other model output variables such as F to ensure that their range was correct.
Monte Carlo studies with a few simulations (100) with selected covariate correlations were executed to confirm the range of AUC 24 with 2.5% and 97.5% percentiles.
A final execution of the TS method with selected covariates and regression models was performed. Coefficients of correlation models were assigned various values with the goal of selecting one of them as an initial guess for the next stage of parameter estimation. Correlations between population variables might have been detected in this step (non-diagonal Ω). However, they were not found, in agreement with other authors (Andreu et al., 2015).
The initial exploration ended with the computation of the dynamic sensitivity of blood TAC concentration (log-transformed) with respect to the parameters of the model, to verify the stability and robustness of the model structure.
The normalised sensitivity of y(t) with respect to the φ(t) parameter is defined as follows: in which y(t) is taken as TAC concentration (log-transformed) at time t and φ(t) is the model's parameter with respect the sensitivity is computed. Model parameters can be time-dependent, as occurs with parameters that have circadian modulation. In general, the sensitivity S φ y (t) varies with time. Absolute values of sensitivity much larger than 1 can induce numerical stability problems. On the other hand, values very close to zero at any time instant t suggest that the measurements of y at t will have scarce influence on the fitting of φ. Typically, S φ y (t) grows just after a new drug administration, indicating the temporal range when measurements have more influence on the fitting of the φ value The second stage of Figure 4 executes the population model fitting using Likelihood maximisation with additive and proportional residual error models for the measured TAC concentration (log-transformed). The First Order Conditional Estimation (FOCE) (see main text) method was applied to the covariate forward stepwise approach, once the residual error model with better behaviour (lower IIV values and parameter standard errors) was selected. The FOCE-I method did not improve the goodness of estimation (parameter accuracy for either IIV values); thus, FOCE results were selected in our study.
Optimisation mathematical techniques related to likelihood maximisation in PhysPK include the possibility of organising parameters to be estimated in groups, in an iterative procedure controlled by the user. Grouping of parameters was based on differences in the sensitivity analyses and in their physiological or pharmacokinetic meaning. This technique helps prevent local extremes or numerical convergence problems, particularly with gradient-based optimisers. PK model parameters Cl, K a , M T , V c , V p and Cl d were assigned to group 1, r chr to group 2 and the coefficient of selected covariate correlation to group 3.
With respect to Ω, the variances ω 2 Cl , ω 2 Ka , and ω 2 M T were assigned to group 1, and ω 2 r chr to group 2. The value of ω 2 t chr was set to 0.01 according to exploration results. An approximate estimation explains this value. Considering Normal distribution and 99.7 % = 1 h, 3 ω t chr ≤ 1/2, and thus ω 2 t chr < 0.02. PBPK model parameters were all assigned to group 1, except for r chr , which was assigned to group 2. The parameter t chr was addressed as it was in the PK model. Concerning the Ω matrix, all variances were in group 1, except for ω 2 r chr , which was in group 2.
Once the parameter population distribution was known (mean and dispersion), the model was customised to each patient at the third stage of Figure 4, maximising its likelihood through the Bayesian theorem and the prior distribution (Kelman et al., 1982). As a result, the set of η values that gave the set of model parameters φ for each patient was obtained The initial values of the model state variables are customised for each patient in the 3 stages of Figure 4, using a PhysPK technique called conditioned simulation. This technique performs simulation loops based on the drug admin-istrations defined in the data set file for a particular patient, until a user-defined condition is fulfilled. Common system conditions are the steady state of the system or the achievement of a particular exposition in a compartment.
After various evaluations, the condition selected was the number of drug cycles (days) prior to the first admission to hospital, where TAC measurements were acquired. Although the patients followed the controlled Prograf administration 6 days prior to admission, it was verified that 2 days were sufficient to achieve the steady state. Therefore, only 2 simulation loops were selected as the condition to be fulfilled, to reduce the computing load. As a consequence, the absolute simulation time at hospital admission started at 48 hours. This time value may be considered a reference.
Adaptive procedures evaluation. Figure 5 presents the 2 procedures designed to evaluate the adaptive techniques and models of the study. Plots are based on the PBPK model and a randomly selected patient (number 5) was taken as an example. The fundamental details of the procedures are described in the main text The top left plot shows the evolution of TAC concentrations (log-transformed) from the beginning of day 1 (first admission day to hospital) to the end of day 8 (hospital re-admission), both according to the PBPK model simulation (orange, 1 to 8 days) and blood TAC measurements (violet, 1 and 8 days). The first sampling absolute time at day 1 occurs at time = 48 hours, as a common reference for all patients indicated in previous Section. Accordingly, the last measurement of blood TAC taken at the end of the hospital re-admission occurs at time = 240 h (48 + 8 · 24). This first plot shows that model predictions a week after TAC formulation conversion (beginning of day 2) from Prograf to Advagraf are inaccurate.
Procedure 1 (bottom block) is composed of a first adaptation using the adaptive WLS, to correct the predicted trough value just before the Advagraf administration at day 8 (first measurement of TAC marked in red colour). The WLS technique calculates the change in M T T in the PK model or LibT in the PBPK model, as described in the theoretical section of the main text. The bottom left plot shows the predicted corrected TAC concentration and the TAC measurement (C0, time = 220 h) after model adaptation. That parametric adjustment is applied at time = 48+24 = 72 h, that is, when the TAC formulation conversion occurs (beginning of day 2). The simulated TAC concentration as a function of time is presented from the beginning of day 2 until beginning of day 8 (72 to 220 h).
A second stage in procedure 1 applies adaptive WLS and Bayesian techniques after the first adaptation. Adaptive methods start from a model with an accurate prediction of TAC concentration at the beginning of the second monitoring day (day 8) and are based on TAC measurements acquired during that day. Any relevant change in the model's parameters at this second stage of procedure 1 suggests that the intrapatient variability is significant. The 2 bottom plots at the right show the corrected TAC concentration predictions versus the actual TAC measurements during the second monitoring day (24 h), for adaptive Bayesian and WLS techniques.
Procedure 2 applies adaptive Bayesian and WLS techniques to correct deviations in TAC predictions during day 8 in a unique stage. The parametric adjustment is applied at time = 72 h (beginning of day 2). This evaluation procedure is compared to the first procedure in the main text. The 2 upper plots at the right of the Figure 5 show the corrected TAC concentration predictions versus actual TAC measurements during the second monitoring day, for adaptive Bayesian and WLS techniques

Results
Fitting and evaluation of PK model (Prograf data)  Exploration phase. Fig. 6 shows the TAC concentrations (log-transformed) and the temporal variation of systemic clearance Cl due to the circadian rhythm (Cl chr (t)) for the randomly selected Patient 1. Plots start at the reference time t=48 h, according to the conditioned simulation described in the Methods section. The PK model was fitted for this patient by the minimisation of least-squared errors between central and blood TAC concentration (log-transformed). The values of the fitted parameters were Cl = 8.1 l/h, K a = 20.8 ml/h, M T T = 1.73 h, Cl d = 49.7 l/h, V c = 11.6 l, V p = 464.7 l, r chr = 0.17 and t chr = 0.49 h. TAC exposition (AUC) 24 hours after drug administration as computed by the PK model was AUC 24 = 205.9 ng/ml · h, whereas the log-trapezoidal value based on the measurements (NCA) was AUC 24 = 207.8 ng/ml·h. Bioavailability as given by the model was F = 18.0%. RMSE between central TAC concentrations and blood measurements (log-transformed) was RMSE c = 0.034.
The least-squared error minimisation was subsequently applied to the full population, by means of the TS technique, to obtain gross estimates of θ i and η-variances ω 2 i of the PK model parameters. Potential correlations between covariates and guess estimations of model parameters were evaluated to detect significant correlation coefficients according to linear and exponential regression models. Two relationships were selected at this phase, BSA (m 2 ) -Cl d (l/h) with the greatest correlation coefficient and Chb (g/dl) -V p (l) with a smaller correlation coefficient. Plots of covariates candidates can be found in the Appendix.
PK model parametric analyses were executed using the linear and exponential correlations, with BSA and Chb as variable parameters (first analysis) and with correlation coefficients as variable parameters (second analysis), using the population parameters estimated in the TS stage. The model confirmed the robustness with respect to correlations parameters, with AUC 24 within the population data range. Exponential correlation models for candidate covariates showed more reduced and proper (more similar to clinical data) values of AUC 24 .
By performing the Monte Carlo method with 100 simulations at this point, using the variances estimated by the TS technique, we confirmed that TAC concentrations and AUC 24 were consistent with the clinical population data, particularly with the exponential regression models.
We performed a second execution of the TS method with selected covariate regression models and initial values for θ i and Ω entries from the previous TS execution. Values for the exponential correlation coefficients were taken as 0.1, 0.3 and 0.5. Table 2 shows the results of the TS method for exponential coefficients equal to 0.1, because this value resulted in the lowest interindividual variance. The smallest estimates of variances occurred in the V c , V p and Cl d variables, which were selected as nonvariable population parameters. These parameters were used as initial guesses in the population fitting phase of the PK model with Prograf administration (see main text).
We completed the initial exploration by calculating the normalised sensitivities of the TAC concentrations (log-transformed) to PK model parameters, as shown in Figure reff:PKsens. The results confirmed the model's structural stability. Plot 7d also confirmed that 48 h were sufficient to achieve a near steady-state situation for TAC concentrations.
The parameters presented in Figure 7 are clustered into variable parameters (Fig. 7a), nonvariable parameters (Fig. 7b) and circadian parameters (Fig.  7c). The first group (θ Cl , θ Ka and θ M T T ) had the strongest influence on TAC concentrations, whereas the second group had less influence, particularly θ V c . The parameter θ t chr exerted a very small influence on TAC concentrations.
The sensitivities of TAC concentrations to the exponential coefficients of the covariate correlation models were also very low. Their behaviour was similar to that of θ r chr in Figure 7c.   Figure 8 presents model TAC predictions, fitted with the Prograf data from our stable paediatric kidney transplant recipients. TAC predictions and measurements were selected for 5 random patients. Two patients were extracted from the top plot and are presented in the bottom plots to clarify the behaviour of the TAC predictions. Patient 10 ( Fig.  8b) showed the most common TAC temporal pattern in our population. Table 3 shows a comparison with parameters obtained by Andreu et al. (2015) model. The M T T parameter from Andreu et al. has been converted through M T T /3 to agree with our MTT definition. We modified the K a from Andreu et al. to K a · V si = K a · 0.13 · 0.475 · 10 3 = K a · 61.7 to agree with our K a definition (see the mathematical development of the PK model). The parameters' values are written as value/BW (kg). The ratio value/BW is presented to facilitate the comparison between models because of differences between adult (Andreu et al.) and paediatric (our PK model) kidney transplant recipients. Comparisons between residual variabilities should consider the different type of error function (proportional in the Andreu et al. model and additive in our PK model). The proportional error variable multiplies the output y = ln(C) in the Andreu et al. model, which in our study has values ranging from 1.1 to 2.3 (blood TAC concentration C in ng/ml, second supplementary file online).    In terms of the normal distribution assumption for parameters and errors, Figure 9 presents the conditional weighted residuals (CWRES) -time for the fitted TAC PK model. The plot shows the time after diurnal and nocturnal Prograf doses. Although the mean CWRES values deviate slightly from zero at several sampling times, the diagnostic plot indicates that the PK model properly fits the simulated data.

Likelihood-based population fitting
We performed a Monte Carlo study of 600 simulations with the PK model for 3 TAC normalised doses to find a dose value as a reference for paediatric kidney transplant recipients. Table 4 shows the AUC 24 results. Each simulation was performed for 48 + 24 hours. Therefore, the time that TAC concentration surpasses MTC (20 ng/dl) (t MTC ) is referred to as 72 hours.  Initial exploration. Figure 10 shows a TAC concentration (log-transformed) and circadian modulation of intrinsic hepatic clearance Cl ih · V l (l/h) predicted with the PBPK model for Patient 1. The model was fitted by minimising the least-squared errors between simulated and blood measurements of TAC concentrations. Fitted parameters were Cl ih = 591.1 1/h, B max = 103.3µg/L, K d = 31.2µg/L, r chr = 0.28 per unit (pu) and t chr = 0.107 h. TAC exposition was AUC 24 = 204.0 ng/ml · h (PBPK model), and AUC 24 = 207.8 ng/ml · h (NCA log-trapezoidal). Bioavailability as computed by the model was F = 16.7 %. The RMSE between the model blood TAC concentrations and blood measurements (log-transformed) was RMSE c = 0.044 for this patient. We applied the least-squared error minimisation technique (TS), parametric analysis and exploratory Monte Carlo in the same manner as in the PK model. However, we did not find additional covariates to those used in the mechanistic equations (see Model development). The final execution of the TS technique for the PBPK model resulted in the estimated mean θ i parameters and gross estimated Ω entries shown in Table  5. The estimated covariances between the parameters were negligible, and thus Ω was set as a diagonal matrix. Values for the previous table were used as initial guesses in the population-fitting phase of the PBPK model with Prograf administration.
(a) Normalised sensitivities of TAC to θ CLih (blue), θ Bmax (red) and θ Kd (orange) with LibT = 100 min.  To finalise the exploration stage, we conducted a dynamic normalised sensitivity study with the population parameters defined by values found through the TS technique. Other model parameters were defined by their physiological or average representative values (e.g., body weight). Figure 11 shows the sensitivity of TAC concentrations (log-transformed) with respect to the mean population parameters of PBPK, clustered in variable parameters (Fig. 11a), circadian parameters (Fig. 11b) and LibT (Fig. 11c). The plots confirmed the structural stability of the PBPK model. Plot 11d showed that steady-state TAC concentrations are reached at approximately 48 h in the PBPK model TAC sensitivities are greater at time points after and near Prograf administrations. Therefore, measurements obtained at those time points have greater influence than others during the model fitting.
The parameter θ Clih has a greater influence on TAC blood concentrations than the RBC binding parameters. The parameter θ t chr scarcely affects TAC concentrations, as occurs in the PK model.
The LibT parameter was included in the sensitivity study, even though it was not defined as a population fitting variable. This variable is related to the TAC formulation and is therefore relevant for understanding the behaviour of the PBPK model when LibT changes. As shown, the average TAC absolute sensitivity to LibT increases with the value of LibT , in addition to presenting abrupt behaviour close to TAC administrations. However, the absolute values for sensitivity are moderated and do not affect the model's stability.
Plot 11d shows that maximum TAC concentrations (C max (t max )) decrease and that t max increases when LibT increases. This drug behaviour is expected when associated with a prolonged-release formulation.  Likelihood-based population fitting. Figure 12 shows the blood TAC concentration (log-transformed) predictions of the PBPK model, fitted with Prograf data from our stable paediatric kidney transplant recipients. We selected plots and measurements for 5 random patients.
Patients 10 and 15 were extracted from the top plot and are presented in the bottom plots to clarify the TAC concentration behaviour. Comparisons of PBPK TAC temporal concentration with PK TAC temporal concentrations ( Fig. 8b and Fig. 8c) show several slight differences in TAC patterns (dynamics). For example, the PBPK model presents faster changes in TAC slopes (derivative) than in the PK model. The TAC patterns are consistent with the circadian rhythms (Barbarino et al., 2013).
Population parameters estimated for our PBPK model can be compared with other published values. The mean intrinsic hepatic clearance (Cl ih · V l ) was 922.4 l/h for a BW of 55 kg. This intrinsic hepatic clearance is much lower than the 10600 l/h calculated by Gérard et al. (2014). However, Gerard et al. determined this parameter with only 6 adult liver transplant recipients and did not consider the circadian modulation effect. In fact, the predictions versus the observed TAC concentrations from the study by Gerard et al. are within the range of 0.33-fold to 3-fold in approximately 90% of the cases. These high TAC deviations could be explained by the fact that the authors' primary intention was to identify the most influential variability factors in TAC trough blood concentrations and not to achieve an accurate predictive model. The predicted intrinsic hepatic clearance from human liver microsomes was 2943.6 l/h, with a range of 1032.9-4851 l/h for a BW of 55 kg (Gertz et al., 2010), which is closer to our estimate.
The estimated mean RBC TAC binding values presented in the main text, B max = 11 µg/l and K d = 1.6µg/l, can be compared with the mean values for adult liver transplant recipients (B max = 418 µg/l and K d = 3.8 µg/l) (Jusko et al., 1995). The discrepancy could be explained by differences in our paediatric population, consistent with the high variability found for these binding parameters for various populations (Jusko et al., 1995).      Figure 14 shows the evolution of the PBPK model for 3 randomly selected patients after adjusting the LibT value (stage 1 of procedure 1) to correct the trough TAC concentrations at the beginning of the first hospital re-admission day. Table 7 shows the adjusted LibT values, with the columns presenting LibT values at Prograf fitting (day 1) and after executing the WLS adaptive technique (WLS C0). The second block of the table presents the C0 prediction deviation (∆ C0) before (Prograf) and after (WLS C0) the adaptation. LibT increased the value of the Advagraf formulation, with the exception of Patient 15. The reduction in LibT in that patient was related to the undervaluation of measured C0 by the model before adaptation, in contrast to the other two patients. Blood TAC AUC 24 was reduced from 155.9 to 110.0 ng/ml · h (Patient 1) and from 206.2 to 140.8 ng/ml · h (patient 20) and increased from 196.4 to 367.4 ng/ml · h for Patient 15, who was the only patient with a reduced LibT after PBPK model adaptation.
The predicted blood TAC concentrations versus blood TAC measurements after the Advagraf administration on day 8, after the PBPK model was readjusted in the second stage of procedure 1, are shown for the same 3 patients as Figure 15. The left plots were obtained with the Bayesian adaptive technique, and the right plot was obtained with the WLS adaptive technique. Both adaptive methods achieved an improvement in the objective function (likelihood in Bayesian and least squares in WLS) and thus in TAC predictions during day 8.
Predictions for TAC concentrations for Patient 20 show that the Bayesian adaptive method might have limitations in achieving accurate predictions through the adapted model in cases in which LibT requires a significant adjustment.
The values of the model parameters after adaptation appear in Table 8. C0instant refers to the beginning of re-admission to hospital. Values of parameters at this time point are equal to those set at the Prograf fitting, with the exception of LibT . The goodness of fit was better for the adaptive WLS than for the adaptive Bayesian technique, according to the RMSE c and AUC 24 values.
This result agrees with the TAC temporal patterns shown in Figure 15. The hepatic intrinsic clearance Cl ih and the absorption constant rate K abs appear to be the parameters that govern the readjustment, whereas those related to  RBC binding B max and K d , and circadian activity reduction r chr , are scarcely changed.
The r chr parameter is not presented in the entire patient table (see main text) because its variation was very close to zero with both adaptive techniques, and it could be extracted from the decision variables without affecting results.  Figure 16 shows the TAC concentration predictions for days 2 to 8 of the PBPK model adjusted with Bayesian and WLS adaptive methods (procedure 2), for randomly selected patients. Plots are centred on days 7 to 8 for Patient 20 to facilitate the inspection of TAC predictions versus measurements. The plots suggest that model TAC concentration predictions are more accurate with the WLS than the Bayesian adaptive method in Patients 5 and 20. Table 9 shows the model parameters at the Prograf fitting and after the adjustment of the PBPK model with adaptive Bayesian and WLS methods. Patients 5 and 20 showed increased LibT from 100 min to 437.4 min and from 100 min to 319.5 min, unlike the smaller increase in Patient 15 (100 min to 222.5 min). The value of the LibT parameter cannot be adjusted by the Bayesian method because it was not defined as a population parameter, which explains the poor performance of the Bayesian method compared with the WLS method. The changes in model parameters confirm the same trend detected in procedure 1 for these patients. The parameter r chr was removed from the final adaptive evaluation because their variation was kept near zero for both methods.
PK model: Adaptive procedure 2 (Advagraf data). Figure 17: Adapted PK model TAC concentrations (blue) (procedure2) versus TAC measurements at day 8 (red) (log-transformed) for randomly selected patients. Time in hours. Figure 17 shows the TAC predictions for days 2 to 8 of the PK model adjusted with Bayesian and WLS adaptive methods (procedure 2) for randomly selected patients. The plots have been centred on days 7 to 8 for Patient 20 to facilitate the inspection of TAC concentrations. The plots suggest that model TAC concentration predictions have similar accuracy for the 2 adaptive techniques. This result agrees with the low correlation between F DT and the release rate of the TAC formulation (see main text).
Model parameter values at the Prograf fitting and after the adjustment of the PK model with adaptive Bayesian and WLS methods are shown in Table  10. The addition of F DT to the WLS adaptive method explains the parameter  Bioavailability as a function of LibT. Figure 18 shows F (t) as a function of LibT , as calculated by the PBPK model starting from a washout condition. After the initial transitory F reaches a near stable value with a temporal average equal to 15% (LibT = 100 min), 9% (LibT = 200 min) and 4% (LibT = 400 min). This reduction in F is associated with prolonged drug release. However, understanding this relationship requires further research.

Appendix: Covariates data
Only one covariate (BSA) was accepted by our PK TAC model, despite an intuitive tendency to include BW and age in a paediatric population. However, BSA is strongly correlated with both the patients' BW and age. The relationship between BSA and BW is included in the well-known Du Bois formula. Figure  19 shows the association between BSA and age in our population. Figure 19: BSA (m 2 ) vs. age (years) in our paediatric population.
The BW-BSA association in our population is presented in Figure 20.