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

1 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports predictive engines based on pharmacokinetics modelling for tacrolimus personalized dosage in paediatric renal transplant patients Manuel prado-Velasco 1 ✉ , Alberto Borobia2 & Antonio carcas-Sansuan 2 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. The immunosuppressant tacrolimus (TAC or FK506) is a calcineurin inhibitor (CNI) employed to reduce the risk of acute rejection and allograft loss of many types of solid organs. The clinical application of TAC has increased in recent decades, unlike other CNIs such as cyclosporine A, due to the lesser efficacy of the latter1–3. TAC is highly lipophilic and barely water soluble, which explains the high variability of its pharmacokinetic properties when orally administered4. TAC is metabolised in liver microsomes mainly by CYP3A4 and CYP3A5 (in humans), to form a major metabolite with negligible pharmacological activity4. A second metabolite with equipotent activity to4 is scarcely produced. The first commercial oral formulation of TAC was a fast-release hard gelatine capsule (Prograf), with a recommended 12 hours between administrations. The bioavailability of Prograf is approximately 20% in adults5; however, this is largely affected by numerous factors, including food intake. Advagraf, a prolonged-release formulation, was subsequently developed to facilitate treatment adherence with the once-daily administration of TAC (usually in the morning). Selecting the optimal TAC dosage for patients who have undergone kidney transplantation (or other solid-organ transplantation) is difficult due to the narrow therapeutic index and large interpatient variability of TAC, the latter of which is associated with oral absorption variability, as well as the influence of CYP3A (cytochrome P450, family 3, subfamily A) and P-glycoprotein (PgP) enzyme activity and their polymorphisms1. TAC pharmacokinetics and toxicity are complex issues, and the need for individualised therapeutic drug 1Department of Graphics Engineering and Multiscale Modelling in Bioengineering Group, University of Seville, Seville, Spain. 2La Paz University Hospital, School of Medicine, Autonomous University of Madrid, Madrid, Spain. ✉e-mail: [email protected] open 2 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ monitoring (TDM) is well known1. Target TAC exposure recommendations maintain the trough blood TAC concentration (C0) as an alternative to blood TAC area under the curve (AUC), and various AUC ranges have been suggested as a function of the C0 target and formulation (twice daily for Prograf; once daily for Advagraf). The recommended C0 target for paediatric kidney transplant recipients at more than 2 months after transplantation is between 5 and 10 ng/ml, and the corresponding target AUC24 (AUC of blood TAC concentration during the last 24 h) is 180–350 ng/ml · h6. The new TDM TAC consensus report includes for the first time the monitoring of intrapatient variability (IPV) exposure as a biomarker to predict treatment outcomes in kidney transplantation, although further research is required to validate the thresholds and the timing for measurements6. TAC IPV is defined as the TAC blood concentrations fluctuations within a subject over a period of time for which the TAC dose is constant7. Several mathematical formulations can be applied to compute the TAC IPV, in which the TAC blood concentrations are corrected for the drug dosage, D, when D is not stable7,8. Borra et al. found that a high TAC IPV was a risk factor for poor long-term outcomes in kidney transplantation8. Subsequent studies have confirmed the prevalence of large TAC IPV in the clinical setting and its influence on graft failure for long-term kidney transplant recipients7,9,10. The main factors that cause a high TAC IPV are considered to be the analytical assay, food, diarrhoea, drug-drug interactions (DDIs), genetic factors, nonadherence and the use of generic TAC formulations7. Several pharmacokinetic TAC models have been developed to more accurately define dosages for target exposure, in line with the model-informed precision dosing (MIPD) approach11. The models include compartmental and noncompartmental approaches, as well as a physiological basis. We use the term “physiologically based pharmacokinetics” (PBPK) to refer both to individual and population-oriented models and employ the term PK to refer to the nonphysiological pharmacokinetic models. The PK models for predicting TAC (Prograf formulation) concentrations in patients (with applications for dose adjustment) include a 1-compartment system with an intravenous (IV) drug administration model for paediatric haematopoietic stem cell transplant recipients12; several 2-compartment system models with first-order oral absorption for adult kidney transplant recipients13, adult liver transplant recipients14 and paediatric kidney transplant recipients15; and a 2-compartment system with absorption governed by a 3-compartment absorption and transit (CAT) subsystem for adult kidney transplant recipients16. Fewer studies with PK models have evaluated TAC concentration predictions with Advagraf, although the results from Benkali et al.17 for adult kidney transplant recipients showed that extended-release TAC can be assessed using a chain of transit regions for absorption. Outcomes from these cited studies have confirmed that PK models offer a suitable approach to achieve personalised TAC dosages, although the authors did not consider the patients’ temporal progression and real-time MIPD. The PBPK approach has demonstrated better predictive capability than nonphysiologically based PK approaches18 and is therefore recommended when developing drug and clinical prescriptions for special populations, such as for paediatric patients19. PBPK models therefore seem more suitable than PK models for managing real-time MIPD11. However, most of the few studies about TAC distribution based on PBPK modelling are not designed for MIPD, but for other issues. These ones include the quantification of DDI mechanisms20 and the analysis of the sensitivity of TAC exposure to covariates21. This explains why those models had significant errors in TAC concentrations predictions, regardless of the use of an extensive whole-body approach20 or a simplified approach (reduced number of tissues)21. In addition, most of the published PBPK models were not designed to assess the patients’ progression. To the best of our knowledge, there are no PK/PBPK models for TAC distribution that describe the modulation of TAC absorption and metabolism due to the circadian activity of CYP3A5 and PgP enzymes, despite the fact that the chronopharmacokinetics of TAC affect TAC exposure22 and is one of the causes of IPV7. Various human studies have shown that the maximum blood concentration (Cmax) and AUC (of blood concentration) after the morning dose are higher than the values after the evening dose23, although there are questions regarding the mechanisms involved. Real-time MIPD in TAC prescription provides a customised TAC dosage (output), calculated as a function of previous TAC measurements and other clinical and physiological data (inputs), which presents a real-time optimisation problem that requires a mathematical model able to predict the dynamics of TAC concentration (and other variables) in the desired tissues. The suitability of real-time MIPD is associated with the predictive capability of the mathematical model, which in turn requires the use of adaptive techniques to correct intrapatient variability24. A well-known application of real-time MIPD is glycaemic control through insulin therapy, particularly in intensive care units25. However, the use of real-time MIPD in TAC requires further research. The goal of this study is to develop 2 TAC models based on the PK and PBPK methodologies and to evaluate and compare the models as predictive engines for real-time MIPD. The PK model is based on the model created by Andreu et al.16, whereas the new PBPK model uses the model created by Gerard et al. as a reference21. The specific objectives of this study are (i) to explore the goodness of fit and predictive capability of both models for twice-daily Prograf administration after the inclusion of the chronopharmacokinetics of TAC, testing the TAC circadian mechanistic proposal of Baraldo et al.23 and the improvement in predictive accuracy of PBPK versus PK models18; (ii) to evaluate the capability of several models’ adaptive techniques based on prior parameter distribution to analyse the conversion from Prograf to Advagraf formulation and to calculate intrapatient variability through the parametric mechanisms involved; and (iii) to analyse the change in absolute bioavailability during formulation conversion according to the models’ predictions. Methods The study was performed in 2 phases. During the first phase, we developed, fitted and evaluated population PK and PBPK TAC models using a retrospective stable paediatric population with kidney transplants who were administered Prograf twice daily (Prograf data). During the second phase, we examined the new models’ capability to attend the patients’ progression and evaluated several personalised optimisation strategies to adapt the 3 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ models to the same paediatric patient group, after transitioning from the Prograf to the Advagraf formulation (Advagraf data). We employed adaptive techniques (such as Bayesian estimation)26 based on the models’ parameter distributions obtained during the first phase to adjust the models to the patients’ progression. This strategy explains why our design uses population PK models instead of individually based models, which are limited to more typical adaptive techniques from the adaptive control theory24 applied in pharmacological controls. The population PK and PBPK mathematical models and adaptive techniques are presented in the next 3 subsections, and the computational procedures for the study’s 2 phases are presented in the last subsection (Model Analysis). A glossary with the names of the variables, definitions and units is presented for each mathematical model. Additional details and data regarding the models can be found in theFirst Supplementary Document. The mathematical models were implemented on multilevel object-oriented PhysPK biosimulation software. An evaluation of the accuracy and underlying methodology of PhysPK with respect to NONMEM software was recently published27. The model fitting, covariate analysis, statistics and optimisation procedures described in the next sections are also implemented in PhysPK. We successfully compared the PhysPK population-fitting techniques against NONMEM28. Further details regarding these computational procedures appear in theFirst Supplementary Document. Pharmacokinetic model development. A 2-compartment structural model with a 3-compartment absorption and transit (CAT3) subsystem was the basis of the PK model for TAC distribution. The CAT3 subsystem is based on the model for oral drug absorption defined in the study by Lawrence et al.29, in which the small intestine is addressed as a series of 3 segments, and absorption in the stomach and colon is ignored30. Similar PK models with more compartments were tested without improvement in either TAC prediction or IPV. The model’s mathematical equations are detailed in theFirst Supplementary Document, while Table1 shows the model variables. TAC elimination was assumed to be linear and given by the central clearance Cl . TAC transit in CAT3 is assumed to be linear, with a mean transit time in the stomach, small intestine and colon equal to MTT . The transit rate constant is =kMTT1/ t . The systemic clearance Cl and absorption constant Ka were modulated with a circadian profile according to the suggestions from the Baraldo and Furlanut study23, based on experimental and clinical studies. Chronomodulated parameters were calculated using Eq. (1). =⋅ − =⋅ − Cl tClftt Kt Kftt () () () (), (1) chrwchr achr a w chr in which ft() w is a 24-h periodic waveform function (Fig.1). As shown in Fig.1, the chronomodulated clearance ( Clchr ) is kept equal to Cl during the diurnal phase of the cycle (0 − t1 ) and reduced with the factor ≤r1 chr during the nocturnal phase ( −t t 23 ). Diurnal-nocturnal ( −t t ) 12 and nocturnal-diurnal ( −t t 34 ) phases are completed in 2 h with a constant slope ramp. The value =time0 is defined as the reference time for the first TAC measurement in the clinical study (Methods section). The diurnal phase could therefore start at a different time. This is modelled with the waveform function −ft t() wchr , in which the diurnal phase starts at =t tchr . The reduction factor rchr and the diurnal starting time reference tchr are variable population parameters. The defined circadian profile is the simplest function that supports the circadian behaviour of TAC Cl according to current knowledge. Kachr was defined with the same function fw . The absolute drug bioavailability ( F ) was calculated using Eq. (2): ∫∫ ⋅⋅=⋅Ft Dt tCltCt t() () d()()d (2) t f t chrc 00 in which Df is the drug mass flow rate. The TAC dosage ( D ) was determined with Eq. (3). Variable Units Description Cl L/h Systemic clearance Cld L/h Distribution clearance Vi LVolume of compartment i† Ci ng/ml Concentration at compartment i† Ka ml/h Absorption constant MTT hMean transit time rchr pu Chrono modulation factor tchr hChrono modulation temporal phase D mg Dose kt 1/h Transit rate constant Vcatj, LFluid volume of zone j into CAT‡ kta 1/h Absorption rate constant Table 1. Definition of the Variables employed in the TAC PK model. †i takes values c (central) or p (peripheral). ‡j takes values s (stomach), si (small intestine), and c (colon). 4 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ ∫ =⋅DDtd (3) FDT f 0 Equations (2) and (3) were calculated during the simulation to provide Ft() bioavailability. Physiologically based pharmacokinetic model development. The PBPK TAC model comprises gastric, liver, kidney and primary blood tissues, together with fat tissue and other tissue. The input variables, such as body weight (BW), sex, age, haematocrit and CYP3A4/5 activity, were considered using mechanistic mathematical equations, which are described in this section. The flow diagram is shown in Fig.2. All tissues are nonhomogeneous for TAC, with tissue-plasma partition coefficients associated with the lipophilic character of this drug. The partition coefficients were calculated using the Poulin and Theil method, based on the solubility of drugs in neutral lipids, phospholipids and water31. Physiological lipid and water fractions in tissues were taken from published databases31–33. Figure 1. 24-h periodic waveform ft() w (pu) with ===ttt10,12, 22 123 and =t24 h 4 . Figure 2. Flow diagram for TAC PBPK model with 4 flow-limited tissues (fat, kidneys, liver and others) and 2 membrane-limited tissues (gut and blood). The blood compartment is defined through the red blood cellplasma component. The gastric system is comprised of a gut lumen where the TAC form is liberated following a zero-order kinetic with sink condition, a one-order absorption membrane and gut tissue perfused with blood. 5 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ The model’s mathematical equations follow usual PBPK assumptions. However, they are written in detail in theFirst Supplementary Document. Table2 shows the variables used in the model. TAC plasma binding is defined as the unbound plasma fraction fup , which was 0.01234, despite some variability. This value is commonly used in other modelling studies21,35. An exploratory Monte Carlo simulation was performed to confirm the model’s robustness. Details are presented in theFirst Supplementary Document. TAC demonstrated extensive binding to red blood cells, with a maximum binding capacity Bmax of 418 ± 258 μg/L and dissociation constant ( Kd ) of .±.3847 μg/l for adult liver transplant recipients36. Due to the importance of TAC binding to red blood cells, its variability with different populations (e.g., healthy individuals have approximately double the value of Bmax 36) and the lack of data for paediatric renal transplant recipients, these parameters have been set as population variables to be fitted. We assumed that TAC elimination is due only to liver metabolism and that renal excretion is negligible2,21. We calculated the total intrinsic hepatic clearance as ⋅ClV ih l (ml/h), with Vl the liver volume (ml). Liberation is governed by the zero-order rate liberation constant ( LibT1/ ) in which LibT is the TAC form liberation time (min). Gut absorption is defined by the one-order absorption rate constant Kabs (1/h), with a value of 4.5 1/h36. The one-order elimination rate constant Kelib (1/h) describes the TAC form that cannot be absorbed from the gut lumen, which is related to LibT through Eq. (4), which states that elimination increases with LibT . The relationship is based on the assumption that higher LibT induces the transit of a larger portion of unabsorbed TAC through the intestine. ==⋅KLibTKLibT LibT () (100) 100 (4) elib elib Model parametric fitting in the exploration phase provided better TAC predictions for Kelib defined by Eq. (4) than with the constant Kelib and other polynomial −KLibT elib relationships. We assumed that LibT depends mainly on the release time of the TAC formulation, which is shown in the additional model data section of theFirst Supplementary Document. Intrinsic hepatic clearance ( Clih ) and the gut absorption rate constant ( Kabs ) were modulated with a circadian profile similar to that employed in the PK model. The chronomodulated values Clt() ihchr and Kt() abschr are defined by Eq. (1), replacing Cl and Ka with Clih and Kabs , respectively. The bioavailability F was computed using Eq. (2), replacing ⋅CltC() chrc with ⋅⋅VCltC() lihchr b . Adaptive techniques. We fitted models to our study patients using a maximising population likelihood approach followed by customisation for each patient (Bayesian optimisation26). We then employed adaptive techniques to adjust the models’ parameters to deviations associated with patient temporal intravariations and to modifications in the external context or therapies (such as changes in drug formulations). The adaptive techniques evaluated in this study are as follows: Variable Units Description Vi ml Volume of tissue i† Ci ng/ml i tissue concentration† Cvi ng/ml Venous i tissue concentration† Cut ng/ml Unbound t tissue concentration‡ Cbt ng/ml Bound t tissue concentration‡ Rt pu Tissue (t) - plasma partition ratio‡ Qi ml/min Venous blood perfusion flow rate in i tissue† fup pu Unbound plasma fraction fut pu Unbound tissue fraction‡ Bmax μ g/L Maximum binding capacity in red blood cells Kd μ g/L Dissociation constant in red blood cells Clih 1/h Intrinsic hepatic clearance per liver volume Kabs 1/h Absorption rate constant LibT min TAC liberation time Kelib 1/h Elimination rate constant in liberation region rchr pu Chrono modulation factor tchr hChrono modulation temporal phase D mg Dose Table 2. Variables employed in the TAC PBPK model. †i takes values l (liver), k (kidney), g (gut tissue), f (fat), o (others), b (blood). ‡t takes the same values as i, except b. 6 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ 1. Adaptive Bayesian optimisation. 2. Adaptive weighted least squares (WLS) optimisation. The clinical protocol is described below. Both techniques benefit from the fact that the target patient belongs to a known prior population, which explains why central (mean) parameters remained constant during the customised adaptation phase. The adaptive Bayesian technique follows the same equations of the well-known Bayesian method26, with the exception of boundary conditions. If η is the interpatient variation of the population variable φ (see Model analysis Section) to be computed, this variable is restricted to fulfilling one of the two following boundary conditions (Eq. (5)): ηη ηηη η −⋅ ≤Δ≤⋅ −≤Δ≤ WW BB min( ,) min( ,) , (5) bmbm00 in which η0 is the value of η at the previous model state (time t0 ), Δ ηηη=− 0 is the increment of η in the readjustment and ηm is the minimal η value that defines the weighted boundary. Wb is the boundary Weight and B is the absolute Boundary. Both equations limit the variation of η with respect to the previous state to a weighted value of η0 or to an absolute boundary B. Selection of a B-based or W-based boundary is based on the knowledge concerning the associated population variable φ . The adaptive WSE optimisation minimises the objective function Fwse given by Eq. (6): ∑ =−  FW Ct Ct(ln( ())ln( ())) ,(6) tbk bkwse 2 k in which W is the weight applied to the squared difference between blood model tacrolimus concentration C t() bk and blood measured tacrolimus concentration  Ct() bk , and =…tkn,1,2, k are the sample times. As opposed to the adaptive Bayesian technique, which can only be applied to adjust parameters pertaining to the prior population set, adaptive WLS optimisation can also be used to adjust other models’ variables that were considered constant in the initial design, which helps extend the model’s scope. Adjusted variables and their boundary equations are fully detailed in the Adaptive Modelling Techniques section in theFirst Supplementary Document. Clinical study. The clinical study was designed as a single-centre, open-label TAC treatment conversion from a Prograf to Advagraf formulation (1:1, mg:mg), which included stable paediatric kidney transplant recipients. It was conducted at La Paz University Hospital (referral hospital for paediatric kidney transplantation in Madrid, Spain). The study protocol was approved by the Ethics Committee of La Paz University Hospital and the Spanish Agency for Medicines and Health Products. The study was registered under European Clinical Trial Register EudraCT 2009-017600-89. All methods were performed in accordance with required guidelines and regulations. Written informed consent was obtained prior to the patients’ enrolment from their parents or guardians, who controlled the patient’s drug administration and food intake according to the instructions provided. A full detailed description of the clinical study design, including the results of the evaluation of the relative bioavailability of Prograf and Advagraf formulations, has been published37. Data collected during the pharmacokinetic phase of the study (2 weeks), which is described in the following paragraphs, support the development of the present research with the aim of reaching the 3 specific objectives stated in the Introduction Section. Patients with a stable renal transplant and younger than 18 years of age were recruited according to the inclusion criteria. The study size was defined as 21 patients, although only 18 patients were needed to achieve statistically significant steady-state AUC24 (AUC of blood TAC concentration during the last 24 h) target differences of 5% with a power of 80%37. The inclusion criteria included a variation in TAC dosage and concentrations lower than 20% during the last 30 days prior to the start of the study, without administering drugs that might affect the TAC PK in the last 15 days. None of the patients dropped out during the pharmacokinetic study period, and none of the patients experienced adverse events during the 2 weeks of the PK study. None of the concomitant treatments were modified during the study. The mean time from transplant was .±.539 325 years, and nearly all patients were white (17). The remaining patients were Hispanic (1), Asian (2) and Arabian (1). The sex distribution was 57% male/43% female. The other anthropometric data (expressed as mean ± SD, range) were age (12.9 ± 4.17 years, 4–17), body weight (42.85 ± 15.42 kg, 15.1–63.8) and height (143.4 ± 18.16 cm, 105–168). The patients were administered a controlled Prograf regimen the week prior to their first admission to La Paz University Hospital (day 1). The patients remained on Prograf during day 1, which was first administered between 8:00 AM and 9:00 AM and then re-administered 12 h later. The drug was administered with 200 ml of water, and the patients ate breakfast, lunch and dinner 11:00 AM, 1:00 PM and 7:00 PM. During the hospital stay, blood TAC was measured just before the morning TAC administration ( =t0 ) and consecutively 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 12.5, 13, 14, 15, 16, 18, 20 and 24 h later (a total of 18 measurements). Blood samples (2 ml) were collected in aliquot duplicate Ethylene Diamine Tetraacetic Acid (EDTA) tubes and analysed with an enzyme immunoassay (Dimension, Siemens Healthcare Diagnostics, Frimley, Camberley, UK), with 2 ng/ml as the lower quantification limit and 30 ng/ml as the upper limit. Measurements above the upper limit were diluted to obtain concentrations in the valid range. Prograf was replaced with Advagraf on the following morning (day 2) with a new controlled protocol. Once-daily Advagraf doses were based on 1:1 total dose equivalence. Patient data acquired during the Prograf 7 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ phase of the study (day 1) were used to fit and evaluate the new population PK and PBPK TAC models. The patients were re-admitted to La Paz University Hospital on day 8, repeating the food intake and blood sampling procedures of the first monitoring day, with the following sampling times: 0, 0.5, 1, 1.5, 2, 2.5, 3, 4, 6, 8, 12, 15 and 24 h (13 measurements). The blood TAC measurement ( =t0 ) of this second series was performed just before administering Advagraf. The sample was the trough blood TAC concentration (C0) at the beginning of the hospital re-admission day. The model’s parameters were readjusted to reduce blood TAC concentration prediction errors during day 8. The adaptive optimisation techniques for the model readjustment are presented in the Adaptive techniques Section. The TAC measurements and data used for the analysis are included in theSecond Supplementary File. Models analysis. Population model structure. 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.16, 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. (7) for the Cl parameter. θ=⋅ η Cl e(7) Cl Cl 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 Cc in PK and ln Cb 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 considered38. 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. Covariates analysis. 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 approach39. 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: θθ θ=⋅ +⋅− ˆ xx(1 ()) (8) ii cr whereas the exponential model is as follows: θθ=⋅             θ ˆx x, (9) ii r c 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. Population model fitting techniques and evaluation. 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, η -shrinkage42 and relative standard errors (RSEs) of estimated parameters. 8 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ 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)43 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)44 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 Keizer45, 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 misspecification46. 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 concentration47 were also added to the calculated exposures. Adaptive procedures evaluation. 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): (i) 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. (ii) 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 RMSEc43. We compared the results of these 2 procedures, along with the procedures’ effect on the models’ suitability as predictive engines for dosage recommendations. Results The following two subsections present the results for the first phase of the study (Prograf data). The last subsection (Adaptive Models) shows the results for the second phase (Advagraf data). Fitting and evaluation of PK model (Prograf data). Our initial examination of the model confirmed its structural stability, robustness and ability to provide initial estimates of the parameters to be fitted. More details of the analysis appear in the Results section of theFirst Supplementary Document. We applied the FOCE technique without covariates to the additive and proportional residual error PK models. The additive residual model obtained fewer RSEs for all mean parameters ( θi ). The two models had similar Ω variance estimates, except for ωKa 2 , whose additive value was one-third of the proportional value. The two models also had a similar accuracy for the TAC estimations. We therefore selected the additive residual error PK model. The MVOF obtained with the PK reduced model (without covariates) was −885.13. We applied the covariate forward stepwise approach. The inclusion of the BSA − Cld covariate regression exponential model reduced the MVOF to −888.08, whereas the MVOF was virtually unchanged when the Chb − Vp covariate regression model was added. In addition, the BSA covariate reduced interindividual variance (IIV) by an average of 70% (CV). Given that the reduction in MVOF (2.95) was near 3.84 (statistic significance level of 5%), we decided to accept BSA as one of the model’s covariate. The value of the estimated exponent θc for exponential correlation (Eq. (9)) was 1.17. 9 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ The mean estimated parameter values θi and associated interindividual variances, together with their RSEs, for the PK TAC model are presented in Table3. CV was approximated as ω i 2 for small η variance ( ωi2 ), resulting in CV values of 31% for Cl , 57% for Ka and 22% for MTT . A first fitting diagnostic was based on shrinkage metric indices42. The shrinkage value related to residual error was ε -shrinkage = 0.014 ± 0.349, and the η -shrinkage parameters were 0.18, −0.08 and 0.17 for Cl , Ka , MTT , respectively, and 0.66 for rchr . Values near zero confirm the goodness of fit. Other relevant values obtained during the Bayesian estimation were =.±.F18 897% , =.±.AUC 200 3407 24 ng/ml · h (model-based), =.±.AUC 201 48 39 27 24 ng/ml · h (NCA log-trapezoidal), and Δ =− .±.AUC 122467 ng/ml · h. The accuracy of TAC concentration (log-transformed) referred to the whole set of samples was =. ±.RMSE 0074 0032 c . These metric values support the goodness of fit of the PK TAC model. Plots of model TAC concentration against TAC measurements for several randomly selected patients and a table that shows a comparison against parameters obtained by the model created by Andreu et al.16 are presented in theFirst Supplementary Document. Parameters have similar and reasonable values, with the exception of the peripheral volume Vp , which is 39.1 l in our model against 526.03 l in the Andreu et al. model. This issue is addressed in the Discussion section. The two models had similar interindividual dispersion, despite the fact that our PK model was fitted for paediatric patients and includes diurnal and nocturnal TAC distribution. Therefore, our PK model demonstrates greater predictive capability. Figure3 shows the pcVPC based on 500 replicas according to the obtained parametric population. TAC concentrations have been normalised with D/BSA mg/m2, according to the wide range of TAC dosages for this paediatric population. As shown, the observation percentiles agree with the predicted CIs obtained from PK model simulations. The results for the IIVs agree with those of other studies16, despite the studies not requiring TAC predictions during 2 sequential 12-h Prograf administrations. Additional studies, showing the fulfilment of a normal distribution assumption and the results of a Monte Carlo simulation with 600 runs for the PK model with 3 TAC normalised doses can be found in theFirst Supplementary Document (PK results). According to the Monte Carlo study, a dose of 0.042 mg/kg achieves an AUC24 exposure of 271.4 ± 196.7 ng/ml · h. The high dispersion of AUC24 emphasises the need to use MIPD to customise the TAC dose. Fitting and evaluation of PBPK model (Prograf data). Our initial exploration of the PBPK model confirmed its structural stability and robustness, which also provided the initial estimates for the parameters to be fitted. More details of the analysis can be found in the Results section of theFirst Supplementary Document. The accuracy of the TAC concentrations was similar with the residual and proportional error, as well as IIV estimates; we therefore selected the additive residual error. Table4 presents the mean estimated parameter values θi and their interindividual variances (along with their RSEs) calculated for the PBPK model using the FOCE technique. The CVs were calculated in the same manner as in PK model, with values 33% for Clih , 6.3% for Bmax and 28.8% for Kd . Shrinkage index for residual error was ε -shrinkage = 0.039 ± 0.30, and η -shrinkages were −0.32, 0.80 and 0.72 for ClB K,, ih maxd , respectively, and 0.12 for rchr . Values near zero confirm the goodness of fit42. Parameter Units Value (RSE %)† Disposition parameters θCl (L/h) 5.41 (5.7) θCLd (L/h) 46.5 (7.5) θVc (L) 23.0 (7.9) θVp (L) 39.1 (7.2) Circadian parameters θrchr (pu) 0.29 (6.4) θtchr (h) 0.7 (−) Absorption parameters θKa (ml/h) 13.8 (7.3) θMTT (h) 0.72 (5.5) Interindividual variabilities IIVCl —0.10 (20) IIVKa —0.33 (18) IIVMT —0.05 (40) IIVchrred —0.003 (30) IIV chrp hi —0.01 (−) Residual error variance ωr2 —0.027 (3.7) Table 3. Population PK model parameters’ mean and IIV estimates: value (RSE), for Prograf fitting. †RSE is calculated through the standard error SE as SE(θ)/θ · 100 for θ and SE ωω⋅ θθ ()/(2)100 22 for variances. 16 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ The conclusions of this study are limited to the controlled and stable renal transplant paediatric population. Mathematical models are limited to the functional scope defined in the Model Development section. The interoccasion variability was not implemented in our models as justified in the Methods Section, despite this term is commonly included in population PK models and it can affect the TDM accuracy. However, our outcomes suggest that circadian rhythms explain a significant portion of IPV associated with IOV. In addition, the exclusion of IOV term may be a preferred approach to maximize precision in the calculation of individual doses using model-based Bayesian forecasting, according to the recent study by Abrantes et al.55. Moreover, adaptive mathematical models have demonstrated their ability to predict customized biochemical levels, such as glycemia, without the inclusion of IOV25. We did not consider DDIs at the present time. Further research is needed to add unstable patients and extend models to consider DDIs. Current knowledge on physiological and pharmacokinetic modelling suggests that these extensions are feasible. Data availability All data generated or analysed during this study are included in this published article (and its Supplementary Information files). Received: 11 November 2019; Accepted: 10 April 2020; Published: xx xx xxxx References 1. Andrews, L. M. et al. Pharmacokinetic considerations related to therapeutic drug monitoring of tacrolimus in kidney transplant patients. Expert Opin Drug Metab Toxicol 13, 1225–1236 (2017). 2. Rahman, Z., Zidan, A. & Khan, M. A. Tacrolimus properties and formulations: potential impact of product quality on safety and efficacy, book section Chapter 1, 1–39 (Nova Science Publishers, Inc., 2013). 3. Chen, H. & Qian, S. Current Immunosuppressive Therapy in Organ Transplantation. (Nova Science Publishers Inc., New York, 2015). 4. Iwasaki, K. Metabolism of tacrolimus (fk506) and recent topics in clinical pharmacokinetics. Drug metabolism and pharmacokinetics 22, 328–335 (2007). 5. DrugBank. Drugbank - tacrolimus (https://www.drugbank.ca/drugs/db00864) (2018). 6. Brunet, M. et al. Therapeutic drug monitoring of tacrolimus-personalized therapy: Second consensus report. Therapeutic Drug Monitoring 41, 261–307 (2019). 7. Shuker, N., van Gelder, T. & Hesselink, D. A. Intra-patient variability in tacrolimus exposure: causes, consequences for clinical management. Transplant Rev. (Orlando) 29, 78–84 (2015). 8. Borra, L. C. P. et al. High within-patient variability in the clearance of tacrolimus is a risk factor for poor long-term outcome after kidney transplantation. Nephrology Dialysis Transplantation 25, 2757–2763 (2010). 9. Süsal, C. & Döhler, B. Late intra-patient tacrolimus trough level variability as a major problem in kidney transplantation: A collaborative transplant study report. American Journal of Transplantation 19, 2805–2813 (2019). 10. Shah, P. B., Ennis, J. L., Cunningham, P. N., Josephson, M. A. & McGill, R. L. The epidemiologic burden of tacrolimus variability among kidney transplant recipients in the united states. American Journal of Nephrology 50, 370–374 (2019). 11. Darwich, A. S. et al. Why has model-informed precision dosing not yet become common clinical reality? lessons from the past and a roadmap for the future. Clin Pharmacol Ther. 101, 646–656 (2017). 12. Wallin, J. E., Friberg, L. E., Fasth, A. & Staatz, C. E. Population pharmacokinetics of tacrolimus in pediatric hematopoietic stem cell transplant recipients: new initial dosage suggestions and a model-based dosage adjustment tool. Therapeutic drug monitoring 31, 457–466 (2009). 13. Størset, E. et al. Improved prediction of tacrolimus concentrations early after kidney transplantation using theory-based pharmacokinetic modelling. British Journal of Clinical Pharmacology 78, 509–523 (2014). 14. Ot e o, I. et al. Tacrolimus pharmacokinetics in the early post-liver transplantation period and clinical applicability via bayesian prediction. European Journal of Clinical Pharmacology 69, 65–74 (2013). 15. Andrews, L. M. et al. A population pharmacokinetic model to predict the individual starting dose of tacrolimus following pediatric renal transplantation. Clin Pharmacokinet 57, 475–489 (2017). 16. Andreu, F. et al. Development of a population pk model of tacrolimus for adaptive dosage control in stable kidney transplant patients. Therapeutic Drug Monitoring 37, 246–255 (2015). 17. Benka li, K. et al. Population pharmacokinetics and bayesian estimation of tacrolimus exposure in renal transplant recipients on a new once-daily formulation. Clinical Pharmacokinetics 49, 683–692 (2010). 18. Bassingthwaighte, J. B., Butterworth, E., Jardine, B. & Raymond, G. M. Compartmental modeling in the analysis of biological systems. Methods Mol. Biol. 929, 391–438 (2012). 19. Barrett, J. S., Della Casa Alberighi, O., Laer, S. & Meibohm, B. Physiologically based pharmacokinetic (pbpk) modeling in children. Clin Pharmacol Ther 92, 40–9 (2012). 20. Zhang, H. et al. Prediction of drug–drug interaction between tacrolimus and principal ingredients of wuzhi capsule in chinese healthy volunteers using physiologically-based pharmacokinetic modelling. Basic & Clinical Pharmacology & Toxicology 122, 331–340 (2018). 21. Gérard, C. et al. Determination of the most influential sources of variability in tacrolimus trough blood concentrations in adult liver transplant recipients: a bottom-up approach. The AAPS journal 16, 379–91 (2014). 22. Barbarino, J. M., Staatz, C. E., Venkataramanan, R., Klein, T. E. & Altman, R. B. Pharmgkb summary: cyclosporine and tacrolimus pathways. Pharmacogenet Genomics 23, 563–85 (2013). 23. Baraldo, M. & Furlanut, M. Chronopharmacokinetics of ciclosporin and tacrolimus. Clin. Pharmacokinet 45, 775–88 (2006). 24. Zhang, P. Chapter 2 - Industrial control engineering (pp. 41–70. William Andrew Publishing, Oxford, 2010). 25. Herpe, T. V. et al. An adaptive input-output modeling approach for predicting the glycemia of critically ill patients. Physiol. Meas. 27, 1057 (2006). 26. Kelman, A. W., Whiting, B. & Bryson, S. M. Parameter optimisation in clinical pharmacokinetics. Computer Programs in Biomedicine 14, 239–248 (1982). 27. Reig-Lopez, J., Merino-Sanjuan, M., Mangas-Sanjuan, V. & Prado-Velasco, M. A multilevel object-oriented modelling methodology for physiologically-based pharmacokinetics (pbpk): Evaluation with a semi-mechanistic pharmacokinetic model. Computer Methods and Programs in Biomedicine 189, 1–11 (2020). 28. Gonzalez-Garcia, I., Prado-Velasco, M., Fernández-Teruel, C. & Fudio, S. Comparison of fo – foce population parameter estimation methods in physpk 2.0 against nonmem 7.3. In PAGE 2017. Abstracts of the Annual Meeting of the Population Approach Group in Europe, 1–2 (2017). 29. Yu, L. X. & Amidon, G. L. A compartmental absorption and transit model for estimating oral drug absorption. International Journal of Pharmaceutics 186, 119–125 (1999). 17 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ 30. Tsunashima, D. et al. Assessment of tacrolimus absorption from the human intestinal tract: Open-label, randomized, 4-way crossover study. Clinical Therapeutics 36, 748–759 (2014). 31. Poulin, P. & Theil, F. P. A priori prediction of tissue: Plasma partition coefficients of drugs to facilitate the use of physiologicallybased pharmacokinetic models in drug discovery. Journal of Pharmaceutical Sciences 89, 16–35 (2000). 32. ILSI. Physiological parameter values for pbpk models. Report, International Life Sciences Institute (1994). 33. Pelekis, M., Poulin, P. & Krishnan, K. An approach for incorporating tissue composition data into physiologically based pharmacokinetic models. Toxicology and Industrial Health 11, 511–522 (1995). 34. Zahir, H., Nand, R. A., Brown, K. F., Tattam, B. N. & McLachlan, A. J. Validation of methods to study the distribution and protein binding of tacrolimus in human blood. Journal of Pharmacological and Toxicological Methods 46, 27–35 (2001). 35. Emoto, C. et al. A theoretical physiologically-based pharmacokinetic approach to ascertain covariates explaining the large interpatient variability in tacrolimus disposition. CPT Pharmacometrics Syst Pharmacol 8, 273–284 (2019). 36. Jusko, W. J. et al. Pharmacokinetics of tacrolimus in liver transplant patients. Clinical Pharmacology & Therapeutics 57, 281–290 (1995). 37. Carcas-Sansuán, A. J. et al. Conversion from prograf to advagraf in stable paediatric renal transplant patients and 1-year follow-up. Pediatric Nephrology 29, 117–123 (2014). 38. Bonate, P. L. Pharmacokinetic-Pharmacodynamic Modeling and Simulation, second edi edn. (Springer, 2011). 39. Jonsson, E. N. & Karlsson, M. O. Automated covariate model building within nonmem. Pharmaceutical Research 15, 1463–1468 (1998). 40. Davidian, M. Introduction to statistical population modeling and analysis for pharmacokinetic data. In Invited white paper for the International Workshop on Uncertainty andVariability in Physiologically Based Pharmacokinetic (PBPK) Models 89. 41. Wang, Y. Derivation of various nonmem estimation methods. Journal of Pharmacokinetics and Pharmacodynamics 34, 575–593 (2007). 42. Savic, R. M. & Karlsson, M. O. Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions. The AAPS journal 11, 558–69 (2009). 43. Barnston, A. G. Correspondence among the correlation, rmse, and heidke forecast verification measures; refinement of the heidke score. Weather and Forecasting 7, 699–709 (1992). 44. Bergstrand, M., Hooker, A. C., Wallin, J. E. & Karlsson, M. O. Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models. The AAPS journal 13, 143–151 (2011). 45. Keizer, R. Create visual predictive checks (https://www.rdocumentation.org/packages/vpc/versions/1.1.0) (2018). 46. Hooker, A. C., Staatz, C. E. & Karlsson, M. O. Conditional weighted residuals (cwres): A model diagnostic for the foce method. Pharmaceutical Research 24, 2187–2197 (2007). 47. Sikma, M. A. et al. Pharmacokinetics and toxicity of tacrolimus early after heart and lung transplantation. Am J. Transplant 15, 2301–13 (2015). 48. CHMP. Advagraf, inn-tacrolimus (https://www.ema.europa.eu/medicines/human/epar/advagraf’) Report, EMA (2018). 49. Cho, E. J. et al. Performance of the dimension tac assay and comparison of multiple platforms for the measurement of tacrolimus. J. Clin. Lab. Anal. 32, e22357 (2018). 50. Shin, T. H. et al. Formulation and in vivo pharmacokinetic evaluation of ethyl cellulose-coated sustained release multiple-unit system of tacrolimus. Int. J. Biol. Macromol. 109, 544–550 (2018). 51. Leino, A. D. et al. Assessment of tacrolimus intrapatient variability in stable adherent transplant recipients: Establishing baseline values. American Journal of Transplantation 19, 1410–1420 (2019). 52. Toutain, P. L. & Bousquet-Mélou, a. Bioavailability and its assessment. Journal of veterinary pharmacology and therapeutics 27, 455–66 (2004). 53. Staatz, C. E. & Tett, S. E. Clinical pharmacokinetics of once-daily tacrolimus in solid-organ transplant patients. Clinical Pharmacokinetics 54, 993–1025 (2015). 54. Prado-Velasco, M. et al. Pbpk versus pk modeling of tacrolimus drug in patients with renal transplant as knowledge engines for personalized posology software: Physpk® development and preliminary results. In PAGE 2017. Abstracts of the Annual Meeting of the Population Approach Group in Europe, 1–2 (2017). 55. Abrantes, J. A., Jönsson, S., Karlsson, M. O. & Nielsen, E. I. Handling interoccasion variability in model-based dose individualization using therapeutic drug monitoring data. British Journal of Clinical Pharmacology 85, 1326–1336 (2019). Acknowledgements This research study was partly supported by a collaboration agreement (2018) with Empresarios Agrupados Internacional (EAI). Author contributions A.B. and A.C.S. designed the clinical study; M.P.V., A.B. and A.C.S. contributed to the concept and design; M.P.V. designed and executed the mathematical models; and M.P.V., A.B. and A.C.S. contributed to the analysis of results. competing interests Dr. Prado-Velasco received compensation as team leader, designer and main developer of PhysPK biosimulation software in EAI (Oct 2015–March 2018). All authors participated in the collaboration agreement (2018) with EAI. Additional information Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-64189-9. Correspondence and requests for materials should be addressed to M.P.-V. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 18 Scientific RepoRtS | (2020) 10:7542 | https://doi.org/10.1038/s41598-020-64189-9 www.nature.com/scientificreports www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2020