Population pharmacokinetic models and Bayesian forecasting can account for inter-patient variability in pharmacokinetics, and so hold particular promise for addressing the inherent complexities of vancomycin dosing in pediatric subpopulations with known alterations in vancomycin pharmacokinetics [5]. Here, we evaluated the predictive performance of two general and several specialized popPK models for vancomycin in large multi-site samples of 371 pediatric oncology patients (Study 1) and 219 pediatric CVICU patients (Study 2). We also refit the best performing models to evaluate if predictive performance could be improved further. A strength of our approach is that predictive performance was evaluated iteratively, where each serum vancomycin level was predicted using Bayesian estimates informed by previous levels. This approach mimics clinical decision making with MIPD in practice—where model predictions are used to inform dosing decisions without knowing future lab results—providing the best translation of model fit-for-purpose for clinical performance [12].

In pediatric CVICU patients we found that models developed on this patient population performed better than general population models, with the refit Shimamoto 2024 model performing best (Fig. 3). In contrast, we found that for pediatric oncology patients, a general popPK model (Colin 2019) outperformed purpose-built models (Fig. 2). Moreover, re-estimating the effect of cancer status on clearance for the Colin 2019 (Onc) model resulted in an estimate statistically indistinguishable from zero, suggesting that cancer status was not predictive of patient PK in this data set (Table 4). Interestingly, both general popPK models (Colin 2019, Le 2014) also performed better in the oncology patients than in matched non-oncology patients (Fig. S5). Together, these results suggest that both well-specified general models and specialized models can achieve suitable clinical performance in distinct pediatric subpopulations, and that this assessment must be made on a case-by-case basis.

Pediatric CVICU patients in our data set tended to be young, with an interquartile post-menstrual age (PMA) of 61–310 weeks and an interquartile chronological age of 0.4–5.2 years (Table 1). Renal function matures rapidly in these early months of life—reaching 90% of adult glomerular filtration rate (GFR) by age 1 and 98% by age 2 in full-term neonates [40]—and varies as a function of both gestational age (GA) at birth and postnatal age, with birth-related differences in GFR between pre-term and full-term neonates with the same PMA becoming negligible beyond 65 weeks [41]. Of the three specialized CVICU models in Study 2, one model (Kamp 2024) lacked a term for GA (Table 2), despite more than a quarter of their patients being under the age of 1 month, with a maximum age of 17 years (Table 3). This model instead relied only on body size (weight) and creatinine clearance (a derived covariate that incorporates age) to scale PK parameters, and it is possible that the lack of maturation function resulted in its relatively poor performance in our data set. The other two CVICU models incorporated PMA as a covariate, as well as weight and creatinine clearance, perhaps allowing them to better capture changes in drug clearance as a function of maturation, serum creatinine, and body size. The oncology patient population was generally older, with 75% of patients having a PMA of 220 weeks or more; however, here too, the best performing models [19, 36] included PMA as a covariate.

All models that included a maturation effect incorporated this effect only on clearance; however, body water composition also varies in newborns and infants compared to older children [42]. Since vancomycin is a hydrophilic drug, it is possible that VD also varies with age. However, we did not observe a relationship between age and volume, and it is possible this effect size is small relative to other sources of variability.

Notably, GA was missing and imputed to a typical full-term birth value of 40 weeks for 362 (97.6%) oncology patients in Study 1, and 174 (79.5%) CVICU patients in Study 2. Previous work has shown that imputing a fixed value for GA can be a reasonable approach when the information needed for more sophisticated methods is unavailable, with the caveat that the true GA may be substantially different from the imputed GA when the same fixed value is used for full-term or pre-term births [43]. For patients with missing GA, we also lacked information on full-term or pre-term birth status, and thus we could only use a single fixed value typical of full-term births; therefore, PMA may have been erroneously estimated for an unknown number of patients in either study. However, in a more general infant population, we found that imputing a GA of 40 weeks had minimal impact on vancomycin model performance for patients over 37 weeks of PMA [5], which was true of all patients in both studies (Table 1), and the present results suggest that selecting a popPK model that accounts for maturation may be preferable to those that do not, even when GA is unknown. In cases where other information such as birth weight or term status are available, more accurate imputation strategies can be used for missing GA [43].

It should also not be taken for granted that a model tailored to a patient population will outperform a literature model. For the CVICU population, re-estimation of model parameters improved predictive performance, while no improvement was observed after re-estimation of the best-performing model for oncology patients.

In the case of the Colin 2019 (Onc) model, the cancer status covariate was a study-specific parameter that came from adult oncology patients who had significantly higher clearance relative to patients from other studies in their meta-analysis [36, 44]. Although previous studies have reported higher clearance in adult [45, 46] and pediatric [31, 47] oncology patients, we are unaware of any published popPK models that have directly estimated the magnitude of additional clearance attributable to cancer status in this population. Instead, these reported differences have been based on either mean differences between oncology and non-oncology patients [45,46,47] or nominal comparisons between studies [31].

Pediatric popPK studies that have tested the effect of cancer status [27], type of hematological disease [30, 31], or primary diagnosis [19] on vancomycin clearance, all discarded these covariates from their final model, suggesting that cancer itself may not be an independent risk factor for augmented renal clearance or altered vancomycin PK. This assertion is supported by Hirai and colleagues [48], who found that febrile neutropenia—but not cancer, other critical illnesses, or surgery—significantly increased the risk of augmented renal clearance in pediatric patients, indirectly influencing vancomycin clearance. Therefore, it is likely that the reported effect of cancer status on vancomycin clearance in the published Colin 2019 (Onc) model was either age-specific or a batch effect (i.e., related to that particular data set and not the underlying biology). Future research should explore covariates that might better explain augmented renal clearance in pediatric oncology patients, such as body temperature [19] or neutropenia status.

In pediatric CVICU patients, the general-purpose Colin 2019 model had adequate predictive performance but was outperformed by all three specialized models a posteriori on all metrics except bias, and this slight edge was eliminated after refitting the Shimamoto 2024 model (Fig. 3). It is interesting that the refit Shimamoto 2024 model, a two-compartment model, outperformed the refit Moffett 2019 model, given that the vast majority of published pediatric vancomycin models use one compartment [2, 10, 11]. However, in contrast to the general pediatric population and other subpopulations, two-compartment models are more often the best structural models in critically ill children [10], and a previous PK analysis of vancomycin in infants undergoing open-heart surgery with cardiopulmonary bypass also found that a two-compartment model performed best in their sample [49]. Together, these findings suggest that an appropriately specified two-compartment vancomycin model may be beneficial in pediatric CVICU patients, as it is in other subpopulations.

Additionally, as in previous studies, we found that MAP Bayesian estimates of individual PK parameters improved predictive performance and attenuated the differences in predictive performance between models [5]. Therefore, clinicians could collect earlier drug levels to improve dose optimization in cases where a pharmacokinetic model appears to be a poor fit for a patient, or uncertainty in covariate values is suspected. For models that continue to fit poorly, selecting another model that better reflects patient covariates or subpopulation, or using clinical judgment to adjust dosing from model predictions, may improve dose optimization for an individual patient [13].

Because both our studies used routine clinical care data, our approach has several limitations related to data quality and availability. These limitations stem from the scope of retrospective data not originally collected for research purposes, which may have a higher risk of incorrect or missing data that can introduce bias. Where possible, we reduced the risk of this bias by maintaining strict inclusion criteria and exclusion criteria for the data used in both studies; however, we were limited in our ability to address bias associated with several forms of uncaptured data.

First, inclusion criteria for either study required manual tagging of patients as an oncology patient (Study 1) or a CVICU patient (Study 2) by the clinician during treatment. This is an optional feature in the InsightRX Nova CDS software tool that uses a free text field, meaning that cancer types were not always known and that patients in the “general” population may also have been oncology or CVICU patients but untagged. Second, body temperature was not recorded for any oncology patients in our sample and had to be imputed/assumed based on other information—specifically, whether or not patients were diagnosed with febrile neutropenia. This diagnosis was also entered using the optional manual tagging feature in the InsightRX Nova CDS software tool, therefore, it is possible that there were additional patients in our sample with febrile neutropenia who were not tagged. Indeed, the improved a priori predictive bias for the Shimamoto 2021 model when assuming all patients had a fever suggests that a positive neutropenia diagnosis may not have been tagged in many cases (Fig. S15). Third, we did not have access to the assays used to measure serum vancomycin concentrations, as this information is not recorded in the InsightRX Nova CDS software tool; although inter-assay and inter-batch variability in vancomycin assays may impact the accuracy of serum vancomycin concentration measurements, these are relatively limited compared to the overall residual error magnitude, so we did not expect a major influence of this factor. Fourth, one oncology-supporting model identified in our literature review (Guilhaumou 2016) was excluded from evaluation, given that it included covariates unavailable in the InsightRX Nova CDS software tool. To the best of our knowledge, both studies otherwise included all relevant oncology-supporting and CVICU-supporting models. Finally, our study only included patients whose data were entered into the InsightRX Nova CDS software tool. Individual sites may have had protocols that dictate alternative approaches for vancomycin therapy management for certain subsets of patients, representing an unknown type of bias.

A further limitation of both our studies is that we excluded certain patient categories from analysis (patients on HD, CRRT, ECMO, and VAD or who were tagged as NICU or PICU patients), as these patient populations were considered pharmacokinetically distinct from our target cohorts of typical oncology and CVICU patients. Specifically, both renal replacement therapy and mechanical circulatory support can further increase the variability of vancomycin PK in critically ill patients [50,51,52,53], adding an additional layer of complexity to MIPD that may not be captured by conventional PK models [52, 54, 55]. Likewise, whereas typical oncology and CVICU patients have predictable PK alterations associated with their conditions [6,7,8], patients admitted to the NICU or PICU present with heterogeneous critical illnesses (e.g., sepsis, trauma, respiratory failure, prematurity, etc.) that introduce fundamentally different pathophysiological states. The varied hemodynamic profiles, organ perfusion patterns, and disease processes in these patients would introduce substantial heterogeneity that would limit the generalizability of our findings to the primary populations of interest. In our data set, ward location was determined based on clinical use of an optional tagging feature and applied to an entire treatment course. It is likely that patients tagged with both CVICU as well as a second ICU tag were moved from one ward to another, suggesting different patterns of care or comorbidities, or reflecting different disease progressions relative to patients tagged CVICU only.

The decision to exclude these patient categories largely reflected the exclusion criteria of the models considered in our comparisons: patients receiving extracorporeal organ support were excluded from the development population of all models, with the exception of the Kamp 2024 model that included a small number of patients on ECMO. General ICU patients were only included in the development populations of the Colin 2019, Hadi 2015, Le 2014, and Kamp 2024 models. Therefore, we caution against generalizing the results of our studies to these patient groups. Development and external evaluation of models that can adequately capture vancomycin PK for oncology and CVICU patients receiving extracorporeal organ support or presenting with concurrent critical illnesses is still needed, and it is unclear whether, for example, a general PK model for patients receiving renal replacement therapy would have better predictive performance than the best performing models identified here. Until then, we recommend performing intensified PK monitoring of vancomycin when applying MIPD to oncology and CVICU patients who are receiving extracorporeal organ support or are admitted to the NICU or PICU to prevent inadequate treatment outcomes, regardless of the chosen model.

Finally, without dense sampling, we were unable to assess the ability of the models to accurately estimate the true area under the concentration–time curve (AUC) experienced by a patient. Although this is the primary drug exposure metric of clinical interest in vancomycin [4], we have previously shown that the ability to accurately predict the next vancomycin serum level in sparsely sampled routine clinical data correlates well with the ability to adequately estimate AUC [14].