Pharmacokinetic Model Based on Stochastic Simulation and Estimation for Therapeutic Drug Monitoring of Teicoplanin in Korean Neutropenic Hematopoietic Stem Cell Transplant Recipients

Introduction

Teicoplanin, a glycopeptide antibiotic, is extensively used for treating serious Gram-positive bacterial infections, particularly methicillin-resistant Staphylococcus aureus (MRSA) and enterococci.¹ Teicoplanin is favored in the clinical setting due to its strong therapeutic efficacy and comparatively lower toxicity profile relative to vancomycin, which makes it particularly valuable for patients requiring complex and prolonged treatments, such as those undergoing hematopoietic stem cell transplantation (HSCT).² HSCT patients are highly vulnerable to severe infections, not only because of their immunocompromised status but also due to the cumulative stress from intensive conditioning regimens. This vulnerability makes accurate antimicrobial therapy essential, with teicoplanin being particularly important for preventing life-threatening infections and supporting patient recovery.^3,4

Teicoplanin exhibits a high protein binding rate (90–95%) and is primarily excreted unchanged by the kidneys, with minimal hepatic metabolism.^5,6 In patients with severe infections, pharmacokinetic (PK) parameters may vary considerably due to pathophysiological changes such as altered renal clearance, hypoalbuminemia, and increased distribution volume.^7,8 As these variations complicate infection management, tailoring therapy through therapeutic drug monitoring (TDM) has become a recognized strategy in clinical practice. In Korea, TDM of teicoplanin is commonly performed by monitoring trough concentrations (C_min) to ensure they remain above a target threshold.^4,9 However, the correlation between C_min and overall exposure may not be reliable in patients with high PK variability, and Bayesian approaches using population PK models are increasingly applied to address this limitation.^4,10

For accurate Bayesian TDM, a population-specific PK model tailored to the target patient group is essential. Only a few teicoplanin PK models have been reported in Korean patients, but none have been developed or validated in the HSCT setting. Consequently, TDM in Korean HSCT patients often depends on models derived from other populations, which may compromise predictive accuracy and dosing precision. Despite this knowledge gap, conducting a prospective full-PK study in this population is challenging, given the high disease severity and the intensive blood sampling already required during hospitalization. As a practical alternative, an approach that leverages C_min data obtained from routine TDM, combined with stochastic simulation and estimation (SSE), provides a feasible way to generate reliable PK information. Developing such a model is essential to improve the reliability of TDM and enable individualized therapy in this high-risk clinical group.

The objective of this study was to develop a PK model specifically for Korean post-HSCT neutropenic patients. Because C_min data alone are insufficient to construct a reliable model, we supplemented them with simulated datasets generated from two recently published Korean teicoplanin PK models.^11,12 However, even with this hybrid dataset, a single estimation may not yield sufficiently stable parameter values. Therefore, we applied an SSE approach, in which repeated simulations and parameter re-estimations allowed us to obtain robust estimates. This strategy provided a reliable foundation for precision TDM of teicoplanin and has the potential to improve infection control and safety outcomes in this vulnerable patient population.

Materials and Methods

Patient Population and Observed Data

TDM data were obtained by retrospectively collecting actual observations of patients treated at the Bone Marrow Transplantation Center (now separated into Catholic Hematology Hospital) at The Catholic University of Korea Seoul St. Mary’s Hospital from 2015 to 2017. The study population was patients aged 18 years and older who were diagnosed with hematologic malignancies, underwent HSCT, were hospitalized for infections, and received teicoplanin. Teicoplanin treatment was selectively administered in the presence of (i) positive gram-positive bacterial cultures, (ii) severe sepsis or shock with pending blood culture results, (iii) history of MRSA infection or colonization, (iv) skin and soft tissue infections, and (v) suspected catheter-related infections. Patients were excluded if it was determined that the PK characteristics of teicoplanin were not representative of the general characteristics of this patient group (eg, receiving massive fluid therapy, renal replacement therapy). No formal prospective sample size calculation was performed; rather, the study included all eligible patients who underwent TDM during the study period, thereby ensuring maximal utilization of available clinical data.

The standard dosing regimen consisted of three loading doses of 400 mg intravenously at 12-hr intervals, followed by a maintenance dose of 400 mg once daily. However, the actual dosing amount was determined at the discretion of physicians, and the standard regimen was not always applied. TDM blood samples were collected approximately one hour before a scheduled dose administration (23 hrs after the last administered dose during the maintenance period), following the general recommendations in Korea as described by Kim et al¹³ after maintaining a stable dosing regimen for at least 48–96 hours. However, minor variations occurred due to practical constraints in routine clinical practice. In this regard, all actual dosing history and sampling time information were obtained and reflected in the analysis. Samples were analyzed with ultra-high performance liquid chromatography tandem mass spectrometry (UHPLC-MS/MS) in a prior study.¹⁴ The calibration range was 3.9–52.9 mg/L. Separately, the lower limit of quantification was validated at 0.72 mg/L, demonstrating the analytical sensitivity of the method and confirming that it can accurately quantify teicoplanin concentrations even below the routine calibration range.¹⁴ Demographic data, including weight, height, sex, age, creatinine, and albumin concentrations, were also captured. Estimated glomerular filtration rate (eGFR) was calculated using the CKD-EPI (Chronic Kidney Disease Epidemiology Collaboration) formula.

This research was reviewed and approved by the institutional review board (IRB) of The Catholic University of Korea Seoul St. Mary’s Hospital (KC17SESI0654). The dataset was obtained retrospectively from electronic medical records, during which all subject-identifiable information was removed. As such, the study constituted a retrospective study using anonymized data and the requirement for informed consent was waived by the IRB. This study was conducted in accordance with the Declaration of Helsinki.

Simulated Dataset

Only two previous studies have reported teicoplanin PK models developed specifically for the Korean population.^11,12 The first model, a two-compartment model, was derived from a study involving 15 critically ill elderly Korean patients aged 60 years or older. Trough and peak concentrations were measured after multiple doses of teicoplanin, with the dosing intervals adjusted according to renal function.¹¹ The second model, a three-compartment model, was developed based on a full-PK study (8 PK samples taken post-infusion) in 12 healthy Korean individuals who received a single 200 mg dose of teicoplanin.¹² Although these populations do not exactly match the patients in the present study, they were considered suitable for generating simulated data for two reasons: i) the simulated dataset was restricted to uninformative timepoints and used only as backbone data for model development, applied consistently across all cases; and ii) the population PK parameter values and covariates were derived from Korean subjects, rather than from disease-specific conditions.

Simulations were conducted to obtain backbone data for model building based on the PK properties of these existing models. Covariates used for calculating PK parameter values included body weight, gender, age, serum creatinine, and eGFR, consistent with those adopted in the published models. To reflect realistic patient characteristics, multivariate normal distributions of these covariates were generated based on the distributions reported in the published Korean PK models. For each simulation run, random sampling was performed to create virtual subjects, ensuring that the number of simulated individuals matched the intended population size.

To keep the resulting model as dependent as possible on the observed TDM data, the number of simulated data points was minimized. Specifically, six early post-dose timepoints (0.3, 1, 1.5, 2.5, 4, and 6 hr) were selected to capture the ^initial distribution and early disposition phases of teicoplanin, which cannot be adequately characterized by trough-only (C_min) data. This strategy ensured that essential PK processes were represented while minimizing the influence of simulated data, despite teicoplanin’s long half-life. Furthermore, 30 subjects were simulated from each model (60 subjects in total), which is the minimum sample size required by the central limit theorem to estimate between-subject variability (ω²). Each dataset therefore contained 360 simulated concentrations (6 timepoints × 60 subjects). In this process, between-subject variability reported in the literature was incorporated into the PK parameters, and residual error was reflected according to the original models.

All random sampling and estimation procedures were performed using R (version 4.4.0) with the “tmvtnorm” package (version 1.6.0), and the simulated datasets were subsequently generated in NONMEM (version 7.5, Icon Development Solutions, Ellicott City, MD).

Base Model Building

The first combined dataset of observed and simulated data was used to develop the base model. Population PK analysis was performed using non-linear mixed-effects modeling with NONMEM. The first-order conditional estimation with interaction (FOCE-I) method was used throughout the modeling process. Two- and three-compartment models with first-order elimination were assessed, as these models were used to generate the simulated dataset. The individual variability of each PK parameter was described using a lognormal distribution.

P_i represents the individual parameter, θ_pop represents the typical parameter value for the population, and η_i is the between-subject variability (BSV), which follows a normal distribution with a mean of zero and a variance of ω². Possible correlations between the individual variabilities were evaluated using OMEGA BLOCK option. Proportional, additive, and combined error models (ε) were assessed for residual variability, which follows a normal distribution with a mean of zero and a variance of σ². To identify the optimal structural model, the objective function value (OFV) was used to compare nested models, and the Akaike information criterion (AIC) was used to compare non-nested models. Additional criteria such as goodness-of-fit plots, parameter stability, and parsimony were also comprehensively considered in the final model selection process. Covariate analysis was conducted to evaluate the impact of the potential covariates on the model. Prior to covariate selection, covariate screening was performed using generalized additive models (GAM) and graphical analysis to identify potential covariates. Stepwise covariate modeling in Perl-speaks-NONMEM (PsN, version 5.3.1)¹⁵ was then used, which included both forward selection (p<0.05, ΔOFV > 3.84) and backward elimination (p<0.01, ΔOFV > 6.63). Tested covariates included height, weight, sex, albumin, serum creatinine level, and estimated creatinine clearance. Biologically plausible and physiologically relevant covariates were selected for the final base model to ensure clinical relevance and improve the predictive performance of the model.

Base Model Evaluation

The bootstrap method and visual predictive check (VPC) were used to evaluate the final base model. 1000 randomly resampled datasets were generated from the combined dataset. The medians and the 5th and 95th percentiles of 1000 parameter estimates were compared with the original base model parameter estimates. For the VPC, simulations were performed 1000 times, and the results were plotted, comparing the observed concentrations and 90% prediction interval of the simulated concentrations. Both bootstrapping and VPC were performed using PsN (version 5.3.1), NONMEM, and R (version 4.4.0).

Stochastic Simulation and Estimation (SSE) to Obtain Final Parameter Estimates

Using the base model, the process of 1) obtaining simulated data, 2) constructing the dataset by merging with observed data, and 3) estimating parameters through model fitting was repeated 1000 times. The FOCE-I method was used for all estimations. In each estimation output file, it was evaluated whether the minimization was successful and whether the parameter values were reliable, and the proportion of estimation runs that satisfied both conditions was reported. Using only the estimation results that satisfied both conditions, the mean and 90% confidence intervals of the parameter estimates were presented along with the relative standard error (RSE).

Predictive Performance Evaluation for the Final Model

To assess the final model and estimated parameters from SSE, a prediction-corrected visual predictive check (pcVPC) was performed.¹⁶ Traditional VPC can be affected by binning across a large variability in dose and/or influential covariates, leading to misinterpretation of the data. Also, it can be inaccurate when used with adaptive designs, such as dose adjustments. Therefore, the presence of different dosing regimens, sampling times, and covariate values in the combined dataset makes pcVPC the appropriate method to evaluate instead of traditional VPC. pcVPC enables correction of the variability coming from binning across the independent variables such as dose, time, and covariate by normalizing the dependent variables based on the typical population prediction. The formula for population prediction correction is as follows:

pcY_ij is the prediction-corrected observation or prediction for the ith individual and jth time point, Y_ij is the observation or prediction for the ith individual and jth time point, PRED_bin is the median of typical population predictions for the specific bin of independent variables, and PRED_ij is the typical population prediction for the ith individual and jth time point. Only observed data were utilized in the pcVPC procedure, not the actual dosing history. Because of the unevenly distributed sampling time points throughout the sampling period, the binning was manually set to 32 bins, considering the dosing and observed time.

Results

Observed Dataset

A total of 417 individuals (222 males and 195 females) who underwent teicoplanin TDM were initially included in the analysis. After excluding 12 patients with missing covariate values and one patient with an erroneously recorded concentration, the final observed dataset comprised 405 patients with 568 valid observations. Sampling times ranged from 19 to 2879 hours after the first dose, yielding 581 teicoplanin concentration measurements. Among these, data records from patients with a regimen maintenance period of less than 48 hours accounted for 7.6%, and the latest sample was collected 47 hours after the final dose. The baseline characteristics of these patients are summarized in Table 1.

Table 1 Demographics of Study Patients

Base Model and Evaluation Outcomes

Because the simulated dataset inherently contained early timepoints that required at least a two-compartment structure, both two- and three-compartment models were evaluated during base model development. The three-compartment model produced superior results in terms of diagnostic criteria, including a lower OFV as well as better normality and homoscedasticity. However, not all patients were adequately described by the three-compartment structure, and a substantial proportion could be sufficiently explained by a two-compartment model alone. In the SSE procedure, where subject-level random sampling is repeated, the inclusion of such individuals increased the likelihood of convergence failure when using the three-compartment model. Therefore, despite acknowledging some degree of structural misspecification, the two-compartment model was selected as the final base structure because it provided more stable and reliable parameter estimation across repeated simulations (Figure 1). Larger weighted residuals were observed at higher concentrations, but this was considered unavoidable given the limited number of such samples and the variability inherent in routine clinical sampling.

Figure 1 Basic goodness-of-fit plots of the base population PK model. Each panel displays four standard diagnostic plots. Top left (A) is a scatterplot of observed concentrations versus population predictions (PRED). Top right (B) is a scatterplot of observed concentrations versus individual predictions (IPRED). Bottom left (C) is a scatterplot of individual weighted residuals (IWRES) versus IPRED. Bottom right (D) is a scatterplot of conditional weighted residuals (CWRES) versus time. The circles indicate the observed or predicted plasma concentration of teicoplanin (mg/L). The black line and red line indicate the line of identity, and locally weighted scatterplot smoothing line, respectively.

Among the two-compartment models, the model with BSV on CL and V2 had the lowest OFV and was selected as the base structural model. As a result of forward selection and backward elimination in stepwise covariate modeling, eGFR and albumin were identified as potential covariates for CL. Considering biologically plausible and physiologically relevant covariates, eGFR and albumin were selected as final covariates. A comparison of the fit statistics (minimized OFV) by model structure is presented in Table 2.

Table 2 Comparison of Structural Model

Bootstrap evaluation showed that most of the estimations were completed successfully, and no particular problems were found with the point estimates and interval estimates of the parameters. In addition, VPC results showed that the model reproduces the trend of the dataset well enough, although there is some variation depending on the density of the data at each time point (Figure 2). This confirmed that the final model satisfied the basic requirements for a base model for SSE analyses. The final base model PK parameter estimates and bootstrap results are summarized in Table 3.

Table 3 Parameter Estimates of Base Pharmacokinetic Model with Bootstrap Results

Figure 2 VPC plots for the base model. (A) Simulated data (B) Observed TDM data. All time points were converted to time after dose. Simulated data and observed trough concentrations are depicted as gray dots. The figure presents the VPC with the prediction intervals for the 5th, 50th, and 95th percentiles of simulations from bottom to top. The solid lines correspond to 5th, 50th, and 95th percentiles. The red shaded areas correspond to the 95th confidence interval of the 5th and 95th percentile. The blue shaded area corresponds to the 95th confidence interval of the 50th percentile.

Final Model From SSE Procedures

Out of 1000 SSE runs, the total number of cases where the estimation was not done properly (minimization failure, parameter near boundary, or rounding error) was 11.9% (88.1% of successful runs). The bias of the mean parameter estimates, calculated from the successful runs relative to the final estimates of the base model, was not significant. The RSEs ranged from 0.522% to 8.88% depending on the parameters. In particular, the RSEs for CL and V2 were 2.09% and 6.81%, respectively, which were deemed to have sufficient reliability to be used for post hoc parameter (empirical Bayesian estimate) estimation. The final PK parameter estimates are summarized in Table 4. The plot of pcVPC results (Figure 3) using the final model parameters obtained from SSE demonstrates that the PK model is capable of accurately reproducing the observed trend, reflecting the actual dosing history and sampling time. Notably, the predictive performance was optimal within the range of 720 hours from the first dose, which coincides with the primary period for TDM procedures. However, the prediction interval exhibited some inflation in the later time ranges, where the data became increasingly sparse.

Table 4 Parameter Estimates of the Final Pharmacokinetic Model Using SSE (N = 1000)

Figure 3 pcVPC plots for the final model. (A) Data up to 720 hours. (B) Data after 720 hours. Solid red line represents the median of prediction corrected concentration (mg/L), and semitransparent red field represents a simulation-based 95% confidence interval for the median. The observed 10% and 90% percentiles are presented with blue dashed lines, and the 95% confidence intervals for the corresponding simulation-based percentiles are shown as semitransparent blue fields. Black dots are the prediction-corrected observed concentrations.

Discussion

In this study, we developed a population PK model of teicoplanin specifically for Korean HSCT patients with neutropenia, providing the first population-specific parameter estimates in this setting. The model was further refined by identifying eGFR and albumin as covariates of clearance, which strengthens its clinical relevance for individualized dosing. By applying the SSE approach, we overcame the limitation of relying solely on C_min data and achieved stable parameter estimation. The adequacy of this methodology was demonstrated through the satisfactory predictive performance of the final two-compartment model for observed concentrations across various dosing histories. Taken together, these results indicate that the model can improve the reliability of TDM in Korean HSCT patients, particularly within the clinically relevant timeframes in which teicoplanin monitoring is usually performed.

The two published Korean studies that informed our simulations employed two- and three-compartment models, respectively, explaining some of the parameter similarities with our results. V1 was estimated at 4.71 L, close to values reported previously (4.79 and 3.96 L), indicating that the simulated dataset adequately fulfilled its intended role in model development. By contrast, clearance was estimated at 1.26 L/h, nearly double the values reported in the earlier studies (0.63 and 0.69 L/h). Because our analysis assumed the same distribution characteristics, the larger clearance estimate indicates that HSCT recipients had lower C_min than prior populations. When a model derived from a reference population is applied to this study population, observed concentrations are therefore likely to be consistently lower than the predicted values. As a result, Bayesian TDM based on such external models may yield theoretical C_min values that are higher than those actually observed. This in turn reduces the likelihood of appropriate upward dose adjustment. Our model instead provides a basis for Bayesian TDM that mitigates this systematic bias and supports more reliable individualized dosing in HSCT patients.

This finding illustrates the significance of population specificity in PK modeling. Even within the Korean population, models developed in patient groups without neutropenia or HSCT are unlikely to capture the PK characteristics of our study cohort. Likewise, although models have been reported in other ethnic groups with hematologic malignancies and neutropenia, the documented inter-ethnic variability in glycopeptide PK raises serious concerns about their applicability to Korean patients.^17,18 Applying such externally developed models without rigorous validation carries a substantial risk of systematic bias and clinically inappropriate predictions, particularly for clearance estimates and covariate–parameter relationships.¹⁹ Therefore, although the use of SSE with simulated data represents only a provisional approach, the development of population-specific PK models in groups where dense PK sampling is rarely feasible due to clinical and ethical constraints can nonetheless be regarded as methodologically sound and contributes to the establishment of clinically relevant TDM strategies.

The covariates identified in our model were consistent with those reported in the literature and were deemed appropriate for representing the known PK properties of teicoplanin. The higher mean eGFR in our population (113 mL/min/1.73 m² compared with 64 and 103 in the literature) may have partially contributed to the greater clearance estimate, plausibly reflecting the reduced muscle mass and lower serum creatinine levels commonly observed in HSCT patients.²⁰ The between-subject variability was larger than that reported previously but remained within a reasonable range, which can be explained by the heterogeneous nature of the HSCT population. These results indicate that the final model was not solely dependent on such data and adequately reflected the characteristics of the TDM patient cohort.

The present study was partially based on existing PK models derived from a Korean population, which may not fully align with the specific characteristics of post-HSCT neutropenic patients targeted in this analysis. This is because reliance on external datasets introduces potential variability that may not be fully representative of this population. In addition, an increased distribution volume may have contributed to the lower trough concentrations observed, but due to the limited information these effects were primarily reflected in clearance, which may have led to its overestimation. Taken together, these factors highlight a degree of uncertainty regarding the model’s ability to fully capture the distinctive PK profile of teicoplanin in this patient group. Future studies incorporating sparse sampling and real-world data on covariates (eg, CKD-EPI-based eGFR and serum albumin levels) from post-HSCT neutropenic patients will be essential to validate the model, improve its accuracy, and enable further optimization of individualized dosing strategies.

Conclusion

This study developed the first population PK model of teicoplanin tailored to Korean HSCT patients with neutropenia. The model provided population-specific parameter and variability estimates and identified eGFR and albumin as clinically relevant covariates of clearance. By applying the SSE approach, we were able to overcome the limitation of trough-only data and obtain stable parameter estimates. The final two-compartment model demonstrated satisfactory predictive performance within the clinically relevant timeframe for TDM, offering a reliable basis for individualized dosing in this vulnerable group. Although our approach relied partly on simulated data, the results underscore the value of developing population-specific models in populations where dense PK sampling is rarely feasible. Future studies incorporating real-world covariates and sparse sampling will be essential to validate and refine this model to further support safe and effective teicoplanin therapy in HSCT patients.