Intravenous (IV) vancomycin has long been the first-line treatment for severe bacterial infections, especially in cases involving Staphylococci species [1]. Despite its well-established efficacy, vancomycin has a narrow therapeutic index and is commonly associated with acute kidney injury (AKI), especially at high levels of exposure [2]. It was shown that even small acute increases in serum creatinine (SCr) could be detrimental to long-term survival in critically ill patients [3]. Therapeutic drug monitoring (TDM) for vancomycin is a recommended practice to balance efficacy and the risk of AKI. This involves the monitoring of the systemic serum vancomycin concentration (VC_s) over time after drug administration and subsequent adjustments of vancomycin dosage as necessary.

Vancomycin is eliminated from the systemic circulation primarily through glomerular filtration in the kidneys. Thus, glomerular filtration rate (GFR) is closely correlated with vancomycin clearance (CL), which is the main factor affecting VC_s [4]. Since GFR is clinically estimated by creatinine clearance, the major determinants of creatinine clearance, including body size and SCr, are among the major covariates of CL among patients from all age groups [5,6]. To improve the prediction of VC_s, the pharmacokinetics of vancomycin has been widely studied to understand the mathematical relationship between CL and these covariates [7-11].

In the neonatal intensive care unit (NICU), the pharmacokinetics of vancomycin among neonates is highly variable due to dynamic patient conditions and interventions [11]. Moreover, for neonates, it is necessary to account for the maturation of renal function, a process unique to the neonatal population that occurs over the first weeks to months postpartum and is associated with postmenstrual age (PMA) [12]. These conditions put NICU patients at a higher risk for suboptimal therapeutic effects of vancomycin and AKI, making accurate TDM of vancomycin indispensable in this population.

A steady-state area under the curve of the VC_s-time profile (AUC) over 24 hours (AUC₂₄) to minimum inhibitory concentration (MIC) ratio (AUC₂₄/MIC) of ≥400 hours has been advocated as the primary predictor of vancomycin efficacy [13]. Nevertheless, since AUC estimation requires measuring multiple VC_s values, which is often impractical in the clinical setting, the American Society of Health-System Pharmacists, the Infectious Diseases Society of America, and the Society of Infectious Diseases Pharmacists published a consensus report in 2009 recommending the steady-state trough VC_s (VC_s,ss,trough) as a surrogate marker for the AUC target (assuming MIC at 1 mg/L) [13]. However, data on the efficacy and safety profile with this trough-based approach are lacking [14]. On the other hand, there is further evidence supporting AUC₂₄/MIC as the pharmacokinetic target. The requirement of multiple VC_s measurements could also be resolved by employing the Bayesian approach as supported by recent research [15,16].

In response, the guideline was updated in 2020 jointly by the 3 societies publishing the 2009 report, together with the Pediatric Infectious Diseases Society, giving new recommendations. First, VC_s,ss,trough is no longer recommended as a pharmacokinetic target; dose optimization should instead target an AUC₂₄/MIC of 400 to 600 hours. Second, the preferred method to estimate individual AUC is to apply Bayesian estimation using 1 trough VC_s (VC_s,trough, presteady-state or steady-state trough) and preferably 1 peak VC_s (VC_s,peak, presteady-state or steady-state peak), based on a population pharmacokinetic (popPK) model for vancomycin. Third, a less preferred method to calculate individual AUC is to use the first-order equations on a set of measured VC_s,ss,trough and steady-state peak VC_s (VC_s,ss,peak) values [14]. Recommendations for initial dosing were also revised. A popPK model–based estimation of individual AUC is preferred over using a universal weight-based dosing scheme [17]. As AUC becomes the basis of dose optimization, ensuring a reliable approach for AUC estimation is a prerequisite of dose optimality. Table 1 summarizes the approaches to AUC estimation and hence dose optimization used in this text.

The rapid development in digital health has made this study and the proposed clinical improvements possible in multiple ways. They are elaborated in the following paragraphs.

To keep up with the current standard of treatment and given the large variability in NICU patients, separate popPK analyses for vancomycin are required for the local NICU population [18]. However, prospective data collection is often costly and burdensome in the clinical environment, while the unstructured collection of retrospective data is prone to errors. Fortunately, as digital records are becoming vital on the clinical frontline, electronic health records (EHRs) now present extractable information for data analyses [19]. It is now feasible to consolidate data retrieved from multiple EHR sources to reconcile a data set suitable for popPK analyses [20].

Establishing a popPK model is the first step to upgrade the TDM practice for IV vancomycin in the local NICU population. To maximize the utility of the popPK model, it is necessary to enable Bayesian estimation for accurate estimations of individual AUC [14]. Unlike conventional strategies to individual dose optimization by equations and nomograms, which can be carried out manually, Bayesian estimation requires numerical approximation processes that can only be performed digitally using computers.

Putting popPK model–based Bayesian estimation into clinical practice is difficult because most clinicians are not experts in this area. To tackle this, a fully automated web-based interface incorporating a popPK model, a numerical approximation solution to Bayesian estimation, and algorithms for dose optimization would be an ideal tool for clinical use. In contrast with a client-based interface, a web-based interface (1) allows remote access with various browser-enabled devices, including clinical computer workstations, tablets, and smartphones; (2) saves installation issues; and (3) is easier to maintain. Such an interface is designed to guide and validate necessary inputs from clinicians, followed by suggestions of dosing regimens, which are expected to help clinicians decide the optimal treatment plan that can enhance clinical outcomes.

The use of a model-based dosing interface is in its pilot stage in Hong Kong. Neonatal vancomycin is among the first drugs being investigated. Experiences gained in this study are expected to improve IV vancomycin treatment significantly and, perhaps more importantly, lay the foundation for the extraction of popPK data from EHRs and web-based dose optimization interfaces for other drugs with narrow therapeutic indices. In support of its implementation, this study was conducted in local medical institutions to develop the popPK model of vancomycin for NICU patients using real-world data from EHR resources. Besides, a previously reported framework of a web-based interface performing Bayesian estimation and individual dose optimization for the use of high-dose methotrexate in local institutions will be adopted to create the dosing interface for neonatal IV vancomycin.

The study data set consists of all Chinese patients within 1 year of postnatal age (PNA) admitted to the NICUs of Prince of Wales Hospital and Queen Elizabeth Hospital in Hong Kong between January 1, 2016, and December 31, 2017. Each potential subject had to be prescribed IV vancomycin, and have at least one VC_s measurement and one SCr measurement in order to be eligible. Subjects with major congenital heart diseases were excluded. Subjects with vancomycin initiated within 7 days of birth were also excluded due to the variable effects of maternal creatinine on the estimation of neonatal renal function. Eligible subjects were identified through the Clinical Data Analysis and Reporting System (CDARS), a database developed and maintained by the Hong Kong Hospital Authority for audit and research purposes. Data of selected subjects were then collected from several in-house EHR platforms, namely the Clinical Information System (CIS), In-Patient Medication Order Entry (IPMOE), electronic Patient Record (ePR), and, whenever necessary, original copies of medical charts. Ethical approval was obtained from the Joint Chinese University of Hong Kong-New Territories East Cluster Clinical Research Ethics Committee (reference number: 2018.094) and Kowloon Central Cluster/Kowloon East Cluster Research Ethics Committee (reference number: KC/KE-18-0096/ER-1) for data collection. Parental consent was not required due to the anonymized and retrospective nature of data collection.

Constant data items collected were sex, birth weight, gestation age, and date of birth. Time-dependent measurements included body weight, VC_s, and SCr. Dosing records were collected in terms of the dose administered and the infusion rate (assuming constant rate). The date and time tags of all time-dependent measurements and dosing records were also collected. Dosing records were available from IPMOE, while other data items were collected from CIS and ePR.

All VC_s records collected over 7 days after the start of the last infusion of vancomycin were removed. All subjects started on vancomycin within the first 7 days of birth were also removed. At each unique time tag of each subject, PMA was calculated as the time difference between the tag and the estimated first day of the last menstrual period of the subject’s mother. SCr was imputed as the previous or next available value, whichever was closer in time. Body weight was imputed by linearly interpolating and extrapolating available values. All VC_s values measured during infusion were removed. SCr values below the lower limit of quantification (LLOQ) (ie, 15 μmol/L) were replaced by 7.5 μmol/L. VC_s values measured had an LLOQ of 1 mg/L, and records below the LLOQ (below the limit of quantification [BLQ]) were flagged.

A popPK model adopts the structure of a nonlinear mixed-effect model [21,22]. A 2-compartment infusion model with first-order elimination was applied, for which the pharmacokinetic parameters CL, central volume (V_c), intercompartmental clearance (Q), and peripheral volume (V_p) of vancomycin were defined [9]. The between-subject variability (BSV) in CL and the between-occasion variability (BOV, variability in CL in the same subject between episodes) were expressed in terms of the coefficient of variance (CV) (ie, CVCL and CVCL_BOV, respectively). Both CVCL and CVCL_BOV were assumed to follow the log-normal distribution [23]. Residual unexplained variability was described by a combined proportional-additive error model [24].

In building the pharmacokinetic parameter model, allometric scaling was applied to describe the association of CL, V_c, Q, and V_p against body weight using the power function, with fixed exponents of 0.75 and 1, respectively [25]. This was tested against freely estimated exponents (one for CL and Q, and another for V_c and V_p) using the likelihood ratio test. The maturation of renal function was described as a function of PMA and tested against the linear, exponential, first-order, and Hill functions [26-28]. The function that returned the best goodness of fit was chosen. The renal function with respect to SCr was described using the power model [29].

Parameter estimation was executed with NONMEM version 7.4 (Icon plc) using first-order conditional estimation with interaction [30]. BLQ data were handled using the M3 method [31]. Perl-speaks-NONMEM was used to coordinate NONMEM execution [32]. Residual plots, predictive plots, and a prediction-corrected visual predictive check (pcVPC) were generated [33]. Bootstrapping using 1000 resampled data sets was performed to assess the stability of parameter estimates [34]. R and the R package ggplot2 were used for graphics generation [35,36].

Model-based approaches to dose optimization rely on the estimation of pharmacokinetic parameters for a subject based on the verified popPK model as described above. The set of pharmacokinetic parameter estimates are then used to approximate the AUC distributions at different doses, such that the dose at which the probability of attaining an AUC₂₄/MIC of 400 to 600 hours (probability of target attainment [PTA]) is maximized (ie, the optimal dose) can be identified by numerical approximation. Practically, the empirical approach helps decide the initial dose, while the Bayesian approach informs dose adjustments afterwards (see Table 1).

The web-based dosing interface in this study is designed to perform the model-based approach to dose optimization. The framework of the interface was replicated from that reported in a previous study for the dose adjustments of single-dose high-dose methotrexate in the pediatric population [37]. On performing the Bayesian approach, the interface demonstrated the ability to generate individual estimates of AUC identical to and more efficiently than NONMEM. The interface was modified to adapt to the popPK model for IV vancomycin estimated and verified in this study, enable empirical dose (the first dose) suggestion, and allow dose optimization at various dosing intervals.

One VC_s value was measured 7 days after the start of the last infusion and thus removed. Forty-five subjects with PNA <7 days when vancomycin was first started were also removed. Data extraction from the EHRs and data exclusion resulted in a final data set consisting of 207 patients and a total of 689 VC_s measurements. The demographics are detailed in Table 2. The time profile of observed VC_s values is shown in Figure 1.

In the final model, a fixed exponent is used for the power functions of the body weight effect and a Hill function is used to describe the PMA-CL relationship. The final model has a minimum objective function value (OFV) of 2329.272. Allowing freely estimated exponents for body weight functions on pharmacokinetic parameters only led to statistically insignificant improvements in goodness of fit (χ² value approximates change in OFV [dOFV]=−5.583, P=.06, at a degree of freedom of 2). Replacing the Hill function with a linear, exponential, or first-order function resulted in worsened goodness of fit (dOFV=+58.440, +83.391, and +140.645, respectively). None of these alternative models led to significantly better goodness of fit.

Parameter estimates for the final model are shown in Table 3. They are very close to the bootstrap means and well within the bootstrap CI and have a condition number of 226 (which is within the usual reference limit of 1000). Accounting for both CVCL and CVCL_BOV, the effects of body weight, PMA, and SCr alone, and all combined explained 43%, 63%, 54%, and 82% of the variability in CL, respectively. Besides, the overall shrinkage of random effects in CL is estimated to be 16.9%, which is within the acceptable range. Predictive and residual plots of the final model are available in Figure 2. The prediction-corrected visual predictive check is shown in Figure 3.

Based on the validated popPK model, the PTAs of different dosing approaches for the subjects in the data set were approximated. The graph on the left in Figure 4 shows that dose adjustments by the steady-state trough approach result in only half of the subjects achieving a PTA over 50%, which is only slightly improved when compared to maintaining the initial doses given. It is also outperformed by the model-based approaches, namely, the empirical approach and the Bayesian approach, which reliably raise the PTA to above 75% for most subjects. Meanwhile, the graph on the right in Figure 4 shows some extreme dose adjustments with the steady-state trough approach when compared to the Bayesian approach, indicating overcorrection of doses with the trough approach without achieving a better PTA profile.

A composite screenshot of the developed interface is available in Figure 5. Detailed screenshots of the developed interface are available in Multimedia Appendix 1. The top panel is always displayed and allows the user to navigate different steps using the interface. By clicking a tab, the corresponding panel will be displayed below the top panel. Step 1 requires user inputs to estimate individual parameters. The user may choose between the model-based approaches (the empirical or Bayesian approach), depending on whether VC_s data are available. If the latter is chosen, then apart from the current body weight, PMA, and SCr, the user also needs to input previous doses administered, measured VC_s, and previous body weight, PMA, and SCr. In step 2, the user may specify the desired therapeutic targets, which defaults to an AUC₂₄/MIC of 400 to 600 hours without constraints by VC_s,ss,trough and VC_s,ss,peak. Step 3 allows the user to specify the range of doses and dosing intervals allowed during optimization, which are by default set according to the usual practices of the hospitals using the interface. In most cases, accepting the defaults for steps 2 and 3 suffices. After inputting the required data, the results of individual dose optimization will be generated in the “optimization” tab, which suggests the dose (and dosing intervals) required to maximize PTA and the graphical illustrations of the steady-state VC_s profile and the expected probability distribution of AUC₂₄.

Real-world data from EHRs were successfully used to develop a popPK model of IV vancomycin for the local NICU population. Based on a previously published dosing interface for high-dose methotrexate, a dosing interface for individual dose optimization of IV vancomycin for the local NICU population was created.

Despite the promising aspects of using EHRs, challenges were present when attempting to reconcile a popPK data set with EHR resources. First, especially when encountering a new EHR source, efforts were required to understand and validate the data structure of the source to ensure the likeliness of generating necessary tables for statistical analyses. To enable popPK analyses, it is essential to ensure that the target information can be reformatted into tables with different row representations (ie, 1 row per, for example, subject, dose, and observed VC_s). Then, since EHRs are primarily archived automatically during clinical operation, there is the issue of unstandardized or ambiguous inputs, especially for manual fields, because different clinics may have different logging practices. For instance, laboratories may run assays with different LLOQs, which could be logged onto the EHR systems using various syntaxes. Other problems encountered were suspected duplicated or missing records. For example, detectable VC_s measured before the first recorded dose or 7 days after the last dose in subjects with normal renal function may indicate missing dosing records.

A major limitation of using EHR data following the above issues is that data errors and ambiguity are often untraceable. To ensure the robustness of the final data set submitted for popPK analyses, it is crucial to remove problematic data that cannot be clarified from the EHR sources while keeping an eye on the possible risk of causing biased estimates (eg, censored data that are missing not at random).

While having to identify unsalvageable data is a downside, using EHRs is a convenient way to obtain a useful volume of data. Under the hectic environment of hospital wards, it is often difficult for clinicians to cater to the collection of study data. Making use of EHRs can ease the data collection process by minimizing the clinical workforce required. Moreover, since most EHR fields are already standardized, organized, and validated to a certain extent, typographical errors are less of a concern when extracting information from EHR sources.

The covariates can explain a significant proportion of BSV as expected. Diagnostic plots and the prediction-corrected visual predictive check show good predictive performance. The agreement between final parameter estimates against bootstrapping results and the relatively small condition number demonstrates the stability of the estimates. The choices of parameter-covariate relationships in the final model structures and the resultant parameter estimates in this study generally agree with previously reported models [26,28,29]. These positive results of evaluations help establish the validity of the model for implementation into the dose individualization interface.

It is anticipated that the implementation of the developed interface can bring about several improvements to the current practice of administering IV vancomycin to treat severe infections in critically ill neonates. First, with the support of the popPK model developed, the interface can estimate individual AUC more accurately and enhance the optimality of the recommended initial dose (with the empirical approach) and maintenance dose (with the Bayesian approach). The recommended dose is also adaptive to significant changes in individual vancomycin PK due to variations in body weight, PMA, and renal functions during treatment. Moreover, with the Bayesian approach, presteady-state VC_s is also usable for estimation, such that waiting until the release of a steady-state VC_s measurement result for dose adjustment is no longer required. Together, these advantages promote the PTA profile and shorten the time to achieve the pharmacokinetic target by reducing the number of dose adjustments required. This is, in turn, expected to improve the treatment outcomes by promoting recovery while mitigating the risk of developing AKI. Apart from that, implementing the interface eliminates the need for manual calculation and thus reduces the risks of arithmetic errors in dose adjustments. The interface is also designed in a user-friendly and foolproof manner to ease its application by clinicians. Furthermore, the interface is developed using open-source software such that accessibility is guaranteed and licensing costs can be saved.

Based on a data set reconciled from real-world data extracted from multiple EHR sources, a popPK model of IV vancomycin has been developed and verified for the local NICU population. Based on the verified model and adoption of a previously published framework, a web-based dosing interface has been built to apply model-based approaches to individual AUC estimation and dose optimization of IV vancomycin. The developed interface is expected to improve clinical outcomes of the treatment of severe infections compared with previously adopted approaches, namely, the weight-based approach for initial dosing and the trough-based approach for dose adjustments. A cohort study will be performed later to show the superiority of using the interface compared with the previous approaches in terms of clinical outcomes. The experiences gained in this study will be valuable for the future use of the data collected from EHR sources for popPK analyses and the development of similar interfaces for other drug entities with narrow therapeutic indices.