External Evaluation of Vancomycin Population Pharmacokinetic Models at Two Clinical Centers

Background: Numerous vancomycin population pharmacokinetic models in neonates have been published; however, their predictive performances remain unknown. This study aims to evaluate their external predictability and explore the factors that might affect model performance. Methods: Published population pharmacokinetic models in neonates were identified from the literature and evaluated using datasets from two clinical centers, including 171 neonates with a total of 319 measurements of vancomycin levels. Predictive performance was assessed by prediction- and simulation-based diagnostics and Bayesian forecasting. Furthermore, the effect of model structure and a number of identified covariates was also investigated. Results: Eighteen published pharmacokinetic models of vancomycin were identified after a systematic literature search. Using prediction-based diagnostics, no model had a median prediction error of ≤ ± 15%, a median absolute prediction error of ≤30%, and a percentage of prediction error that fell within ±30% of >50%. A simulation-based visual predictive check of most models showed there were large deviations between observations and simulations. After Bayesian forecasting with one or two prior observations, the predicted performance improved significantly. Weight, age, and serum creatinine were identified as the most important covariates. Moreover, employing a maturation model based on weight and age as well as nonlinear model to incorporate serum creatinine level significantly improved predictive performance. Conclusion: The predictability of the pharmacokinetic models for vancomycin is closely related to the approach used for modeling covariates. Bayesian forecasting can significantly improve the predictive performance of models.


INTRODUCTION
Vancomycin is a glycopeptide antibiotic used as the gold standard treatment for serious infections in adults caused by Grampositive bacteria, especially methicillin-resistant Staphylococcus aureus (Pacifici and Allegaert, 2012). Vancomycin is also effective in infants with serious Gram-positive infections. However, the therapy window of vancomycin is narrow, and differences in neonatal development and pathophysiology result in high interindividual variability in vancomycin pharmacokinetics (Stockmann et al., 2015). Although excessive exposure to vancomycin can lead to side effects including ototoxicity and nephrotoxicity , under-dosing is often associated with treatment failure and patient mortality (Rybak et al., 2020). Therefore, despite the challenges, it is imperative to optimize vancomycin regimens in neonates.
Therapeutic drug monitoring is an applicable approach for the pharmaceutical care of vancomycin. According to the American Society of Health-System Pharmacists consensus (2020), the administration target for vancomycin is an area under the concentration-time curve (AUC)/minimum inhibitory concentration of ≥400 h in neonates and infants (Rybak et al., 2020). Although obtaining a sufficiently large number of samples to estimate the AUC is difficult in clinical practice, especially for neonates, a population pharmacokinetic analysis could provide sufficient pharmacokinetic parameters to estimate the AUC through sparse sampling. It is possible to model vancomycin dosing in neonates through reliable individual pharmacokinetic characteristics using Bayesian approaches.
Choosing the appropriate population PK model to estimate the initial and maintenance dosage for vancomycin is essential in clinical practice; however, the performance of most of the published pharmacokinetic models is still unknown. Zhao et al. (2013a) conducted an external evaluation of six models in neonates and found that the analytical method used for serum creatinine (SCR) is a crucial factor in explaining the variability of predictions among different studies. However, more than ten population PK studies have been conducted since then, using several new modeling strategies. Therefore, it is still worth evaluating all the published population pharmacokinetic models for vancomycin in neonates.
Our research aimed to systematically evaluate the published population pharmacokinetic models of vancomycin in neonates, using data from independent cohorts collected from two clinical centers. Moreover, factors that may influence model predictability were also investigated, such as structural model selection and covariate screening approaches, to provide informed guidance for future studies.

Review of Published popPK Studies
The PubMed, Scopus, and Web of Science databases were systematically searched for population pharmacokinetic analyses of vancomycin published up to October 2020. The key words "vancomycin," "pharmacokinetic" or "pharmacokinetics" or "model" or "nonlinear mixed effect model" were used in the search strategy. The publications were included if 1) the study was a population pharmacokinetic analysis of vancomycin in neonates and 2) the article was written in English. The publications were excluded if 1) the model was not created using a nonlinear mixed-effects modeling approach, or 2) the model could not be recreated using the published information, or 3) the modeling populations overlapped or the articles were duplicated.

External Evaluation Cohort
Patients Datasets were derived from published population PK studies conducted in neonates who received vancomycin at Shanghai Children's Hospital between January 2013 and December 2016 , and Suzhou Hospital Affiliated to Nanjing Medical University between September 2011 and March 2016 (Li et al., 2017). Patients included in these two studies were preterm neonates with a postmenstrual age (PMA) of ≤48 weeks and term neonates with a postnatal age (PNA) of ≤28 days. All patients were treated with vancomycin for at least 3 days, and at least one vancomycin level was determined based on routine therapeutic drug monitoring. Patients with extracorporeal membrane oxygenation or who were on continuous renal replacement therapy were excluded from this study.
The following information was collected in each study: gestational age (GA), PMA, PNA, current weight (WT), birth weight, dosing records, measurements of vancomycin levels, and SCR level.
The doses of vancomycin ranged from 10 to 15 mg/kg, administered every 8 h or every 12 h with a 1 h or 2 h infusion duration. Peak samples were collected 1 h after completion of drug infusion, and trough samples were collected half an hour before vancomycin administration in each neonate. Trough and peak levels were determined after at least four repeated doses.

Bioassay
Vancomycin levels were determined using a fluorescence polarization immunoassay with an Architect i2000SR (Abbott Laboratories, Chicago, IL, UNITED STATES). The limit of detection was 1 mg/L, and the calibration range was 3-50 mg/ L. The intra-day and inter-day coefficients of variation were <20%. SCR assays were performed at the Shanghai Children's Hospital using the enzymatic method and were analyzed with a 7,180 automatic analyzer (Hitachi High-Tech Science Systems Corporation, Tokyo, Japan). The calibration range was from 3 to 100 mg/L. SCR assays were performed at the Suzhou Hospital Affiliated to Nanjing Medical University using the enzymatic method and were analyzed with a 7,600 Automatic Analyzer (Hitachi High-Tech Science Systems Corporation, Tokyo, Japan). The calibration range was 0.08-100 mg/L.
Creatinine clearance was calculated using the Schwartz formula as in Eq. 1 (Schwartz et al., 1984): where CLcr represents creatinine clearance, HT (cm) represents height, SCR (umol/L) represents serum creatinine, k is 0.45 for term neonates, and 0.33 for preterm neonates. SCR was standardized to the enzymatic method (SCR*) if the Jeff method (SCR † ) was employed in the external model by Eq. 2 (Srivastava et al., 2009) If the method was not clarified in the report, the enzymatic method was used.
The reported population pharmacokinetic model was reconstructed based on information extracted from the original articles. The NONMEM code for each model was determined by a double check. The predictabilities of all candidate models were externally evaluated by prediction-and simulation-based diagnosis and Bayesian forecasting (Zhao et al., 2016;Mao et al., 2018;Cai et al., 2020).

Prediction-Based Diagnostics
The predicted population concentrations (PRED) were estimated and compared with the corresponding observations (OBS) by estimating the relative prediction error (PE%) using Eq. 3: The median prediction error (MDPE) was used to evaluate predictive accuracy, whereas the median absolute prediction error (MAPE) was used to evaluate predictive precision. F 20 and F 30 were also calculated as combination indexes of both accuracy and precision, and indicate the percentage of PE that fell within the ±20% and ±30% ranges, respectively. When the standards of MDPE ≤ ± 15%, MAPE ≤30%, F 20 > 35%, and F 30 > 50% were reached, the model could be determined as satisfactory and clinically acceptable.

Simulation-Based Diagnostics
A prediction-and variability-corrected VPC (pvcVPC) (Bergstrand et al., 2011) was conducted for simulation-based diagnostics. The pvcVPC takes into account typical population predictions and typical population variabilities compared with the traditional VPC, accounting for the different expected variabilities within individuals. The pvcVPC was conducted with 1,000 simulated datasets generated using the models to be evaluated. The pvcVPC was performed using the Perl speak NONMEM toolkit (PsN, version 4.7.0).
Maximum a posteriori Bayesian (MAPB) forecasting was conducted to assess the influence of prior observations on model predictability. Patients with ≥1, 2, three observations were included in the analysis for Bayesian forecasting using zero, one and two previous observations, respectively. For a patient, the individual prediction (IPRED) of the third observation was predicted using the first and second observations, the second observation was predicted using the first observation, and then compared with the corresponding observations. The relative differences denoted by the individual prediction error (IPE%) were calculated using Eq. 4 below: To evaluate the predictability of the candidate models when prior information is increased, the standards of an IPE% (MDIPE) ≤ ± 15%, an IPE% (MAIPE) ≤ 30%, an IF 20 > 35%, and an IF 30 > 50% were used for MAPB forecasting.

The Impact of Modeling Approaches
Different modeling strategies were used in previous studies, which may affect the predictive performance of the model. To explore the impact of these different modeling strategies, we evaluated the predictability of various structural models and covariate models employed in previous studies. The assessment methods include the aforementioned predictionbased diagnostics and Bayesian forecasting methods.

Review of Published popPK Analysis on Vancomycin
After a systematic literature search, 18 neonatal vancomycin models (Seay et al., 1994;Grimsley and Thomson, 1999;Capparelli et al., 2001;Kimura et al., 2004;Mulla and Pooboni, 2005;Marqués Minñana et al., 2010;Mehrotra et al., 2012;Zhao et al., 2013b;Frymoyer et al., 2014;Li et al., 2017;Sheng et al., 2017;Song et al., 2017;Chen et al., 2018;Li et al., 2018;Moffett et al., 2018;Colin et al., 2019;Germovsek et al., 2019;Moffett et al., 2019) were included in this study. The literature search procedure is shown in Supplementary Figure S1. Of the enrolled studies, six were from UNITED STATES, three from the UNITED KINGDOM, five from China, one each from Japan and Spain, and one study enrolled patients from France, Greece, France, and Malaysia. Only 10 models described the analytical method used for the determination of SCR.
Most studies employed sparse sampling strategies. Six models were established with a two-compartment model (2CMT), whereas 12 models were established with a one-compartment (1CMT) model.
Weight, age, and renal function were the most important covariates for clearance identified in the previous studies. Maturation models were employed in 11 studies and could be described by Eq. 5 (Holford et al., 2013): where CLstd represents the baseline clearance, Fsize refers to the body size factor, and Fmat refers to the maturation factor. Weight (current weight and birth weight) and age (postmenstrual age, PMA; postnatal age, PNA, and GA) were regarded as the main factors for body size and maturation, respectively. The current weight was used in most of the reported studies. PMA was chosen over GA and PNA as it presented the most parsimonious way to account for both antenatal and postnatal maturation, which can be incorporated as a sigmoid Emax model and asymptotic exponential model, as shown in Eq. 6 and Eq. 7 (Salman et al., 2019). Among the included studies, seven out of 18 models applied the sigmoid Emax model: where TM 50 is the value of PMA when clearance maturation reaches 50% of adults; Hill is the slope parameter for the sigmoid E max maturation model.
Renal function was often presented by SCR in reported studies and was included in nonlinear manner.
The characteristics of each study are summarized in Table 1, and the information on population pharmacokinetic models were shown in Table 2.

External Evaluation Cohort
The study population consisted of 171 neonates in whom there were 319 assessments of vancomycin levels. Of these, 80 neonates with 165 vancomycin levels were from Shanghai Children's Hospital (SH), and 91 neonates with 154    Frontiers in Pharmacology | www.frontiersin.org March 2021 | Volume 12 | Article 623907 vancomycin levels were from Suzhou Hospital Affiliated to Nanjing Medical University (SZ). The proportion of preterm neonates was larger in the SZ group than in the SH group. The patient demographics of both cohorts are shown in Table 3.

Prediction-Based Diagnostics
There were large differences in predictability of the different models. As can be seen from the results of the prediction-based diagnostics shown in Figure 1 and Supplementary  Table S1, no model satisfied the standards of MDPE ≤ ± 15%, MAPE ≤30%, F 20 > 35%, and F 30 > 50%. Nine models showed good predictive accuracy, with an MDPE of less than ±15%. However, MAPE was more than 30% in all models, indicating a poor predictive precision for all models. Of note, the model reported by Moffett et al. (2019) (Moffett et al., 2019) reached the criteria of F 20 > 30% and F 30 > 45%, showing better accuracy and precision of predictability than other models. The boxplot for PE% for each model is shown in Figure 1.
FIGURE 2 | Box plots of individual prediction error (IPE%) in external data with Bayesian forecasting for published population pharmacokinetic models in different scenarios (0 represents predictions without prior information and one to two represents with prior one to two observations, respectively).

Simulation-Based Diagnostics
In the case of simulation-based diagnosis, the pvcVPC showed significant differences between observations and model simulations in all reported studies (Supplementary Figure S2). A clear tendency of either over-or under-prediction was observed for all models.

Bayesian Forecasting
In total, 171 neonates and 171 observations, 112 neonates and 224 observations, 24 neonates and 72 observations with zero, one, two previous samples, respectively, were included in the Bayesian forecasting. After Bayesian forecasting with one or two prior observations for all models, the mean values of MDPE, MAPE, F 20 , and F 30 compared favorably with the prediction-based diagnostics, indicating that the predictive performances had improved, as shown in Supplementary Table S2. Furthermore, 12 of the published models showed a median IPE <20%, a median absolute IPE <30%, an IF 20 > 35%, and an IF 30 > 50% after Bayesian forecasting with one or two prior observations. The box plots for predictability are shown in Figure 2, and the predicted indices are listed in Supplementary Table S2.

The Impact of Modeling Approaches
The structural model employed in the previous studies included the 1CMT model or 2CMT model. We evaluated the predictive performance of these two models by establishing 1CMT and 2CMT base models using the external data. Because trough concentrations could not fully describe a twocompartment model, the volume of distribution of the central compartment was fixed at 1.27 L, and the inter-compartmental clearance was fixed at 1.161 L/h in the 2CMT model, according to the study by Song et al. (Song et al., 2017). The MDPE of the 2CMT model was less than that of the 1CMT model (3.49% vs. 10.33%), indicating that the predictive accuracy was better for the 2CMT model ( Figure 3 and Supplementary Table S3).
To assess the predictive performance of the various covariate models, we developed models with the identified covariates (WT, PMA, and SCR) and corresponding formulas (a maturation model, a nonlinear model, and a linear model) based on the 1CMT or 2CMT structural models. The maturation model was used by Eq. 6, nonlinear model used by Eq. 8, and linear model used by Eq. 9.
After incorporating body size into the maturation model, the F 20 and F 30 were improved significantly compared to the base 1CMT model (F 20 : 27.65% vs. 14.12% and F 30 : 44.71% vs. 19.41%), or that of the 2CMT model (F 20 : 30.0% vs. 14.12% and F 30 : 43.53% vs. 19.41%). Moreover, the predictive performance of the model with WT and SCR included in a nonlinear fashion was much better than in the model where they were included in a linear fashion (Figure 3 and Supplementary Table S2). With one prior observation, the IF 20 and IF 30 values of the base model after Bayesian forecasting were all >35% and 50%, respectively, demonstrating an obvious improvement using Bayesian forecasting. Box plots for the predictive performance in each model are presented in Figure 3.

DISCUSSION
This study performed a comprehensive external evaluation of the published population pharmacokinetic models of vancomycin in neonates. Based on prediction-and simulation-based diagnostics, none of the published models had a good predictive performance according to pre-specified standards. However, after Bayesian forecasting with one or two prior measurements of vancomycin levels, the predicted performance improved significantly. This finding is consistent with external evaluation studies of other antibiotics such as rifampicin, voriconazole, and tobramycin (Cheng et al., 2020), and immunosuppressive agents (Zhao et al., 2016;Mao et al., 2018;Cai et al., 2020).
Body size is a pivotal index for the CL and V of vancomycin in neonates. Comparing different covariate modeling approaches, nonlinear models, especially the maturation model, showed much better predictability than the linear model. When the maturation model was adopted, F 20 and F 30 improved by 30%-50% compared with the linear model.
For drugs with narrow therapeutic windows, weight-based dosing is most commonly used for neonates because it is easy to perform. However, some studies have reported that there are adverse drug reactions related to weight-based dosage regimens for children, especially for drugs with narrow therapeutic ranges, leading to ineffective treatment outcomes and even fatalities in some cases (Koren et al., 1988;Konstan et al., 1991;Back et al., 2019).
As body weight does not fully describe organ maturity, the maturation model could better explain the physiological status of early, slower growth and subsequent faster growth in neonates (Back et al., 2019). It also allows for a quantification of the relationship between the mass/ structure of organs and size (Fsize), which is known to exhibit a nonlinear pattern of growth in neonates (Holford et al., 2013) (Back et al., 2019). Moreover, published population PK analyses that included body size in the maturation models report better forecasting and better clinical use (Andrews et al., 2018;Béranger et al., 2019).
Renal function is also a very important factor affecting the pharmacokinetics of vancomycin, since vancomycin is mainly eliminated via the kidney. Creatinine clearance is used as index to describe renal function in adult patients; however, in neonates, SCR levels are a more reliable indicator of renal function, which is consistent with the findings of previous population PK studies. It has previously been shown that incorporating SCR in CL in a nonlinear fashion is better than incorporating it in a linear fashion. This finding also supports the fact that renal function matures in a nonlinear manner in neonates.
Our study has several limitations. As mentioned previously, the creatinine determination method (Jeff and enzymatic method) has been shown to have a large impact on external predictability (Zhao et al., 2013a). In this study, although we used correction equations to reduce variation between the two methods, several of the previous studies did not clearly state the method used for creatinine determination; therefore, this should be noted in future studies. Moreover, only peak and trough data were collected, and the comparison between the 1CMT and 2CMT models was not fully assessed and so may require further investigation.

CONCLUSION
Based on our study, the published models performed poorly in prediction-based and simulation-based diagnostics. The maturation model based on weight, age, and nonlinear incorporation of SCR had better predictability than other modeling approaches. Moreover, the Bayesian method significantly improved the predictive performance of the published models, and could thus play an important role in vancomycin dosing recommendations and guiding clinical practice.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Shanghai Children's Hospital and Suzhou Hospital Affiliated to Nanjing Medical University. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
Y-XL, ZJ, J-JL, and Z-LL participated in research design. J-JL and Z-LL contributed to acquisition of the evaluation data. Y-XL and ZJ performed the research and analyzed the data, W-JN reviewed the model selected and information verification. Y-XL, ZJ, and HW drafted the manuscript. All authors revised and approved the manuscript.

FUNDING
The work was funded by the Wu Jieping Medical Foundation (2019). The project number is 320.6750.19090-1.