Optimizing Vancomycin Therapy in Critically Ill Children: A Population Pharmacokinetics Study to Inform Vancomycin Area under the Curve Estimation Using Novel Biomarkers

Area under the curve (AUC)-directed vancomycin therapy is recommended, but Bayesian AUC estimation in critically ill children is difficult due to inadequate methods for estimating kidney function. We prospectively enrolled 50 critically ill children receiving IV vancomycin for suspected infection and divided them into model training (n = 30) and testing (n = 20) groups. We performed nonparametric population PK modeling in the training group using Pmetrics, evaluating novel urinary and plasma kidney biomarkers as covariates on vancomycin clearance. In this group, a two-compartment model best described the data. During covariate testing, cystatin C-based estimated glomerular filtration rate (eGFR) and urinary neutrophil gelatinase-associated lipocalin (NGAL; full model) improved model likelihood when included as covariates on clearance. We then used multiple-model optimization to define the optimal sampling times to estimate AUC24 for each subject in the model testing group and compared the Bayesian posterior AUC24 to AUC24 calculated using noncompartmental analysis from all measured concentrations for each subject. Our full model provided accurate and precise estimates of vancomycin AUC (bias 2.3%, imprecision 6.2%). However, AUC prediction was similar when using reduced models with only cystatin C-based eGFR (bias 1.8%, imprecision 7.0%) or creatinine-based eGFR (bias −2.4%, imprecision 6.2%) as covariates on clearance. All three model(s) facilitated accurate and precise estimation of vancomycin AUC in critically ill children.


Introduction
Vancomycin is the drug of choice for treatment of serious Gram-positive infections in children and is one of the most frequently administered drugs in the pediatric intensive care unit [1]. Its efficacy and toxicity are most closely related to an individual's 24 h area under the curve (AUC 24 ) [2][3][4][5][6]. Traditionally, vancomycin therapeutic drug monitoring (TDM) relied upon measurement of trough concentrations (C min ), which were used as a surrogate for AUC 24 [7]. However, with the continued maturation of Bayesian dosing software programs, AUC-based dosing using population models and limited sampling is becoming more routine [8]. Unfortunately, Bayesian AUC estimation from a single trough measurement can be inaccurate, and validated vancomycin population pharmacokinetic (popPK) models to inform Bayesian dosing and AUC estimation in critically ill children are lacking.

Study Population
We performed a prospective, observational study in the Pediatric Intensive Care Unit (PICU) at the Children's Hospital of Philadelphia (CHOP) from August 2018 to July 2021. Patients aged 1-17 years old receiving intermittent dosing of intravenous (IV) vancomycin for a suspected infection, defined as performance of a microbiological culture within 24 h of vancomycin initiation, were eligible for inclusion. Eligible patients were identified as soon after initiation of vancomycin as possible. Those receiving renal replacement therapy, plasmapheresis, or extracorporeal membrane oxygenation were ineligible. To be deemed evaluable, subjects had to have ≥3 PK samples collected, as well as ≥1 urine sample and 1 plasma sample collected for biomarker measurement (see below).
After enrollment was completed, evaluable subjects were divided into model training (n = 30) and testing (n = 20) groups. To perform this, we randomly selected 20 subjects who had ≥4 PK samples obtained to serve as the model testing group, while the remainder were assigned to the model training group. This approach was taken to facilitate AUC 24 estimation using noncompartmental methods within the model testing group, as described below.
The CHOP Institutional Review Board (IRB) approved the study protocol (IRB 18-014851, Approved: 20 March 2018) with a waiver of documented assent; verbal assent was obtained, as appropriate. Documented informed consent was obtained from the patient's parent(s)/legal guardian(s).

Dosing, PK Sampling, and Biomarker Measurement
Vancomycin was ordered for clinical care in all patients, with dosages and infusion rates determined by the clinical team. Typical initial dosages were 10-15 mg/kg/dose every 6-8 h, depending on age, weight, and estimated renal function. All decisions about dosing and duration of therapy were at the discretion of the clinical team. During the study, routine TDM included collection of trough (C min ) concentrations, with doses adjusted to achieve a goal of 5-15 µg/mL, as appropriate. In 2020, AUC-based TDM was implemented at our hospital, such that patients requiring >48 h of vancomycin therapy had two blood samples collected (at 1 h after the end of the infusion and at C min ). AUC was calculated using log-linear methods and doses were adjusted to achieve an AUC of 400-600 mg h/L.
For each participant, five vancomycin concentrations (PK sampling) were obtained after ≥4 vancomycin doses during a single dosing interval at the goal times shown in Table S1; since clinical care could interfere with the precise timing of sample collection, samples were accepted outside of these windows. Sampling took place within 48 h of enrollment for all participants. Samples could be collected prior to the fourth dose if a patient had impaired renal function such that he/she did not receive vancomycin at regular dosing intervals, so long as he/she had received ≥24 h of treatment; this applied to only one participant. All samples were collected via arterial catheter, peripheral venipuncture, or venous catheter, if not used for administration of vancomycin. In addition to PK sampling, the clinical team performed therapeutic drug monitoring (TDM) according to institutional standard practice. Results of TDM samples were recorded and included in our study if collected prior to PK sampling; individuals only participated through the time of PK sampling, so dosing information and TDM results that took place after PK sampling for the study were not recorded.
All vancomycin concentrations were measured by chemiluminescent microparticle immunoassay (Abbott, Abbott Park, IL, USA) in the CHOP Chemistry Laboratory; the lower limit of quantification (LLOQ) of this assay was 3.0 µg/mL. Results of both PK sampling and TDM were made available in the subject's electronic medical record.
After enrollment, urine and blood samples were obtained for biomarker measurement. Two urine samples were collected from subjects for measurement of the following biomarkers: NGAL, KIM-1, cystatin C, osteopontin, and creatinine. These samples were obtained starting the day of enrollment in the evening (3-7 pm) and morning (6-10 am) prior to PK sampling. Collection of urine samples outside of the above windows was permitted, based on the condition and clinical care needs of the patient. Urine samples were collected via indwelling urinary catheter, clean intermittent catheterization (if performed for clinical care), cotton balls, urine cup, or urine bag. In addition to prospective urine collection, we also identified and obtained any available residual urine samples that had been collected clinically within 24 h prior to initiation of IV vancomycin. These samples served as baseline biomarker measurements when available. A single blood sample was drawn within 24 h of PK sampling for measurement of plasma biomarkers (cystatin C and NGAL), as well as creatinine if not performed clinically. We similarly identified any available residual plasma samples that had been collected clinically up to 24 h preceding initiation of IV vancomycin to serve as baseline measurement of biomarkers.
KIM-1, osteopontin, NGAL (urine and plasma), and cystatin C (urine and plasma) were measured using Quantikine ® ELISA (R & D Systems, Inc., Minneapolis, MN, USA). Urine creatinine (uCr) was measured via two-point end enzymatic method on a Cobas ® system (Roche Diagnostics, Basel, Switzerland). These tests were performed in the CHOP Translational Core Laboratory. Plasma creatinine was measured by two-point rate spectrophotometric method (Vitros5600™ analyzer, Ortho Clinical Diagnostics, Raritan, NJ, USA) in the CHOP Chemistry Laboratory. Figure 1 displays a flow diagram of our approach to population PK model training and testing. Nonparametric population PK modeling was performed using the Pmetrics package (version 1.9.7; Laboratory of Applied Pharmacokinetics and Bioinformatics, Los Angeles, CA, USA) [28] for R (version 3.6.3; R Foundation for Statistical Computing, Vienna, Austria) [29] in RStudio (v1.2.5033; RStudio, Inc., Boston, MA, USA) [30]. Oneand two-compartment models were constructed using the nonparametric adaptive grid (NPAG) algorithm [31]. Parameters included V d (volume of the distribution) and CL (total body clearance) for one-compartment models and CL (clearance), Q (intercompartmental clearance), V1 (central volume), and V2 (peripheral volume) for two-compartment models. Clearance parameters were allometrically scaled for standardized weight to a power of 0.75 and volume parameters were scaled by standardized weight (power 1); weight was standardized by the median of the subjects' weights (27 kg). The weighting function on observations was 1/(gamma * SD 2 ), where SD (standard deviation) was a combined additive and multiplicative function of the assay imprecision as a polynomial equation: SD = C 0 + [C 1 × observed concentration]; C 0 had a value of 1.5 (half the LLOQ) and C 1 a value of 0.1, i.e., an assay with 10% coefficient of variation (CV%). Gamma was initially set to 1 and fitted to estimate the residual model error. Following determination of the base model structure, covariate model selection was then performed in a two-stage process. Since vancomycin is renally eliminated, we first sought to identify the best renal function marker to include as a covariate on CL. This included plasma biomarkers (creatinine, cystatin C (pCysC), and NGAL (pNGAL)), as well as eGFR based on these biomarkers. GFR was estimated using creatinine alone according to the bedside Schwartz equation (Schwartz bed ) [11], pCysC alone based on the Hoek equation [32], and both creatinine and pCysC based on the full age spectrum equation [33]. These cystatin C equations were chosen based on our previous evaluation of their ability to inform vancomycin CL during parametric population PK modeling [20].

Population PK Model Training
Next, we evaluated additional covariates on CL and V1 using a forward selection approach. All covariates were selected based on physiologic plausibility. Binary covariates included sex, receipt of vasopressor medications, and presence of augmented renal clearance (defined as eGFR > 130 mL/min/1.73 m 2 ; tested on CL only) [34], while continuous covariates included age, Pediatric Index of Mortality 3 (PIM3) score [35], and the above listed plasma and urinary biomarkers (on CL only). Urinary biomarkers were normalized to urine creatinine (i.e., [biomarker]/[uCr]) to account for urine volume. Continuous covariates were normalized to the median population value with the covariate effect evaluated using a power function; age was also evaluated as a Hill function on CL [36]. Meanwhile, binary covariates were parameterized using a linear proportional approach such that the covariate effect reflects a proportion increase/decrease of the typical parameter value in the absence of the covariate.
Covariate selection was guided by the principle of parsimony and by measures of goodness-of-fit. Models were evaluated at each step by inspection of observed-versuspredicted concentration plots, as well as examination of the model's bias, imprecision, regression coefficient, log-likelihood ratio (−2 × LL), and Akaike's information criterion (AIC) value. With the addition of a covariate to a model, a difference in the log-likelihood ratio of >3.84 was considered significantly improved fit (corresponding to p < 0.05 for one additional degree of freedom). When comparing models with the same degrees of freedom, lower AIC, bias, and imprecision of the observed-versus-predicted concentrations guided model selection.

Model Testing via Area under the Curve Comparisons
After model training, we evaluated the ability of the full model to predict vancomycin AUC for each subject in the model testing group. We first calculated the AUC for the dosing interval during which PK sampling was performed using his/her observed concentrations and linear up-log down trapezoidal noncompartmental analysis (NCA) [37]. Because of the timing of PK sampling, the drug was assumed to be at steady state, and therefore the same minimum concentration (C min ) was used prior to and after a given dose in these calculations. The 24 h "true" or observed AUC (AUC obs ) was then computed based on the subject's anticipated number of doses in a 24 h period.
We then used Bayesian estimation to evaluate how well our full model predicted AUC obs . The population joint density of the full model was employed as a Bayesian prior for each subject in the model testing group. Simulated concentration-time profiles were generated with predicted outputs each minute. The predicted concentration at the time of each measured PK sample, as well as at the start and end of his/her infusion, was recorded. The trapezoidal method was then taken, using all of the subject's predicted concentrations to calculate a predicted AUC 24 (AUC full ).
Recognizing that urinary biomarkers would not be routinely available in clinical practice, we also performed Bayesian estimation using two simpler models: a model incorporating only CysC-based eGFR (from Hoek equation) on CL and a model incorporating only creatinine-based eGFR (from bedside Schwartz equation) on CL. These models were chosen as comparators because they are more parsimonious and contain covariates that are readily available at all (creatinine) or many (CysC) pediatric institutions, which we felt would allow for easier clinical implementation. As above, AUC 24 was calculated for both the Hoek (AUC Hoek ) and Schwartz (AUC Schwartz ) models.
To assess the predictive performance of our models for AUC estimation, we determined the bias and imprecision of the AUC predictions from each model for each subject. Bias was calculated as (AUC pred − AUC obs )/AUC obs × 100 and imprecision as |AUC pred − AUC obs |/AUC obs × 100, where AUC pred is the generic AUC 24 predicted by our models: AUC full , AUC Hoek , and AUC Schwartz . Data were then summarized as the median bias (median percentage predictive error) and imprecision (median absolute percentage predictive error) for each model, as well as the fraction of subjects whose AUC pred was within 5%, 10%, 15%, and 20% of his/her AUC obs . Spearman correlation between AUC pred and AUC obs was also determined for each model. The a priori acceptance criteria were that AUC pred was within 20% of AUC obs in 85% of subjects and that the correlation between AUC pred and AUC obs was >0.9.

Area under the Curve Estimation from Optimal Sampling
The goal of Bayesian AUC estimation in clinical practice is to accurately estimate vancomycin AUC 24 using limited PK sampling. As such, we also sought to evaluate the ability of our model(s) to estimate AUC 24 using one and two optimally timed PK samples per subject. To achieve this, we utilized the multiple-model optimization algorithm (MMopt) in Pmetrics, which finds the optimal times based on when all the PK curves generated by the support points in the nonparametric model are most separated (e.g., the time points that are most informative), minimizing the Bayesian risk of misclassifying an individual as the wrong set of support points [38]. For each subject in the model testing group, MMopt was used identify the one and two most informative sampling time points for Bayesian estimation of AUC. We then created a reduced dataset for each subject that contained only the one or two observed concentrations closest to the optimal sampling time(s) identified by MMopt. We repeated the processes above in 2.d. with the population joint density of each model employed as a Bayesian prior and simulated concentrationtime profiles generated using only these reduced datasets. Predicted concentrations were again recorded and AUC 24 was calculated via the trapezoidal method. Bias, imprecision, and correlation were determined for each model, this time based on predictions from limited sampling.

Study Population
In total, 85 patients provided consent to participate. Five were deemed screen failures following consent and were excluded prior to performance of any study procedures. Meanwhile, 29 subjects were unevaluable because vancomycin was discontinued by the clinical team prior to PK sampling (n = 28) or could not undergo PK sampling (n = 1). One additional subject was excluded due to laboratory processing issues of the PK samples, making results uninterpretable (results reported out of order). As a result, 50 individuals participated in our study and were fully evaluable. Of these, 44 were eligible to be in the testing group (i.e., had ≥4 PK samples available), of whom we randomly chose 20. The characteristics of the model training and test groups are shown in Table 1. In general, the groups were similar, although the test group were less often on vasopressors at the start of vancomycin and received slightly larger vancomycin doses at the time of PK sampling.

Population PK Model Training
Thirty subjects contributed 150 vancomycin concentrations (14 clinical samples and 136 research PK samples) towards model training. The range of vancomycin concentrations was 3.9 to 67.8 µg/mL. Two concentrations were below the limit of quantification and coded as LLOQ/2 (i.e., 1.5 µg/mL). Observed concentrations vs. time after dose are plotted in Supplementary Figure S1. Table S2 displays the model training steps. A two-compartment model best described the data (AIC 820.5 vs. 879.6 for one-compartment). When evaluating renal function covariates on vancomycin CL, cystatin C outperformed SCr (AIC difference of −14.6) and pNGAL (AIC difference of −14.3). Meanwhile, eGFR based on the Hoek equation had a lower AIC (793.3) compared to eGFR based on the Schwartz equation (807.6) or the full age spectrum equation (801.5). Incorporation of eGFR Hoek or pCysC on CL resulted in comparable −2 × LL, AIC, and population predicted vs. observed R 2 values. We proceeded with further model training using the eGFR Hoek model since clinical dosing guidance is typically based on a patient's eGFR rather than a direct biomarker result.
When evaluating the addition of other covariates on CL and V1, numerous covariates, including urinary biomarkers on CL, improved model fitness based on a reduction in −2 × LL and AIC (Table S2). However, the incorporation of the uNGAL on CL, parameterized as the natural logarithm of uNGAL/uCr, led to the largest −2 × LL and AIC reductions. The correlation between log uNGAL and eGFR Hoek was low (−0.27). No additional covariates were informative on vancomycin CL or V1; thus, this model constituted the full model for testing. Observed versus population and individual predicted concentration plots of the full model are shown in the Supplemental Files ( Figure S2), and the population PK parameter estimates from this model are shown in Table 2.

Area under the Curve Comparisons
Twenty subjects comprised the model testing group. These subjects were similar to the model training group in terms of age, weight, PIM3 scores at PK sampling, timing of PK sampling, and vancomycin dosages received (Table 1), although they less often received vasopressors and had higher eGFR at vancomycin initiation. Biomarker concentrations at the time of PK sampling were also similar between the two groups (Table S3).
The median AUC 24 among the 20 subjects calculated using the trapezoidal method and observed concentrations (AUC obs ) was 456 mg h/L (range: 367-885). When using estimated concentrations from the individual Bayesian posteriors and our full model (AUC full ), the median AUC 24 was 475 mg h/L (range: 327-907). The median bias and imprecision of AUC full compared to AUC obs were 2.3% and 6.2%, respectively, and the correlation between the AUC estimates was 0.939. Meanwhile, 90% of subjects' AUC full estimates were with 20% of their AUC obs .
The PK parameter estimates for the two simplified models (based on eGFR Hoek and eGFR Schwartz ) are shown in Table S4. The median bias and imprecision of AUC Hoek compared to AUC obs were 1.8% and 7.0%, respectively, and the correlation between AUC Hoek and AUC obs was 0.926. The median bias and imprecision of AUC Schwartz compared to AUC obs were -2.4% and 6.2%, respectively, and the correlation between AUC Schwartz and AUC obs was 0.893. Additionally, 85% and 95% of subjects' AUC estimates were with 20% of their AUC obs using the Hoek and Schwartz models, respectively.
Observed versus individual predicted concentrations (i.e., Bayesian posteriors) for each of the three models are shown in Figure 2, while Table 3 displays the performance parameters of these models to estimate AUC 24 . The optimal sampling times relative to the end of the infusion for each of the three models are shown in Table S5. When fitting all available PK samples, the three models performed similarly (with AUC obs as the reference). However, when fitting the one or two PK samples closest to the MMopt optimal times, the performance of our full model declined. The imprecision of AUC full more than doubled when using one or two PK samples, and fewer subjects' AUC 24 were within 20% of AUC obs . However, ≥85% of subjects' AUC 24 from the Hoek and Schwartz models were within 20% of AUC obs when using limited sampling.

Discussion
In this population PK study, cystatin C was superior to the traditional renal function biomarker of serum creatinine as a marker of vancomycin clearance in critically ill children. This is consistent with other studies that described better correlation between vancomycin CL and CysC-based eGFR than creatinine-based eGFR in noncritically ill and critically ill pediatric patients [20,39]. Contrary to our previous work [20], however, we found that eGFR based on CysC alone (using the Hoek equation) outperformed eGFR based on both CysC and creatinine (full age spectrum equation). This highlights the inadequacy of creatinine as a renal function marker in critically ill children and suggests that routine measurement of CysC when administering vancomycin may be a more reliable approach. As the availability of cystatin C at pediatric institutions increases, and experience with cystatin C-based eGFR equations mounts, cystatin C may replace creatinine as the biomarker to inform dosing of other renally eliminated drugs in pediatrics as well. With increased appreciation of the negative ramifications of AKI on clinical outcomes, it is crucial to provide safe doses of nephrotoxic medications such as vancomycin. Given the well-recognized limitations of creatinine in pediatric patients, particularly in the ICU setting, it may be time to move towards a better biomarker in our most vulnerable children.
A goal of this work was to explore the potential value of novel urinary biomarkers to describe vancomycin disposition in critically ill children using a popPK modeling approach. Since changes in serum creatinine and other blood markers of kidney function may be delayed in the setting of AKI [40], we hypothesized that urinary biomarkers could facilitate detection of fluctuations in vancomycin CL before blood biomarkers. In fact, urinary NGAL was an influential covariate on vancomycin CL during our full model training; other urinary biomarkers (KIM-1, CysC, osteopontin) also led to large reductions in the AIC during model training steps, although not to the same extent as NGAL. We believe that these findings are important and warrant further investigation. Rather than solely relying on blood biomarkers to signify kidney function, bedside measurement of urinary biomarkers could provide insight into which patients have subclinical kidney injury (i.e., not detected via blood biomarkers) and may require early, pre-emptive dose adjustments. Although biomarkers did not improve estimation of AUC 24 estimation in our testing group, they could potentially be used as screening tests to identify patients at highest risk for impaired kidney function and/or toxicity.
Urinary biomarkers can detect subclinical kidney injury, but the majority of patients in our study actually had augmented renal clearance (ARC). This phenomenon describes a state of hyperfiltration and increased drug clearance, which can be detrimental in critically ill patients with severe infections, and it is unclear how reliable urinary biomarkers are in that clinical situation. Although the precise definition of ARC in children has been debated, an eGFR >130 mL/min/1.73 m 2 has been utilized in other vancomycin studies [41]. In our study, half of each of the model training and testing groups had ARC by Hoek equation calculation. This may explain why our full model failed to outperform the Schwartz and Hoek model in terms of AUC 24 estimation, as urinary biomarkers may not perform as well in patients with ARC as they do in those with AKI. Similarly, eGFR equations are generally more accurate in patients with impaired renal function rather than in ARC. Given the small sample size of our model testing group and limited number of subjects with impaired renal function, we may have been unable to fully demonstrate the value of novel biomarkers (both urinary and plasma) when it comes to AUC estimation.
It is also possible that urinary NGAL measurement is not precise enough to substantially influence AUC 24 estimation, particularly when using limited PK sampling. Urinary NGAL ranged from 1.4 to 2809 ng/mg creatinine in the model training group and from 7.8 to 10,034 ng/mg creatinine in the model testing group. Although measurements were not significantly different between the groups using Wilcoxon rank sum testing, this large variability in urinary biomarker values may preclude precise AUC estimation at an individual level using Bayesian methods. While we did not explore utilizing cut-points to categorize biomarker values (e.g., low, medium, high values), that is a potential avenue for future investigation.
Another goal of this study was to develop a popPK model that could be used to inform AUC-based vancomycin dosing using Bayesian estimation and limited sampling. We utilized a nonparametric popPK modeling approach, which differs from the more widely used parametric popPK methods in that it assumes that the distribution of parameter values is not necessarily described by an a priori defined continuous function (e.g., normal distribution) [42]. As a result, more commonly used statistical measures of variability in parametric popPK models (e.g., mean, standard deviation, coefficient of variation) may not fully describe the shape and structure of a nonparametric distribution and, as such, the idea of typical parameter values and interindividual variability around them does not apply. However, because nonparametric approaches allow for parameter probability distributions to take any shape, subpopulations and outliers can be detected, which is ideal in a popPK study of a highly dynamic population such as critically ill children, and methods for Bayesian estimation of individual PK parameters are robust. A more thorough review of nonparametric and parametric popPK approaches has been published [43].
At the time that this study started, AUC estimation was not standard of care. However, AUC-based dosing is now routine at our institution for anyone receiving >48 h of treatment. Clinical pharmacists implement log-linear methods to calculate AUC based on two vancomycin concentrations. To provide a distinct advantage over this approach, we felt it was important to specifically evaluate the utility of our models to inform AUC estimation based on a single, optimally timed vancomycin measurement. Both the Hoek and Schwartz models accurately estimated AUC 24 (within 20% of AUC NCA ) in 85% of subjects from a single sample with a median bias of 2.5% and −0.9%, respectively. An ongoing, prospective observational study (NCT05691309) at our institution is evaluating how well Bayesian AUC using our Hoek model aligns with clinical AUC calculations (using log-linear calculations) and assessing the ways in which this would influence dosing recommendations.
There are limitations to our study. First, we did not enroll children younger than one year of age and so our model(s) cannot be applied to infants. Cystatin C values are affected by age in this group, likely due to renal maturation occurring over the first year of life; thus, infants were specifically excluded from our study. A recent popPK study in critically ill neonates found that a model including creatinine on vancomycin CL performed similarly to a model including cystatin C and age [43]; thus, cystatin C may not be a superior renal function biomarker in the neonatal population. Second, because we required that a subject have at least four PK samples to be included in the model testing group, it is possible that our process of group assignment introduced bias; hence, the groups were different. We did not detect any statistically significant differences in the blood or urinary biomarkers between these groups, which were the covariates included in our models, but it is possible that other factors not considered covariates during model development in the training group were influential on vancomycin PK and AUC 24 estimation within the smaller model testing group. Third, recently published data suggest that urine chemistry analytes can be affected by collection of urine using cotton balls [44]. We are unaware of studies that have evaluated this for the urinary biomarkers included in our study, but this should be explored in future studies as urine collection methods could be a potential source of variability that affect the association between biomarkers and vancomycin clearance. Lastly, urinary biomarkers are sensitive for detection of AKI and, thus, may fluctuate over time. Although we collected ≥2 samples per subject, it is possible that changes in biomarker concentrations occurred outside of the windows of our urine sample collection. Although understanding the fluctuations of urinary biomarkers in critically ill children would be of interest, the need to collect urine at very specific times or to collect samples more often would further limit the clinical applicability of these biomarkers to inform drug dosing.

Conclusions
The present study developed a two-compartment model to describe vancomycin PK in critically ill children. Plasma cystatin C and urinary NGAL were informative on vancomycin clearance in our full popPK model. However, the full popPK model was not superior for estimation of Bayesian AUC 24 compared to simpler popPK models that included only cystatin C-or creatinine-based eGFR on clearance. Future studies will evaluate the utility of our models to inform vancomycin dosing using Bayesian estimation compared to two-point log-linear regression calculations.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/pharmaceutics15051336/s1, Table S1: Goal PK sampling times relative to end of the vancomycin infusion; Table S2: Model training steps; Table S3: Biomarker concentrations at time of PK sampling; Table S4: Population PK parameter estimates for the full model, eGFR Hoek model, and eGFR Schwartz ; Table S5: Summary of optimal sampling times relative to the end of the vancomycin infusion among model testing group; Figure   Informed Consent Statement: Informed consent was obtained from all subjects' parents/legal guardians involved in the study.
Data Availability Statement: Data can be made available upon request.