This multicenter observational study was conducted in six major hospitals in Saudi Arabia. Data were retrospectively collected from the electronic medical records of the following hospitals: King Saud University Medical City in Riyadh, King Faisal Specialist Hospital in Riyadh, National Guard Hospital in Riyadh, King Faisal Specialist Hospital in Jeddah, Prince Sultan Military Medical City in Riyadh, and Armed Forces Hospital Southern Region in Abha. Data were retrieved for VLBW neonates (0–30 days old) who were admitted to the NICU and received vancomycin for proven or suspected MRSA or MRCONS infections. Patients were included if they received vancomycin for more than 48 h and had at least one vancomycin serum trough concentration measurement in a steady state. The collected data included PNA, PMA, weight, Scr, vancomycin dose, frequency and concentrations. For vancomycin initial dosing recommendations, online pediatric and neonatal Lexi-Drugs® (Lexicomp, Ohio, USA) (20) was the main neonatal dosing reference at National Guard Hospital, whereas all other participating hospitals relied mainly on Micromedex NeoFax (IBM Corp., Armonk, NY, United States) (21) (Table 1). Subsequent dose adjustments were guided by the vancomycin trough concentrations. In King Saud University Medical City, both vancomycin peak and trough concentrations were obtained for therapeutic drug monitoring. In all other hospitals, the practice was to collect trough concentrations only. Peak samples were collected 1 h after the end of infusions, and trough concentrations were collected 30 min before the next dose in a steady state. Exclusion criteria were: renal failure, receiving hemodialysis, missing dosing information and those sampled before the steady state were excluded. The study was IRB approved in all six participating centers.

Three different validated methods were used for the analytical assay, as each hospital had a different system. One hospital used a chemiluminescent microparticle immunoassay with an Architect i4000SR immunoassay analyzer (Abbott), two hospitals used the kinetic interaction of microparticles in solution (KIMS) with a Cobas c system (Roche/Hitachi), and three hospitals used the homogeneous particle-enhanced turbidimetric inhibition immunoassay (PETINIA) with an Alinity c analyzer (Abbott). The assay range for the Architect immunoassay is from 0.5 to 83 µg/ml, the range for the KIMS assay is 4–80 µg/ml, and the range for the PETINIA is 1.4–100 µg/ml. For the Scr assays, four centers used the Jaffe method, and two centers used the enzymatic method. The following formula was used to convert the enzymatic assay to the Jaffe method (22):Jaffeconcentration=0.122+Enzymaticconcentration1.05

We split the data randomly into the training (70%) and validation (30%) data sets. Data splitting was performed using R statistical software. Pharmacokinetic modeling was performed using the Monolix software (2020R1) with the stochastic approximation expectation maximization algorithm.

First, we developed a structural model for vancomycin. This included testing one- and two-compartment models with linear and nonlinear elimination. For the between-subject variability in pharmacokinetic parameters, we used the following lognormal distribution:pi=p.exp⁡(ηi),where pi is the individual pharmacokinetic estimate for the i^th individual, p is the median value of the pharmacokinetic parameter, and ηi is the standard deviation for the i^th individual with a mean of 0 and variance ω². For residual unexplained variability, we tested different statistical models, including the constant, proportional, and combined error models. As different assays were used for vancomycin, we tested having different residual variability models accordingly.

After we developed the base model, we evaluated the effect of covariates on pharmacokinetic parameters. We tested the effects of postnatal age, gestational age, PMA, weight, height, sex, Scr, presence of congenital heart disease, and VLBW vs. extremely low birth weight (ELBW) on all pharmacokinetic parameters. Covariates were tested by plotting the post hoc individual pharmacokinetic parameters (empirical Bayesian estimates) against covariates to identify possible correlations. Testing was performed stepwise using a likelihood ratio rest. The significance level was set at p < 0.05. The effect of continuous covariates was modeled using a power function as follows:Pi=P∗(XiX)β,where Pi is the individual pharmacokinetic parameter for the i^th individual, P is the typical value of the pharmacokinetic parameter in the population, Xi is the covariate value for the i^th individual, and Xi is the population mean of the covariate. The exponents for the effects of weight on total body clearance (Cl) and the volume of distribution (V) were fixed at 0.75 and 1, respectively (23–25). To determine the effect of age on Cl, we tested both power, linear, and sigmoidal functions to account for maturation (25, 26). The sigmoidal function is expressed as follows:FPMA=PMAhillPMAhill+TMA50hill,where hill is the Hill coefficient for clearance, and TMA₅₀ is the PMA, where Cl is 50% of the mature value.

Model evaluation and selection were guided by the objective function value, physiological plausibility, and standard diagnostic plots. The relative standard error of the pharmacokinetic estimates was also assessed, with an optimal value of <30%.

The final population pharmacokinetic model was evaluated using the validation dataset. Model performance was assessed using simulation diagnostics, Visual Predictive Check (VPC), and Normalized Prediction Distribution Error (NPDE) (27). The null hypothesis for the NPDE is that it follows a normal distribution, with a mean of 0 and a standard deviation of 1. The VPC is based on 1,000 dataset simulations, followed by a comparison between the observed data and the simulated concentrations.

The final model was used to perform Monte Carlo simulations to determine the optimal vancomycin dose. We evaluated various doses ranging from to 10–20 mg/kg in 2.5 mg/kg increments every 6, 8, 12, 18, and 24 h. The final population pharmacokinetic model with the distribution of the pharmacokinetic parameters was embedded into the “Simulx” package in R. Simulx is part of the “mlxr” package in R and can be used to simulate the time concentration profile for different pharmacokinetic models (Lixoft, Antony, France, http://simulx.webpopix.org/). To ensure that the distribution of the significant covariates was identical to our patient population, we replicated our dataset 40 times (40 × 236 = 9,440). We simulated the time-concentration profiles for all 9,440 cases. We evaluated two targets in the simulations: (1) AUC₀₋₂₄ of 400–600 µg.h/ml, as recommended in the IDSA treatment guidelines, (10) and (2) the less conservative target of AUC₀₋₂₄ of 400–800 µg.h/ml, based on a prior study in pediatrics demonstrating increased nephrotoxicity for AUC₀₋₂₄ > 800 µg.h/ml (28). We also estimated the probability of both trough levels >5 µg/ml and 20 µg/ml for each dosing regimen. Here, we considered troughs as the secondary marker for nephrotoxicity, as higher trough concentrations are also associated with nephrotoxicity (29–31). Prior studies in adults indicate risk of nephrotoxicity is approximately 15% for troughs between 15 and 20 µg/ml and increases to approximately 30% for troughs >20 µg/ml (30). Therefore, troughs >20 µg/ml should be always avoided, while troughs >15 µg/ml should be minimized if possible, as some patients can achieve a vancomycin therapeutic target of an AUC₀₋₂₄ 400–600 µg.h/ml at troughs <15 µg/ml (12, 30–32).

Vancomycin pharmacokinetic data were best described using a one-compartment model with linear elimination and a proportional error model. Significant covariates identified were weight for volume of distribution (V) and weight, Scr, and PMA for clearance (Cl). Bodyweight was scaled to the median bodyweight of 0.93 kg, and the exponents for the effect of weight on V and Cl were fixed at 1 and 0.75 (23, 25, 26). The effect of PMA was best described using a sigmoidal function. The typical Cl value for a VLBW neonate weighing 0.93 kg, PMA equal to 26 weeks, and Scr of 0.6 mg/dl was 0.09 L/h, while the typical V value for a VLBW neonate weighing 0.93 kg was 0.81 L. The formulas below describe the predicted values of Cl and V:V=0.81∗(Weight0.93)Cl=0.09∗(Weight0.93)0.75∗(0.6serumcreatinine)0.48∗PMA4.42PMA4.42+26.34.42Figure 1 shows the correlation between Cl normalized by body weight vs. Scr and PMA. The results of the final model are presented in Table 3. Diagnostic plots are shown in Figures 2–4.

The validation dataset contained 97 observations from 74 participants. The diagnostic plots are shown in Figures 5 (VPC) and Figure 6 (NPDE). For the VPC, the observed and predicted medians and percentiles were similar. The mean NPDE was 0.043, and the standard deviation was 1.1, indicating that the model had limited bias and captured most of the data variability.

Simulations were performed for various dosing regimens based on the final population pharmacokinetic model. The simulations were performed for different age and Scr groups as follows: PMA above and below 29 weeks and Scr < 0.6, from 0.6–0.9, and >0.9–1.2 mg/dl, for a total of six groups. There were only five patients with Scr > 1.2 mg/dl; therefore, the simulations were capped at 1.2 mg/dl. The range for Scr, PMA and weight in the simulation data set was: Scr from 0.2–1.2 mg/dl, PMA from 22 to 42 weeks and for bodyweight was from 0.46 to 2.2 kg. Table 4 shows the simulations for selected dosing regimens; we only presented the doses that meet these criteria: probability of AUC₀₋₂₄ > 400 µg/h/ml is >60%, and probability of AUC₀₋₂₄ > 800 µg.h/ml or trough >20 µg/ml is <10% to avoid nephrotoxicity. For each group, we identified 2–4 dosing regimens that meet these criteria. Dose selection would depend on the selected target: the narrower target of AUC₀₋₂₄ which is equal to 400–600 µg.h/ml, or the wider target of AUC₀₋₂₄ which is equal to 400–800 µg.h/ml.

Table 5 shows the recommended optimal dose depending on the selected target (400–600 or 400–800 µg.h/ml). The recommend doses to maintain an AUC₀₋₂₄ of 400–600 µg.h/ml are expected to achieve the therapeutic target in more than 60% of patients. If we select the wider range of 400–800 µg.h/ml, the recommended doses would achieve the target in more than 80% of patients. However, the 400–800 µg.h/ml target would require higher doses and may increase the probability of elevated trough concentrations and AUCs. The only exception was for neonates with PMA > 29 weeks and Scr < 0.6 mg/dl, the highest achieved probabilities were 56% for the target AUC₀₋₂₄ of 400–600 µg.h/ml and 74% for the target AUC₀₋₂₄ of 400–800 µg.h/ml with the selected dosing regimen.