Pharmacokinetic Model Based on Stochastic Simulation and Estimation for Therapeutic Drug Monitoring of Tacrolimus in Korean Adult Transplant Recipients : Therapeutic Drug Monitoring

Secondary Logo

Journal Logo

Original Article

Pharmacokinetic Model Based on Stochastic Simulation and Estimation for Therapeutic Drug Monitoring of Tacrolimus in Korean Adult Transplant Recipients

Choi, Suein MD*,†; Hong, Yunjeong MS*,†; Jung, Sook-Hyun MS; Kang, Gaeun MD, PhD§; Ghim, Jong-Ryul MD, PhD; Han, Seunghoon MD, PhD*,†

Author Information
doi: 10.1097/FTD.0000000000001006
  • Open



Tacrolimus is a calcineurin inhibitor widely used as an immunosuppressant; it is also used as the first line of treatment after kidney and liver transplantation in Korea. Tacrolimus is primarily metabolized in the liver before elimination and is rarely excreted through the urine. Although its half-life is reportedly approximately 12 hours, it is highly variable from 3.5 to 40.5 hours depending on patient demographics and conditions.1 The absorption rate is also variable, reaching peak blood concentrations in 0.5–6 hours with approximately 25% bioavailability.1 Tacrolimus has a narrow therapeutic index,1 and the correlation between blood concentrations and doses is poor, indicating highly variable inter- and intraindividual pharmacokinetics (PK).2–4 Various factors, including ethnicity, sex, age, transplant organ, disease, steroid usage, and time after treatment initiation, contribute to the variability in tacrolimus PK.5 Therefore, therapeutic drug monitoring (TDM) is a general clinical practice in patients with a transplant to maintain effective therapeutic levels and to avoid toxicity.3

In standard TDM procedures, the selection of tacrolimus dosage based on “a priori model-based individual dosing” and optimizing individual PK parameters using the Bayesian estimation method are critical for improving targeted exposure in patients. The population PK model plays the most important role.6 This technique focuses on estimating individual PK parameters that maximize the likelihood of an observed concentration considering both fixed and random effects. Thus, for desirable TDM results, the PK model should include clinically relevant covariates and the proper magnitude of random effect parameters for the population subjected to TDM. However, for Korean adult transplant recipients, there are only a few full-PK studies that can be used to build a population PK model.7,8 These were based on the data from a single hospital, which could not be considered a representative population. Moreover, the models require covariates that are not routinely captured in Korean clinical practice (eg, genotype and hematocrit in some centers during later periods after transplant). These factors limited the utilization of TDM in many transplantation centers.

Therefore, in this study, we aimed to establish a population PK model with enhanced clinical utility regarding covariate application and individual PK parameter distributions. A tacrolimus population PK model was built from a combined multicenter data set containing simulated data from the existing literature and data collected from several institutions. Because most of the actual data in the TDM procedure were trough concentrations and did not provide sufficient information to build a structural model, the simulated data sets were used as the backbone for model development. Furthermore, instead of using a single simulated data set, the stochastic simulation and estimation (SSE) method was used to obtain a more generalized outcome with the precision of the estimated PK parameter values.


Strategy for Model Development

The overall strategy of this analysis is shown in Figure 1. This study comprised 2 steps. In the first step, we integrated an observed data set with a simulated data set from previously published models.7,8 We developed the final model from an integrated data set using a previously built structural model7 and performed covariate analysis. In the second step, the identical integration and parameter estimation procedure was repeated 1000 times using the final model built in sep 1. The final parameter was acquired by summarizing 1000 sets of estimated parameters. Finally, the tacrolimus PK model was developed using the final structural model built in step 1 and the summarized final parameters estimated in step 2. The ratio of individual numbers in the observed and simulated data sets was 9:1.

Overall strategy of the analysis.

Observed Data set

Tacrolimus concentration data of Korean patients older than 18 years, who underwent kidney or liver transplantation at the Inje University Busan Paik Hospital (IUBPH) and Chonnam National University Hospital (CNUH) between 2005 and 2019, were collected. The data acquisition was approved by the institutional review board of each hospital (approval number: 12-194 for IUBPH and CNUH-2018-099 for CNUH). Blood samples were collected before the administration of the morning dose. The dosage was changed, if needed, based on the tacrolimus trough concentration and the patient's target therapeutic range. Generally, the target therapeutic range was 8–12 ng/mL throughout the first postoperative month, 8–10 ng/mL for up to 3 months, 6–8 ng/mL until 6 months, and 4–6 ng/mL thereafter. However, the target therapeutic range could differ depending on the type of transplant organ, patient status, institution, combination drug regimen, and clinician's decision. Sample analysis was performed using affinity column-mediated immunoassay on the Dimension EXL 200 system (Siemens, Berlin, Germany) and chemiluminescence microparticle immunoassay using the ARCHITECT i2000 system (Abbott Laboratories, Abbott Park, IL) at IUBPH and CNUH, respectively. The range of quantification was 2–30 ng/mL. Because the concentration results from these methods are generally accepted and used for clinical decision making, the data from the 2 hospitals were considered compatible despite possible quality differences for accuracy and precision. Age, sex, weight, height, drug name, transplant organ, and the number of postoperative days (PODs) were considered potential covariates, and missing values were imputed using the missForest package (ver1.4)9 in R (ver4.0.3).

Simulated Data set

Two published PK models were used in the simulation.7,8 These studies were performed using data from patients in Korea who had undergone kidney transplantation and received an immunosuppressive regimen, including tacrolimus (Prograf, Astellas Pharma Korea, Seoul, Korea). One model was based on prospectively collected trough concentrations from 122 patients during the first month after transplantation and densely sampled concentrations from 55 patients (predose, 0.5, 1, 2, 4, 6, 8, and 12 hours after administration) 10–15 days after transplantation7 (PK model for an early posttransplant period). The second model was based on retrospectively collected trough concentrations from 80 patients during the follow-up period of 400 days after transplantation (PK model for the late posttransplant period).8 Simulations were performed to reproduce the data by reflecting their sampling points (9 points for the early period and 5 points for the late period in each virtual individual). Concentrations outside the quantification limit for the observed concentrations were omitted. Covariate values included in the models were randomly sampled based on the distribution of corresponding variables in the previous reports. In each simulated data set, 30 individuals were included for both early and late posttransplant periods with 420 simulated concentrations (270 from the early period and 150 for the late period). This number of simulated subjects was considered the minimum number to maintain the general trend of the observed data and to secure the condition for estimating the magnitude of between-subject variability (ω2), which is assumed to be normally distributed (according to the central limit theorem).

Base Model Development

The base model was developed using the combined data of an observed and a simulated data set using population PK analysis. Population PK analysis was performed using the first-order conditional estimation method with an interaction (FOCE + I) in nonlinear mixed-effect modeling (NONMEM, version 7.5, ICON Development Solution, Ellicott City, MD). R, Rstudio (ver1.4.1), and Xpose4 package (ver4.7.1)10 were used for graphical analyses and model diagnostics.11

Based on the population PK models from original articles, which were used as the backbone in this study, one-compartment PK models with first-order and lag time were selected as the structural model.7,8 The interindividual variability of each PK parameter was assumed to follow the log-normal distribution and described using an exponential model as follows:(1)Pi=PTV×exp(ηi),where Pi denotes the individual parameter (eg, CL/F or V/F) for the i-th individual, PTV is the typical value of the model parameter for the population, and ηi is the interindividual random effect following a normal distribution with a mean of zero and variance of ω2, accounting for the deviation of the i-th individual from the typical value of PTV. The OMEGA BLOCK option was used to evaluate the covariance between the random effects. Residual variability for each institution was evaluated using proportional, additive, and combined models. The models were assessed based on objective function value and goodness-of-fit plots.

Covariate Analysis

Covariate analysis was performed on variables (sex, age, body weight, transplant organ type, POD, and drug name) that could be acquired in the clinical TDM practice and expected to explain the interindividual and interoccasion variabilities. The effect of the data source (literature, IUBPH, or CNUH) was also considered. A generalized additive model (GAM) and graphical analysis were used for covariate screening, and a forward selection-backward elimination method with the NONMEM minimization process was used for covariate selection (forward selection P value < 0.05, backward elimination P value < 0.01).

Evaluation of the PK Model Performance for SSE

The final base model was assessed by bootstrapping and visual predictive check validation. The medians and fifth–95th percentiles of the results from 1000 random samplings were compared with the original parameters in bootstrapping. The observed concentrations and 90% prediction intervals of the simulated concentrations versus time were plotted, and any bias or trend was evaluated. Bootstrapping and visual predictive check (VPC) validation processes were performed using R, NONMEM, and Wings for NONMEM (WFN, ver750,

Stochastic Simulation and Estimation

As shown in step 2 (Fig. 1), 1000 data sets consisting of 30 subjects for each model were simulated using published PK models with sampled covariates from their raw data. Subsequently, the simulated data were combined with the observed data set. Finally, the population PK parameters were estimated based on the final PK model for each of the 1000 combined data sets using the FOCE-I method in NONMEM. Thus, 1000 sets of population PK parameters were estimated, and the median of the successfully estimated values of each parameter was selected as the final parameter.

All statistical analyses were performed using the R software (ver4.0.3). The results of the SSE analysis are presented as mean ± 90% confidence interval and relative standard error (RSE) to assess the robustness of the estimated parameters.

Final Model and Parameter Evaluation

Visual predictive check validation (VPC) was performed to evaluate the accuracy and robustness of the final model and estimated parameters from the SSE analysis. Because each patient included in the data set had a different dosage regimen, sampling time, and covariate values, a prediction-corrected VPC (pcVPC) suggested by Bergstrand et al12 was conducted instead of a traditional VPC. Traditional VPC is known for deceiving the analysis of TDM data because its dose adjustments and the observed dependent variable are inherently correlated.12 The pcVPC method is useful, as it helps avoid misleading diagnoses because the variability caused by independent variables (time, dosage, and other covariate values) within a bin can be corrected by normalizing the dependent variable based on the population prediction.12 Population prediction correction was performed on log-transformed dependent variables, as shown in Equation (2):12(2)ln(pcYij)=ln(Yij)+(ln(PREDbin)ln(PREDij)),where pvYij is the prediction-corrected observation of the prediction for the ith individual at the jth time point, Yij is the observation or prediction for the ith individual and jth time point, PREDbin is the median of typical population predictions for the specific bin of independent variables, and PREDij represents the typical population predictions for the ith individual at the jth time point. In total, 1000 simulated replicates of the original design were used to determine a 90% prediction interval. Perl-speaks-NONMEM (PsN, ver4.8.1, and Rstudio with R were used for the simulation and pcVPC.


Population PK Modeling for the Final Structural Model

For the analysis, 1829 tacrolimus observations (1535 from IUBPH and 294 from CNUH) from 269 patients (122 from IUBPH and 147 from CNUH) who underwent tacrolimus TDM were included as an observed data set and were combined with the simulated data set. In this process, observations that were not within the quantification limit were excluded (4 observations). Patient demographic information from the institutions is presented in Table 1. The characteristics of the patients at each institution were comparable. As the original articles suggested, tacrolimus PK data were well-described by a one-compartment, first-order elimination model with lag time. The final parameter values from the structural model were within a reasonable range; POD, age, and organ transplantation were selected as significant covariates. The goodness-of-fit plots and VPC of the PK model with all covariates suggest that the final model adequately describes the observed data. The final estimated PK parameters are presented in Table 2.

TABLE 1. - Patient Demographics
Characteristics Number of Patients
 No. of patients 122
 Male/female 64/58
 Age (yr) 46.6 ± 12.3
 Body weight (kg) 63.2 ± 17.0
 Height (cm) 165 ± 8.90
Type of drug
 Tacrobell 56 (48.7%)
 Prograf 59 (51.3%)
Type of transplantation
 Kidney 114 (93.4%)
 Liver 4 (3.28%)
 Number of patients 147
 Male/female 89/58
 Age (yr) 49.2 ± 11.2
 Body weight (kg) 67.1 ± 14.1
 Height (cm) 164 ± 8.23
Type of drug
 Tacrobell 0 (0.00%)
 Prograf 147 (100%)
Type of transplantation
 Kidney 125 (85.0%)
 Liver 22 (15.0%)

TABLE 2. - Parameter Estimates of Final PK Model
Parameter Description Estimate %RSE Bootstrap Median (95% CI)
Structural model
   Oral clearance (CL/F) before surgery (L/h) in kidney transplant patients 21.9 12.2 22.3 (12.5,28.8)
   Baseline CL/F after surgery in kidney transplant patients (L/h) 20.5 5.9 20.3 (18.1,22.9)
   Baseline CL/F after surgery in liver transplant patients (L/h) 12.9 17.9 12.7 (7.74,20.2)
   Age effect −0.397 27.7 −0.443 (−0.767, −0.191)
   POD effect in early period in kidney transplant patient 0.112 13.3 0.107 (0.0748, 0.138)
   POD effect in liver transplant patient 0.0557 63.2 0.0589 (0.00161, 0.195)
   POD effect in late period in kidney transplant patient −0.023 13.3 −0.0237 (−0.0506, −0.00148)
   Apparent volume (L) 499 12.2 482 (355, 647)
   Body weight effect 1.30 34.6 1.26 (0.0756, 2.29)
   Absorption rate constant 2.59 14.3 2.60 (1.96, 3.47)
   Lag time of first-order absorption 0.25 (fix)
Interindividual variability (CV%)
 ωCL/F Interindividual variability of CL/F 43.4 6.4 50.8 (39.7, 67.5)
 ωV/F Interindividual variability of V/F 84.2 11.9 89.6 (30.0, 118)
Correlation coefficients
 ρCL/F-V/F Correlation coefficients between CL/F and V/F 0.221 40.5 0.198 (0.0186, 0.281)
Residual error
 σadd, 1 * Additive error 0.0001 (fix)
 σprop, 1 * Proportional error 0.244 3.6 0.244 (0.227, 0.262)
 σadd, 2 * Additive error −2.77 13.7 −2.91 (−3.5, −2.1)
 σadd, 3 * Additive error 2.78 10.3 2.71 (2.34, 3.03)
*Number refers to the hospital where data were collected (1: Seoul National University Hospital, 2: Inje University Busan Paik Hospital, 3: Chonnam National University Hospital).
PK, pharmacokinetics; POD, postoperative days; CI, confidential interval; CL, clearance; F, bioavailability; V, volume of distribution; WT, body weight.

POD, age, and transplantation organ (liver and kidney) were selected as statistically significant covariates of apparent clearance (CL/F), and weight was included as a significant covariate of the apparent volume of distribution (V/F). POD had a significant effect on CL/F in both kidney and liver transplant patients, but the effect differed with the transplant organs. In kidney transplant patients, after an immediate decrease post transplantation, CL/F increased as POD increased by 40% within 4 weeks after transplantation. However, the CL/F started to decrease gradually after 1 month until it reached a plateau. In liver transplant patients, baseline CL was lower than that of kidney transplant patients and recovered to a similar value in the late period. The changes in CL/F by POD in both kidney and liver transplant patients are shown in Figure 2. The age and body weight of the patients affected the CL/F and V/F, respectively. As age increased, CL/F decreased, and as body weight increased, V/F increased in both kidney and liver transplant patients.

Clearance (CL/F) profile of 1000 simulated patients for the 1000-day posttransplantation by transplant organ type. Red line represents the mean CL/F in liver transplant patients. Blue line represents the mean CL/F in kidney transplant patients Semitransplant red field represents prediction interval (PI) between 10% and 90% in kidney transplant patients. Semitransplant blue field represents PI between 10% and 90% in liver transplant patients.

SSE Analysis and Final Parameters

The results of the SSE analysis of 1000 simulation sets are presented in Table 3. The success rate of parameter estimation, including rounding errors with significant digits ≥2, was 72.4% (724 of 1000 simulation sets). Thus, the confidence intervals of the mean values and RSEs were calculated from 724 sets of successfully estimated parameters. The RSEs of the PK parameters ranged from 0.75% to 27.0%, and the sizes of interindividual variability of PK parameters ranged from 0.223 to 0.826. Among the major PK parameters, the apparent central volume (V/F) had the largest RSE value, and the absorption rate constant (ka) had the largest interindividual variability. Based on the 90% confidence intervals of the means and RSEs, the final parameter estimated from the SSE analysis was dependable and robust. The final estimated values of CL/F and V/F were close to those reported in other studies (21.2 L/h and 510 L, respectively).7,8 However, considering the distribution patterns of PK parameters (skewed distributions), medians were found to be more appropriate than means as representative values, and therefore, they were selected as the final parameters of the PK model.

TABLE 3. - Final Parameter Estimates Using SSE Method (n = 1000)
Parameter Description Median Mean (±95% CI) %RSE
Structural model
   Oral clearance (CL/F) before surgery (L/h) in kidney transplant patient 21.2 21.1 ± 0.0322 1.55
   Baseline CL/F after surgery in kidney transplant patients (L/h) 18.8 18.8 ± 0.0224 1.22
   Baseline CL/F after surgery in liver transplant patients (L/h) 11.4 11.1 ± 0.0159 1.31
   Age effect −0.382 −0.379 ± 0.00418 11.2
   POD effect in early period in kidney transplant patient 0.118 0.118 ± 0.000191 1.65
   POD effect in liver transplant patient 0.0582 0.0585 ± 0.000184 3.03
   POD effect in late period in kidney transplant patient −0.0100 −0.0101 ± 0.000139 14.0
   Apparent volume (L) 510 503 ± 1.95 4.38
   Body weight effect 1.30 1.37 ± 0.0363 27.0
   Absorption rate constant 3.10 3.21 ± 0.0427 16.4
   Lag time of first-order absorption 0.25 (fix) 0.25 (fix)
Interindividual variability
 ωCL/F Interindividual variability of CL/F 0.223 0.223 ± 0.00155 7.06
 ωV/F Interindividual variability of V/F 0.826 0.829 ± 0.00699 8.57
Correlation coefficients
 ρCL/F-V/F Correlation coefficients between CL/F and V/F 0.191 0.190 ± 0.00349 18.7
Residual error
 σadd, 1 * Additive error 0.0001 (fix) 0.0001 (fix)
 σprop, 1 * Proportional error 0.0653 0.0657 ± 0.000754 11.7
 σadd, 2 * Additive error −2.48 −2.44 ± 0.0418 17.4
 σadd, 3 * Additive error 2.86 2.85 ± 0.0021 0.75
*Number refers to the hospital where data were collected (1: Seoul National University Hospital, 2: Inje University Busan Paik Hospital, 3: Chonnam National University Hospital).
CI, confidential interval; CL, clearance; F, bioavailability; V, volume of distribution; WT, body weight.

Final Model Assessment

To assess the appropriateness of the final parameters estimated from the SSE analysis, pcVPC was performed with 13 bins for each institution per period. Because the sampling points were not spread evenly throughout the sampling period, the time interval of each bin was set according to the number of observations in each bin instead of the automated binning process provided by the PsN. The results of pcVPC of institutions based on PODs are presented in Figure 3. Based on the pcVPC results, the central trend and variability in the observed data could be reproduced appropriately by simulating the final PK model, indicating that the predictive performance of the final PK model was acceptable.

Prediction-corrected VPC (pcVPC) plot based on institution data and postoperative days (A: All data in early period, B: All data in late period, C: IUBPH data in early period, D: IUBPH in late period, E: CNUH data). Solid red line represents the median of prediction-corrected blood concentration (ng/mL), and semitransparent red field represents a simulation-based 95% confidence interval (CI) for the median. The observed 10% and 90% percentiles are presented with dashed red lines, and the 95% CIs for the corresponding simulation-based percentiles are presented as semitransparent blue fields. The prediction-corrected observed blood concentrations are presented as blue circles.


In this study, we developed a population PK model of tacrolimus that can play the role of a backbone for practical TDM procedures in Korean populations by reflecting the clinically relevant covariates and suggesting the feasible magnitude of PK parameters in the target population. The final estimated values of CL/F and V/F (21.2 L/h and 510 L, respectively) were close to those reported in other studies.7,8 Considering previous reports,4,7,8,14–16 the population values and variabilities of the parameters were not overestimated, and the residual error was within the known intrapatient variability of tacrolimus and conventional precision range of analysis. The final model appropriately described the observed data with an acceptable prediction interval.

We managed to overcome the limitations owing to the nature of TDM data in PK model development and demonstrated the possibility of using the SSE method as a valuable tool to integrate the published models with acquired data, especially when the raw data of published models could not be obtained. SSE can reduce the potential bias of a single simulated data set produced by the errors of generated random number. After summarizing the outcome from repeated estimation results, the variability of those errors can be referred as the precision of the parameters while the central tendency remains.17 Multiple studies have shown that stochastic modeling can reduce the errors and provide a better estimate than the deterministic method because this rigorous statistical approach can decompose the error into variability and system errors.17–21

As the present model was developed using real-world data, several clinically meaningful covariates were identified during the process. POD, age, and transplant organ (liver and kidney) were found to have a significant effect on CL/F, whereas only weight was associated with V/F. During the early period, CL/F gradually increased with time after immediate reduction owing to transplantation surgery, regardless of the transplant organ type. This finding is consistent with that of previous reports.7,22,23 The immediate decrease in CL/F after transplantation may be due to the decreased CYP3A activity from hyperacute postsurgical changes, such as acute inflammation and perfusion injury,24–26 and also by increased bioavailability because of low gastrointestinal motility and decreased gut metabolism.7,27 An increase in CL/F during the early period could be explained by several factors, including the recovery of gastrointestinal motility and CYP3A activity, coadministration of corticosteroids, and postoperative hypermetabolism.25 Corticosteroids are known to induce both CYP3A and P-glycoprotein, responsible for the metabolism and absorption of tacrolimus, respectively27; therefore, they can increase CL and decrease F of tacrolimus.28–31 However, the CL/F of liver transplant patients was mostly lower than that of kidney transplant patients during the early period, probably because it takes up to 6 months for the liver, which is responsible for the elimination of tacrolimus from the body, to recover to its normal functions after transplantation.23,32

After the early period, CL/F decreased with time and reached a plateau after 1 year in kidney transplant patients. This was related to a reduction in corticosteroid administration and multiple other factors including a reduction in postoperative hypermetabolism, ischemia–reperfusion injury, and increased hematocrit levels as the kidney function recovered.8,29,33 When tacrolimus binds to hematocrit extensively, whole-blood clearance decreases with increasing hematocrit levels.34–36 However, unlike in kidney transplant patients, a decrease in CL/F was not observed in liver transplant patients during the late period. Instead, it plateaued after a gradual increase, possibly due to the recovery and stabilization of liver function. The plateau level values of CL/F in both kidney and liver transplant patients were similar.

Patient's age was also identified as a major factor affecting CL/F. A possible explanation for the decrease in CL/F owing to aging may include age-related differences in P-glycoprotein and CYP3A enzyme activities as organ functions deteriorate in old age.37,38 Furthermore, in the case of V/F, the body weight of the patient was another significant covariate, and this finding was also consistent with that of previous reports.8


A tacrolimus PK model that can describe published PK models and newly collected data from the Korean population was developed using the SSE method. Although there are limitations in the model development owing to the nature of TDM data, the SSE method was found useful for obtaining the complete information from TDM data by integrating published PK models while maintaining the variability of the model.


The authors acknowledge the financial support of the Catholic Medical Center Research Foundation made in the program year of 2018.


1. Venkataramanan R, Swaminathan A, Prasad T, et al. Clinical pharmacokinetics of tacrolimus. Clin Pharmacokinet. 1995;29:404–430.
2. Kershner RP, Fitzsimmons WE. Relationship of FK506 whole blood concentrations and efficacy and toxicity after liver and kidney transplantation. Transplantation. 1996;62:920–926.
3. Holt DW, Armstrong VW, Griesmacher A, et al. International federation of clinical chemistry/international association of therapeutic drug monitoring and clinical toxicology working group on immunosuppressive drug monitoring. Ther Drug Monit. 2002;24:59–67.
4. Staatz CE, Willis C, Taylor PJ, et al. Population pharmacokinetics of tacrolimus in adult kidney transplant recipients. Clin Pharmacol Ther. 2002;72:660–669.
5. Burton M, Shaw L, Schentag J, et al. Applied Pharmacokinetics and Pharmacodynamics. Philadelphia: Lippincott Williams & Wilkins; 2006.
6. Allegaert K, Flint R, Smits A. Pharmacokinetic modelling and Bayesian estimation-assisted decision tools to optimize vancomycin dosage in neonates: only one piece of the puzzle. Expert Opin Drug Metab Toxicol. 2019;15:735–749.
7. Han N, Ha S, Yun HY, et al. Population pharmacokinetic-pharmacogenetic model of tacrolimus in the early period after kidney transplantation. Basic Clin Pharmacol Toxicol. 2014;114:400–406.
8. Han N, Yun HY, Hong JY, et al. Prediction of the tacrolimus population pharmacokinetic parameters according to CYP3A5 genotype and clinical factors using NONMEM in adult kidney transplant recipients. Eur J Clin Pharmacol. 2013;69:53–63.
9. Daniel JS, Peter B. MissForest—non-parametric missing value imputation for mixed-type data. Bioinformatics. 2012;28:112–118.
10. Jonsson EN, Karlsson MO. Xpose--an S-PLUS based population pharmacokinetic/pharmacodynamic model building aid for NONMEM. Comput Methods Programs Biomed. 1999;58:51–64.
11. Keizer RJ, Karlsson MO, Hooker A. Modeling and simulation workbench for NONMEM: tutorial on pirana, PsN, and Xpose. CPT Pharmacometrics Syst Pharmacol. 2013;2:50.
12. Bergstrand M, Hooker AC, Wallin JE, et al. Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models. AAPS J. 2011;13:143–151.
13. Lindbom L, Pihlgren P, Jonsson EN. PsN-Toolkit--a collection of computer intensive statistical methods for non-linear mixed effect modeling using NONMEM. Comput Methods Programs Biomed. 2005;79:241–257.
14. Prado-Velasco M, Borobia A, Carcas-Sansuan A. Predictive engines based on pharmacokinetics modelling for tacrolimus personalized dosage in paediatric renal transplant patients. Sci Rep. 2020;10:1–18.
15. Shuker N, Van Gelder T, Hesselink DA. Intra-patient variability in tacrolimus exposure: causes, consequences for clinical management. Transpl Rev. 2015;29:78–84.
16. Leino AD, King EC, Jiang W, et al. Assessment of tacrolimus intrapatient variability in stable adherent transplant recipients: establishing baseline values. Am J Transpl. 2019;19:1410–1420.
17. Paláncz B, Stewart K, Homlok J, et al. Stochastic simulation and parameter estimation of the ICING model. IFAC-PapersOnLine. 2016;49:218–223.
18. Tornøe CW, Jacobsen JL, Pedersen O, et al. Grey-box modelling of pharmacokinetic/pharmacodynamic systems. J Pharmacokinet Pharmacodyn. 2004;31:401–417.
19. Duun-Henriksen AK, Schmidt S, Røge RM, et al. Model identification using stochastic differential equation grey-box models in diabetes. J Diabetes Sci Technol. 2013;7:431–440.
20. Donnet S, Samson A. A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models. Adv Drug Deliv Rev. 2013;65:929–939.
21. Leander J, Almquist J, Ahlström C, et al. Mixed effects modeling using stochastic differential equations: illustrated by pharmacokinetic data of nicotinic acid in obese Zucker rats. AAPS J. 2015;17:586–596.
22. Antignac M, Barrou B, Farinotti R, et al. Population pharmacokinetics and bioavailability of tacrolimus in kidney transplant patients. Br J Clin Pharmacol. 2007;64:750–757.
23. Antignac M, Hulot JS, Boleslawski E, et al. Population pharmacokinetics of tacrolimus in full liver transplant patients: modelling of the post-operative clearance. Eur J Clin Pharmacol. 2005;61:409–416.
24. McKindley DS, Hanes S, Boucher BA. Hepatic drug metabolism in critical illness. Pharmacotherapy. 1998;18:759–778.
25. Boucher BA, Wood GC, Swanson JM. Pharmacokinetic changes in critical illness. Crit Care Clin. 2006;22:255–271.
26. Haas CE, Kaufman DC, Jones CE, et al. Cytochrome P450 3A4 activity after surgical stress. Crit Care Med. 2003;31:1338–1346.
27. Tuteja S, Alloway RR, Johnson JA, et al. The effect of gut metabolism on tacrolimus bioavailability in renal transplant recipients. Transplantation. 2001;71:1303–1307.
28. Plosker GL, Foster RH. Tacrolimus: a further update of its pharmacology and therapeutic use in the management of organ transplantation. Drugs. 2000;59:323–389.
29. Undre NA, Schäfer A. Factors affecting the pharmacokinetics of tacrolimus in the first year after renal transplantation. Transpl Proc. 1998;30:1261–1263.
30. Hesselink DA, Ngyuen H, Wabbijn M, et al. Tacrolimus dose requirement in renal transplant recipients is significantly higher when used in combination with corticosteroids. Br J Clin Pharmacol. 2003;56:327–330.
31. van Duijnhoven EM, Boots JMM, Christiaans MHL, et al. Increase in tacrolimus trough levels after steroid withdrawal. Transpl Int. 2003;16:721–725.
32. Li M, Zhao Y, Humar A, et al. Pharmacokinetics of drugs in adult living donor liver transplant patients: regulatory factors and observations based on studies in animals and humans. Expert Opin Drug Metab Toxicol. 2016;12:231–243.
33. Vlahakos DV, Marathias KP, Agroyannis B, et al. Posttransplant erythrocytosis. Kidney Int. 2003;63:1187–1194.
34. Beysens AJ, Wijnen RM, Beuman GH, et al. 506: monitoring in plasma or in whole blood? Transpl Proc. 1991;23:2745–2747.
35. Schijvens AM, van Hesteren FHS, Cornelissen EAM, et al. The potential impact of hematocrit correction on evaluation of tacrolimus target exposure in pediatric kidney transplant patients. Pediatr Nephrol. 2019;34:507–515.
36. Staatz CE, Willis C, Taylor PJ, et al. Toward better outcomes with tacrolimus therapy: population pharmacokinetics and individualized dosage prediction in adult liver transplantation. Liver Transpl. 2003;9:130–137.
37. Tanaka E. In vivo age-related changes in hepatic drug-oxidizing capacity in humans. J Clin Pharm Ther. 1998;23:247–255.
38. Gupta S. P-glycoprotein expression and regulation. Age-related changes and potential effects on drug therapy. Drugs Aging. 1995;7:19–29.

therapeutic drug monitoring; Korean population; tacrolimus; pharmacokinetics; nonlinear mixed-effect modeling

Copyright © 2022 The Author(s). Published by Wolters Kluwer Health, Inc. on behalf of the International Association of Therapeutic Drug Monitoring and Clinical Toxicology.