Population Pharmacokinetics of Vancomycin in Chinese ICU Neonates: Initial Dosage Recommendations

The main goal of our study was to characterize the population pharmacokinetics of vancomycin in critically ill Chinese neonates to develop a pharmacokinetic model and investigate factors that have significant influences on the pharmacokinetics of vancomycin in this population. The study population consisted of 80 neonates in the neonatal intensive care unit (ICU) from which 165 trough and peak concentrations of vancomycin were obtained. Nonlinear mixed effect modeling was used to develop a population pharmacokinetic model for vancomycin. The stability and predictive ability of the final model were evaluated based on diagnostic plots, normalized prediction distribution errors and the bootstrap method. Serum creatinine (Scr) and body weight were significant covariates on the clearance of vancomycin. The average clearance was 0.309 L/h for a neonate with Scr of 23.3 μmol/L and body weight of 2.9 kg. No obvious ethnic differences in the clearance of vancomycin were found relative to the earlier studies of Caucasian neonates. Moreover, the established model indicated that in patients with a greater renal clearance status, especially Scr < 15 μmol/L, current guideline recommendations would likely not achieve therapeutic area under the concentration-time curve over 24 h/minimum inhibitory concentration (AUC24h/MIC) ≥ 400. The exceptions to this are British National Formulary (2016–2017), Blue Book (2016) and Neofax (2017). Recommended dose regimens for neonates with different Scr levels and postmenstrual ages were estimated based on Monte Carlo simulations and the established model. These findings will be valuable for developing individualized dosage regimens in the neonatal ICU setting.


INTRODUCTION
After more than 60 years of widespread clinical use, vancomycin remains the gold standard antibiotic prescribed for the treatment of sepsis caused by coagulase-negative Staphylococci and methicillin-resistant Staphylococcus aureus in neonatal intensive care units (ICU) (Tong et al., 2016). In 2011, guidelines issued by the Infectious Disease Society of America and pediatricspecific guidance recommended targeting trough concentrations of >15 mg/L for critically ill children and >10 mg/L for all other pediatric patients (Liu et al., 2011). Low concentrations may result in less-effective therapy and an increased propensity for the development of bacterial resistance, whereas high concentrations are reported to be associated with nephrotoxicity and ototoxicity, although these toxic effects are less common in neonates .
Vancomycin is mainly eliminated by the kidneys and is highly correlated with creatinine clearance. In neonates, the pharmacokinetics of vancomycin are highly variable because of developmental and pathophysiological changes (Stockmann et al., 2015).It is both challenging and imperative to optimize the vancomycin dosage regimen to achieve adequate exposure within a short period of time.
The maximum a posteriori Bayesian estimation method has already been used to support vancomycin dosing decisions in adults (Deng et al., 2013;Jacqz-Aigrain et al., 2015) and children (Le et al., 2014). To obtain accurate estimation with this method for individualized therapy, it is crucial that reliable population pharmacokinetic characteristics are known for the target patients. Several pharmacokinetic studies have been conducted in different ethnic groups of adults and showed different clearance (CL) for vancomycin both in Chinese patients [6.05 (2.38-6.06) L/h and 6.06 ± 2.46 L/h] (He et al., 2014;Lin et al., 2015) and in Caucasian patients (4.03 ± 1.7 and 5.83 ± 2.39 L/h) (Guay et al., 1993;Sánchez et al., 2010). Moreover, several population pharmacokinetic studies have been conducted in Caucasian (Seay et al., 1994;Grimsley and Thomson, 1999;Capparelli et al., 2001;Mulla and Pooboni, 2005;Allegaert et al., 2007;Anderson et al., 2007;Marques Minana et al., 2010;Oudin et al., 2011;Mehrotra et al., 2012;Zhao et al., 2013;Frymoyer et al., 2014), Japanese (Kimura, 2002), and Malaysian (Lo et al., 2010) neonates, but a large variability in clearance was observed among groups. The average clearance across all three populations ranged from 0.08 to 0.50 L/h for a neonatal patient with a weight of 3 kg and postmenstrual age (PMA) of 40 weeks.
As little information about vancomycin population pharmacokinetics is known for Chinese patients, the goal of the present study was to establish a population pharmacokinetic model for critically ill Chinese neonates receiving vancomycin therapy and to provide a rational dosage regimen for Chinese neonates.

Patients
Chinese neonates who received vancomycin in the neonatal ICU at Shanghai Children's Hospital (Shanghai, China) between January 2013 and December 2016 were enrolled. Patients included in this study were neonates with PMA≤ 48 weeks for preterm neonates and with postnatal age (PNA) ≤ 28 days for term neonates. All patients were treated with vancomycin for at least 3 days, and at least one vancomycin concentration was assayed. Patients with extracorporeal membrane oxygenation or continuous renal replacement therapy were excluded from the current study. The following information was retrospectively collected from electronic medical records: gestational age (GA), PMA, PNA, body weight (WT), dosing history and concentration of vancomycin, serum creatinine (Scr) levels, clinical laboratory tests of other renal and hepatic functions, and co-administered medications.
This study was carried out in accordance with the recommendations of the Declaration of Helsinki (2000). The protocol was approved by the Ethics Committee of Shanghai Children's Hospital. Parents or guardians of patients gave their informed consent before enrollment.
The dose of vancomycin (Vancocin, Lilly, S.A, Suzhou, China) was 10-15 mg/kg, which was administrated every 8 h (q8h) or every 12 h (q12h) with a 2-h infusion, in accordance with local protocols. Blood samples were collected 1 h after completion of drug infusion (peak concentration) or half an hour before the start of vancomycin administration (trough concentration) for each patient. Trough and peak concentrations were obtained after at least four repeated doses. Potential outlier data points (observations) were identified by employing conditional weighted residual (CWRES) results outside a range of ± 6 ( Byon et al., 2013).

Bioassay
Concentrations of vancomycin were determined using the fluorescence polarization immunoassay with an ARCHITECT i2000SR (Abbott Laboratories, Chicago, IL, USA). The limit of detection was 1 mg/L, and the calibration range of this assay was 3 to 50 mg/L. The intra-day and inter-day coefficients of variation (CVs) were all < 20%.
Scr assays were performed with the enzymatic method and were analyzed with a 7180 Automatic Analyzer (Hitachi High-Tech Science Systems Corporation, Tokyo, Japan). The calibration range was from 3 to 100 mg/L. The intra-day and inter-day CVs were all < 3.75%. Creatinine can pass through the placenta and many endogenous factors can influence creatinine determination in neonates, which can limit the concentration of creatinine detected.

Base Model
The vancomycin concentration data were fitted by the one or two compartment with a first-order elimination model. NONMEM subroutines were specified as ADVAN1-TRANS2 or ADVAN3-TRANS4, respectively.
An exponential model (Equation 1) was used to account for inter-individual variability (IIV): where TV(P) is the typical value of the pharmacokinetic parameter, and P i refers to the pharmacokinetic parameter of the ith patient with random variable η i , which is normally distributed with a mean of zero and variance of ω 2 . Residual unexplained variability was tested by an additive model (Equation 2), exponential model (Equation 3), or a combined additive and exponential model (Equation 4): where Y represents the observation, IPRED is the individual predicted concentration, and ε is a symmetrically distributed variable, with a mean of zero and variance of σ 2 .

Covariate Models
The continuous covariates, including WT, GA, PNA, PMA, Scr, blood urea nitrogen, serum albumin concentration, aspartate transaminase and alanine transaminase, and categorical covariates, including gender, concomitant drug (ceftriaxone, meropenem, gentamicin, furosemide, ibuprofen, and dexamethasone), were screened for their influence on clearance and the volume of distribution. Covariate screening was performed according to the following four steps: Step 1 Body weight and age have significant impacts on the pharmacokinetics of vancomycin in neonates (Wallis and Williamson, 2011;Jacqz-Aigrain et al., 2015) and physical maturation is a time-dependent characteristic that must be considered in neonates (Schmidt and Derendorf, 2014). Therefore, body weight (WT) and age (PNA, GA, and PMA) were screened first. Four different models based on allometric scaling were tested using Equation (5): where COV median is the median of the covariate, MF is the maturation factor that is defined as the process of becoming mature. The model displaying the best fit was selected for further analysis.
Model I: In the simplest exponent model, the exponent θ was estimated, and MF was fixed to 1, indicating that maturation was not considered. This model is shown as Equation (6): Model II: For the maturation model (Holford et al., 2013), the exponent θ was assigned a fixed value of 0.75, and MF was calculated according to Equation (7): where TM 50 is the age (in terms of GA, PMA, or PNA) at which clearance maturation reaches 50% of that of adults, and Hill is the slope parameter for the sigmoid E max maturation model. Model III: This is referred to as the WT-dependent exponent model (Holford et al., 2013): Model IV: This is referred to as the age-dependent exponent model (Ding et al., 2015): where θ 0 is the value of the exponent at a theoretical WT of zero (Equation 8) or at birth (0 years) (Equation 9), k max is the maximum decrease of the exponent, k 50 is the WT (Equation 8) or age (Equation 9) at which a 50% decrease relative to the maximum decrease is attained, and the Hill coefficient is used to determine the steepness of the sigmoid decline.
Step 2 In previous studies, renal function has been identified as an important covariate (Grimsley and Thomson, 1999;Capparelli et al., 2001;Kimura, 2002;Oudin et al., 2011;Mehrotra et al., 2012;Zhao et al., 2013;Derschmills et al., 2014;Frymoyer et al., 2014). We thus investigated and Scr by using the exponential model (Equation 6) and a linear model (Equation 10): The one displaying the best fit was used for further analysis.
Step 3. The remaining covariates were then accessed sequentially by forward inclusion and backward elimination approaches using the exponential model (Equation 8) or linear model (Equation 10) for continuous variables and a proportional model (Equation 11) for categorical variables, such as gender in the P i of vancomycin: Where θ is the fractional change in TV(P) for males.
Step 4 Taking into account the rapid variation in the physical status of neonates, the covariates identified above could be defined as time-varying covariates to illustrate IIV. The timevarying covariates model splits the individual covariate effects into baseline and change-from-baseline effects. Two different models based on covariate scaling of the pharmacokinetic parameters were tested using Equations (12) and (13): where BCOV is the baseline value of the covariate; BCOV median is the median of the baseline value of covariate; θ BCOV describes the effect of IIV. DCOV is equal to the current covariate value minus BCOV at each time point and corresponds to the fractional change in the typical value with each unit difference in BCOV relative to BCOV median. θ DCOV describes the effect of covariate variation within an individual and is the fractional change in the typical value with individual changes in COV, and (η , P i ) refers to variable with a mean of zero and variance of ω 2 that describes the random effect of P i . If θ BCOV and θ DCOV are different, DCOV is fixed to zero, and an additional variance parameter (η COV, P i ) that accounts for IIV to influence the covariates for the population parameter estimates is included (Equation 14).
θ COV is the parameter estimation value of the covariate, COV is the value of the covariate, COV median is the median of the covariate, η COV is random variable (with zero mean and variance ω 2 ).
Nested models in covariate screening were compared statistically using a likelihood ratio test on the differences in the objective function value (OFV). The change was considered significant if the decrease in OFV was >3.83 (χ 2 , df = 1, P < 0.05) for the forward inclusion step, and the increase in OFV was >6.63 (χ 2 , df = 1, P < 0.01) for the backward elimination step.
The covariates in the model were selected based on physiological plausibility of parameter estimates, goodness-of-fit plots, and statistical significance.

Model Evaluation
The performance of the final model was first evaluated by visual inspection of diagnostic goodness-of-fit plots. Goodness-of-fit plots included the following scatterplots: observation (DV) vs. individual prediction (IPRED), DV vs. population prediction (PRED), conditional weighted residual errors (CWRES) vs. time, and CWRES vs. PRED (Hooker et al., 2007).
The robustness of the model was assessed using a nonparametric bootstrap (Ette, 1997), with repetition of 2000 NONMEM runs of the final model. The bootstrap median parameter values and the percentile bootstrap 95% confidence intervals were compared with the respective values estimated from the final model.
Normalized prediction distribution error (NPDE) (Comets et al., 2008) was used to evaluate the predictive performance of the model based on a Monte Carlo simulation with the R package (version 2.0, http://www.npde.biostat.fr/). NPDE results were summarized graphically using (1) quantile-quantile plot of the NPDE, (2) a histogram of the NPDE, (3) scatterplot of NPDE vs. time, and (4) scatterplot of NPDE vs. PRED. If the predictive performance is satisfied, the NPDE will follow a normal distribution (Shapiro-Wilk test) with a mean value of zero (t-test) and a variance of one (Fisher's test).

Dosing Regimen Design
The final established population pharmacokinetic model was used to obtain dosing regimens for vancomycin to reach AUC 24h /MIC ≥ 400 which is known to produce an effective therapeutic outcome (Jacqz-Aigrain et al., 2015).When MIC = 1 mg/L, the daily dose can be calculated by the final model for determining the clearance and Equation 14: Simulations were performed for virtual patients with various levels of renal function and ages, to determine the most appropriate scheme to satisfy the therapeutic criteria. For this purpose, 1,000 replicates of each scenario were simulated by the Monte Carlo method and were completed by the $SIMULATION modules in NONMEM software. Virtual patients were designated as having PMAs covering a 2week window between 28 and 44 weeks, and their corresponding WTs were calculated according to the World Health Organization growth chart for infants (Centers for Disease Control Prevention, 2009). Patients were designated a Scr level of 15, 20, or 35 µmol/L and with an age of PNA 7 days (≤ 1 week) or PNA 15 days (>2 weeks).

Patients
Data from 165 vancomycin measurements, with a trough concentration for 75 of the measurements and a peak concentration for 90 of the measurements, from 80 subjects were    Frontiers in Pharmacology | www.frontiersin.org available to perform the population pharmacokinetics analysis. For each subject, an average of two samples were obtained .The GA range was 25.7 weeks to 41.1 weeks, with a mean WT of 2.87 kg. Of these patients, 59% were preterm infants and 57% had respiratory tract infections. All the observation were included, and no outlier records were identified. Clinical characteristics of the neonates included in the analysis are summarized in Table 1.

Model Building
A one-compartment model with first order elimination described the pharmacokinetic of vancomycin. The residual unexplained variability was described best by a proportional model. As only peak and trough samples were collected, the relative standard error (RSE) of the IIV of the volume of distribution was poor (>50%) and was not estimated.
For the first step of covariate screening, several models were tested, and the results are shown in Table 2. Among the four models examined, the simple exponent model (Model I) and maturation model (Model II) had lower Akaike information criteria and Bayesian information criteria than the ADE model (Model III) and BDE model (Model IV). The maturation model had a condition value of 23,082, much greater than 1,000, indicating model instability (Byon et al., 2013). Moreover, the RSE of most PK parameters in the maturation model were more than 100%, implying inaccuracy of the model parameter estimates. Therefore, the simple exponent model was employed in further analyses.
Second, Scr and WT were identified as significant covariates and were thus retained in the model. Further incorporation of time-varying covariates did not improve the model performance,   (16) where CL is vancomycin clearance, WT is current body weight in neonates, Scr is serum creatinine, V is volume of distribution for vancomycin. The final model parameter estimates are shown in Table 3.

Model Evaluation
Goodness-of-fit plots for the base model and final model are shown in Figure 1. Compared with the base model, the final model showed no obvious bias or significant trends within these scatterplots. Moreover, the data fitting for the final model was much improved relative to that of the base model. The results from the bootstrap procedure are shown in Table 3. The median values from the bootstrap procedure were close to the parameter estimates from the NONMEM, with < 5% bias. In addition, more than 99% of the bootstrap runs were successful, indicating that the model was stable.
The NPDE distribution and histogram are presented in Figure 2. The assumption of a normal distribution for the differences between predictions and observations was acceptable. The quantile-quantile plots and histogram also confirmed the normality of the NPDE (Figure 2).

Dosing Regimen Design
The dosage regimens recommended by current guidelines and estimated by the established model are displayed in Figure 3 and Table 4. The guideline schemes from the BNF for children, the Blue Book and the Neonatal Formulary are nearly in accordance with the 15% to 85% dosage interval from the present model that achieves a therapeutic target of 400≤ AUC 24h /MIC < 800. The targeted vancomycin concentration based on the FDA labeled dosage, the Red Book, the Pediatric & Neonatal Dosage Handbook, and Neofax was not likely sufficient, especially for neonates with Scr of 15 µmol/L. Moreover, Neonatal Formulary shows an overdose for a subpopulation of neonates with Scr of 35 µmol/L, which might indicate an increased risk of toxicity.

DISCUSSION
In the present study, we first developed a population pharmacokinetic model for vancomycin in Chinese neonates in the ICU and found that WT and Scr levels have significant influences on clearance. Little obvious ethnic difference of vancomycin clearance was shown in Asian and Caucasian neonates from our study.
Differences in vancomycin pharmacokinetics have been noted between Asian and Caucasian populations based on two recently published studies (Lin et al., , 2016. To expand this analysis, we first looked at a group of 12 previous studies that examined neonatal populations from various countries. The population within some of these studies were heterogeneous. Patients with extracorporeal membrane oxygenation were included in the model of Mulla et al. (Mulla and Pooboni, 2005), and the majority of neonates used to develop the models of Allegaert et al. (2007) and Capparelli et al. (2001) were preterm. Due to differences in the physiologic development of neonates within these populations, models from these three studies were excluded from the comparison. Additionally, neonates within 6 other studies (Seay et al., 1994;Grimsley and Thomson, 1999;Lo et al., 2010;Marques Minana et al., 2010;Oudin et al., 2011;Zhao et al., 2013) were smaller, with a mean WT <2 kg. Given the relative weight of neonates in our study and the impact of WT on vancomycin pharmacokinetic parameters, these six studies were excluded from comparisons. The three remaining studies (Kimura, 2002;Mehrotra et al., 2012;Frymoyer et al., 2014) were included in our analysis of ethnic differences regarding vancomycin clearance.
The clearance of a standardized patient, as determined by the different models, was calculated to investigate differences in vancomycin clearance relative to ethnicity. The standardized patient had a WT of 2.8 kg, PMA of 37 weeks with different Scr levels.
As shown in Table 5, across varying Scr levels ranging from 20 µmol/L to 50 µmol/L. there was 27∼39 % difference in vancomycin clearance between the current study and the study by Mehrotra (Mehrotra et al., 2012), but comparable to another study (<13%) which was also conducted in US (Frymoyer et al., 2014). This information does not support the conclusion that there are obvious differences in vancomycin clearance between Chinese and Caucasians. However, the estimated clearance was much higher than the study conducted in 19 Japanese neonates (Kimura, 2002). The reason was unclear. The present study might be under-power to conclure ethnic impact on vancomycin PK in neonates. Factors in the current study, including analytical methods, and disease progression could affect the assessment of ethnic differences. Dosage recommendations by the label and reference books are variable as shown in Table 6. The variabilities were attributed to the different covariates considered within these recommendations. Body weight is the most notable covariate for vancomycin dosing found within all references. Age-based (as PMA, PNA, and GA) dosing is also included in most of the references, such as the FDA label recommendation, BNF for children, Blue Book, Neonatal Formulary, Pediatric and Neonatal Dosage Handbook and Neofax. Dosing based on Scr is only included in 2 references, the Red Book, and Pediatric & Neonatal Dosage Handbook and only covered Scr level > 61.9 µmol/L. However, Scr was identified to have large impacts on the CL of vancomycin in all previous population pharmacokinetic studies as well as the present study.
Based off data from our final model, we found that the current recommended doses of vancomycin from FDA labeled dosage, Red Book, Pediatric & Neonatal Dosage Handbook, and Neofax may be inadequate to meet a treatment target of AUC 24h /MIC ≥ 400, especially for patients with a greater renal clearance status, especially Scr < 15 µmol/L. This finding is consistent with several previous studies (Krivoy et al., 1998;Liu et al., 2011;Abdel et al., 2015;Zhao et al., 2015). They found that the usual recommended dose of 60 mg/kg/day  Frontiers in Pharmacology | www.frontiersin.org did not achieve vancomycin pharmacodynamic targets in most patients. Silva, D C et al (Silva et al., 2012) reported that a vancomycin dose of 81 mg/kg/day was required to achieve an AUC 24h /MIC > 400 in 56% of patients. Doses as high as 120 mg/kg/day were also recommended to improve the therapeutic pharmacodynamic targets (Abdel et al., 2015). Higher than usual vancomycin doses may be required to treat patients with severe Gram-positive infections.
There are several limitations to this study. The current guidelines or consensus for therapeutic drug management (TDM) of vancomycin in the United States (Rybak et al., 2009), Japan (Matsumoto et al., 2013) and China (Ye et al., 2016) recommends that trough concentrations should be collected at regular intervals in clinical settings. Therefore, only peak and trough concentrations were collected in our study and a one compartment model was applied for the structural model even though vancomycin is more normally modeled with a twocompartment model. The simplification to a one compartment model may lead to deviation of clearance estimation. However, the bias usually is <20% and does not obviously affect the estimation of the area under the curve (AUC). (Ling et al., 2014) (Kowalski and Hutmacher, 2001). The recommended regimen was based on a study population from our hospital. Therefore, generalization to other ICU neonates treated with vancomycin, especially if their covariate characteristics lie outside the range of our study population, would require additional investigations.
Moreover, several cofactors may affect the clinical outcomes of patients, such as disease progression, baseline weight, gestational age, and medication interactions. Therefore, with a limited number of patients enrolled in the study, we did not compare the clinical outcomes between those who had different exposure levels. Furthermore, as a pharmacodynamics indicator, the MIC value was obtained for only a few patients, such that we could not build a population pharmacokinetic-pharmacodynamics model, which would have been a better predictor of vancomycin's therapeutic effect.
This study offers initial vancomycin dosing regimen with varying degrees of PMA, WT, and serum creatinine for neonates. For patients with complex disease conditions, Bayesian approaches might be used to provide individualized dose recommendations instead of look-up tables or nomograms. Dose calculators and other decision support tools based on population pharmacokinetic (PPK) models (Fuchs et al., 2013) could contribute to simplifying the complex Bayes calculations and making them more intuitive to the user in clinical practice.
Recently, a nonparametric (NP) population modeling approach was reported to have advantages in patient's individual dosing adjustment, which permits development of dosage regimens to hit desired therapeutic targets with maximum precision (Jelliffe et al., 2011;Neely et al., 2018). This will be further investigated.

CONCLUSION
In summary, this study has built a population pharmacokinetic model of vancomycin for Chinese neonates. WT and Scr levels were the important covariates, which affect clearance. Moreover, this study found no obvious differences in the clearance of vancomycin comparing Caucasian and Chinese neonates. For patients with normal renal function, the dosing recommendations are likely not sufficient based on the target of AUC 24h /MIC ≥400, with the exception of British National Formulary (2016-2017), Blue Book (2016) and Neofax (2017). However, these sources provide little information on dosing adjustments based on patient renal function, thus our model provides a method for adjusting the vancomycin dose accordingly.