This study used the same dataset as [10], consisting of 5483 patients with 8689 encounters in total, which was obtained from an EHR data warehouse containing information regarding encounter-level administrative data collected from Memorial Hermann Health System (MHHS), a large health care system based in Houston, Texas, United States. Patients aged more than 18 years who received at least 1 dose of vancomycin during the period from August 2019 to March 2020 were included in the study. This study excluded patients who underwent renal replacement therapy and patients with inappropriate timing of vancomycin levels. Deidentification of the extracted cohort was conducted to protect the privacy and security of the patient data. Similar to a study by Nigo et al [10], the model was trained and evaluated only using encounters with at least 1 serum vancomycin level measured after the first vancomycin dose.

The PKRNN-2CM model is designed to use variables that are routinely collected in EHR. The EHR data for this study were extracted at the encounter level, including 30 laboratory tests, 5 vital signs, 324 types of medications, and demographic information. For medication features, we represented drugs using their simple generic names rather than pharmacologic classes. This choice was based on our prior work [17], which showed that finer clinical granularity can improve model performance in high-dimensional EHR models. Using generic names preserves clinically meaningful distinctions without introducing the excessive detail associated with brand names or dose-specific entries. In addition, medication variables were encoded using embedding representations, which help mitigate data sparsity while retaining relevant clinical information. The variables included in this study are commonly measured in clinical practice, particularly complete blood counts and basic metabolic panels.

Out of the 40 continuous variables, some laboratory tests, such as serum calcium level, phosphorus, total bilirubin, total protein, serum magnesium, and nucleated red blood cells, may occasionally have missing values. However, these laboratory tests are typically measured at least once during hospitalization when bacterial infections are suspected and vancomycin is used. Other variables included in the study are demographic data and vital signs, which are also routinely collected in hospitals. The 8 categorical variables are used to create an embedding vector for medications. In the same manner as described by Nigo et al [10], for each encounter, the start time was determined by the earliest record time, while the end time was determined by the time stamp of the last vancomycin concentration. Vancomycin and other medication administration, vancomycin and other laboratory measurements, and the end of the day were used as events to update the state of the PKRNN-2CM model. Within the MHHS dataset, the infusion duration of vancomycin is dosage-dependent: 1000 mg is infused over 1 hour, 1001 to 1500 mg over 1.5 hours, 1501 to 2000 mg over 2 hours, and doses exceeding 2000 mg over 2.5 hours. For simplicity, a uniform infusion rate of 1 g/h was adopted. This standard rate is used in MHHS as a default rate to prevent acute reactions such as red-man syndrome. While individual infusion rates were not available in the MHHS dataset, the PKRNN-2CM model can be adjusted to incorporate various infusion rates, making it adaptable to different clinical settings. We used forward filling and mean imputation to handle missing covariates. Continuous variables were standardized using z scores.

The PKRNN-2CM model, similar to the PKRNN-1CM model [10], is an autoregressive RNN built on EHR data, requiring minimal manual input, and providing real-time concentration predictions by automatically adjusting to reflect subtle changes. The PKRNN-1CM model consists of 3 components that enable the vancomycin predictions to be evaluated at any time point. First, a code embedding layer is used to convert information from the EHR into the PKRNN-1CM model for every time step. Next, an RNN layer is applied to predict vancomycin elimination rates and compartment volumes, and finally, based on the output from the RNN layer, a 1-compartment PK layer uses the PK equation, which is an ODE, to compute the predicted vancomycin concentration. For the PKRNN-2CM model, we extend the 1-compartment PK model to a 2-compartment PK model. A schematic representation of the structure of the PKRNN-2CM model can be found in Figure 1. The input to the model is the time-series EHR data after data preprocessing, which contains 8-dimensional embedding derived from categorical data that includes all 911 different medication codes of 324 medications (Multimedia Appendix 1 shows the generic names of medications used in the model) and 40D continuous data that includes time, time difference, vancomycin dose, vancomycin serum concentration measurements, laboratory tests, and demographic information (Multimedia Appendix 1 shows all the continuous variables included in the model). For every time step, the RNN layer predicts 4 parameters η₁, η₂, η₃, and η₄ [18] that are related to the PK parameters. After that, the PK layer based on the 2-compartment PK model uses an ODE to calculate the estimated PK parameters CL (clearance for the central compartment), Q (intercompartmental clearance), V₁ (volume of the central compartment), and V₂ (volume of the peripheral compartment), then uses the estimated PK parameters to predict the vancomycin concentration.

PKRNN-2CM differs from PKRNN-1CM in two aspects: (1) it uses additional PK parameters and different initial values; and (2) it uses a first-order ODE system with 2 ODEs to calculate vancomycin concentrations. As Figure 2 shows, there are 2 PK parameters in the PKRNN-1CM model, elimination rate k and volume of distribution V, while the PKRNN-2CM model has 2 clearances, central (CL) and intercompartmental (Q), and 2 volumes, V₁ and V₂ representing the central and peripheral compartments, respectively. Furthermore, empirically we found that the performance improved by removing the constraints of mass conservation as in PKRNN-1CM. Instead, in PKRNN-2CM we make the concentration continuous, which we believe is better because the model is allowed to correct its own prediction error in a smoother way.

Based on a study by Lim et al [18], the 4 PK parameters CL, Q, V₁, and V₂ can be described by 4 related parameters η₁, η₂, η₃, and η₄ that follow a multivariate normal distribution (MVN). The RNN layer in the PKRNN-2CM predicts η₁, η₂, η₃, and η₄ where the initial values are determined according to Lim et al’s [18] study:

as opposed to PKRNN-1CM which uses RNN to predict k and V directly. In the PK layer, η₁, η₂, η₃, and η₄ are used to calculate PK parameters following the equations below, where the content in brackets indicates the unit:

where CrCL is the creatinine clearance calculated from the Cockcroft-Gault equation:

The vancomycin concentration can then be calculated using the ODE system below:

Where C_A and C_B are the concentrations in central and peripheral compartment (mg/L).

A closed-form solution to this ODE can be found in Multimedia Appendix 2.

The simulation was designed to be a valuable tool to guide the development and deployment of our PKRNN-2CM model by bridging the gap between the constraints of real-world data and the requirements of a comprehensive and reliable predictive system. Through this approach, we can address the limitations of sparse real-world datasets and gain deeper insights into the behavior of the PKRNN-2CM model in comparison to the PKRNN-1CM model. To ensure that the simulated dataset was comparable to the original real-world dataset, the simulation was based on actual patient information such as medications and laboratory results. As much patient information as possible was used in the simulation to guarantee the similarity between the simulated datasets and the original real-world dataset, with the only simulated data being the vancomycin measurements. We first trained a PKRNN-2CM model using real-world data, resulting in the PKRNN-2CM–derived model, which was then used to simulate vancomycin concentrations. This PKRNN-2CM–derived model was set as the “ground truth” to generate measurements for the simulated datasets. We then applied either the PKRNN-1CM or PKRNN-2CM model to the simulated data to learn the parameters, resulting in the PKRNN-1CM–estimated or PKRNN-2CM–estimated models. Although the PKRNN-2CM–derived and PKRNN-2CM–estimated models share the same architecture, they are trained independently without exchanging information. In summary, the PKRNN-2CM–derived and PKRNN-2CM–estimated models operate independently, with the only connection being that the PKRNN-2CM–estimated model uses the predicted concentrations from the PKRNN-2CM–derived model as simulated measurements.

The simulation used in this study involves the use of vancomycin concentrations predicted by the PKRNN-2CM–derived model as labels for testing the performance of the PKRNN-1CM–estimated and PKRNN-2CM–estimated models based on different types of simulation.

Figure 3 is a schematic diagram showing the detailed simulation and evaluation process. In step 1, the real-world time-series EHR data are loaded into the system. In step 2, the PKRNN-2CM–derived model is trained on real data. In step 3, simulation is set to false and model performance is assessed through the calculation of root mean square error (RMSE), comparing observed and predicted concentrations. In step 4, simulation is set to true, and the real data are removed, and predicted concentrations from the PKRNN-2CM–derived model are used to generate a simulated dataset based on the simulation options. In step 5, the PKRNN-1CM–estimated or PKRNN-2CM–estimated model is trained on the simulated data. Step 6 evaluates the estimated model’s performance using RMSE. Finally, a comparative analysis of RMSEs between PKRNN-1CM–estimated and PKRNN-2CM–estimated models measures their efficacy in different scenarios. A detailed simulation process overview can be found in Multimedia Appendix 1.

To establish the framework for our model implementation, we meticulously configured the code embedding and the RNN layers, and also chose training hyperparameters and optimization techniques for optimal performance. In configuring the code embedding layer, categorical data are embedded into 8D vectors. For the initialization of the embedding layer, weights were established using a Gaussian distribution. Each time step involves the input of a 48D vector (40+8) into the RNN layer, comprising embedded categorical data (8) and normalized continuous data (40). The RNN layer used a single-layer gated recurrent unit (GRU) with a hidden size of 64. The output layer is characterized by a linear layer of size (64, 4), mapping the GRU’s hidden layer to the parameters η₁, η₂, η₃, and η₄ at each time step. For model training, the Adamax optimizer is used with a learning rate and weight decay set at 1 × 10⁻² and 0.2, respectively. The training minibatch size is configured at 50, and an early stopping mechanism, with “patience” set to 10, is implemented to mitigate the risk of overfitting. The mean squared error serves as the loss function for training the model. Additionally, 2 regularization mechanisms were incorporated. First, a penalty was imposed on the deviation of the predicted parameters η₁, η₂, η₃, and η₄ from the prior multivariate Gaussian distribution. The second regularization term uses the L2 norm of the first-order difference to discourage abrupt changes in the output of the RNN layer.

For the simulation implementation, several datasets were generated from the PKRNN-2CM–derived model based on different simulation options. In the “peak” datasets, a measurement was simulated after the vancomycin administration. The time intervals between administration and measurement were dose-dependent; the simulated measurements were made after 2 hours if the dose was less than or equal to 1000 mg, or 3 hours if the dose was greater than 1000 mg. The “trough” datasets simulated each measurement 1 hour before the vancomycin administration.

For consistency and comparability, all models were evaluated using the same dataset and criteria. Depending on patient identification, the data were divided into training, test, and validation sets in a ratio of 50:30:20. The performance of PKRNN-2CM and PKRNN-1CM was compared on both the original MHHS dataset and the simulated datasets using RMSE. A 2-sample 2-tailed t test was used to compare the model performance between PKRNN-1CM and PKRNN-2CM. Additionally, we ran the PKRNN-1CM–estimated or PKRNN-2CM–estimated model with different simulation options to evaluate how the estimated models captured the area under the concentration curve. This study was conducted with Python 3.8 (Python Software Foundation) using the PyTorch 1.9.0 library. All experiments were run on an Nvidia A100 80GB GPU. The training process takes less than 60 seconds per epoch, and evaluation and testing take less than 30 seconds per epoch.

A detailed code repository for reproducing the experiment can be found on GitHub [19].

This study was approved by the institutional review board at The University of Texas Health Science Center at Houston (approval: HSC-MS-19‐1011) and was considered exempt from consent due to the retrospective nature of the study and use of deidentified electronic health record data. All data were handled in compliance with institutional policies and applicable regulations for the protection of patient privacy and confidentiality. No identifiable patient information was included in the analysis or reported in this study

Table 1 summarizes the baseline characteristics of the study cohort. A total of 5483 patients were included in the analysis.

Table 2 shows that the PKRNN-2CM model exhibited better performance compared to the PKRNN-1CM model (a 2-sample 2-tailed t test; P=.01) on the MHHS datasets. A Bayesian vancomycin therapeutic drug monitoring (VTDM) model [18] was implemented as the baseline here, with PK parameters updated after each measurement using maximum a posteriori estimation. The VTDM model performed worse than any of the PKRNN models.

Simulation results were reported using average RMSEs from a common test set. For each dataset simulated from the PKRNN-2CM–derived model, we tested and compared the model performance of the PKRNN-1CM–estimated and PKRNN-2CM–estimated models. Table 3 presents the summarized simulation results.

Table 3 shows the performance of the PKRNN-1CM and PKRNN-2CM on the simulated dataset. Overall, models trained on simulated measurements achieved lower RMSE than those trained on real-world data, and performance improved with increased measurement frequency. Given a fixed measurement budget, combining peak and trough observations consistently outperformed sampling only peak or trough levels, highlighting the benefit of balanced sampling strategies. In most settings, PKRNN-2CM showed slightly lower variability than PKRNN-1CM, indicating more stable performance. The VTDM model performed poorly on the simulation data, likely due to its inflexibility in adapting to the dynamic PK parameters generated by PKRNN-2CM–derived data. The results of VTDM on the simulated data are shown in Multimedia Appendix 1.

Figures 4-6 present 3 randomly selected patient scenarios to show model behavior across different clinical settings. Each figure compares the PKRNN-1CM and PKRNN-2CM models using both real data (top panels) and simulated data (bottom panels). In Figure 4, the PKRNN-2CM model shows a clear advantage, achieving lower RMSEs in both real and simulated settings (real: 27.18 vs 21.32; simulated: 17.54 vs 15.20). In Figure 5, the 2 models perform similarly on real data, while the PKRNN-2CM–estimated model provides better accuracy on the simulated dataset (13.24 vs 11.46). In Figure 6, the PKRNN-1CM model slightly outperforms the PKRNN-2CM model on the real data (14.57 vs 16.19), whereas the PKRNN-2CM–estimated model performs better in the simulated scenario (11.14 vs 9.74).

Across Figures 4-6, predicted area under the concentration-time curves (AUCs) from the PKRNN-1CM model exhibited more abrupt transitions, particularly during the postinfusion exponential decay phase. These arise because the 1CM corrected its trajectory abruptly when new observations deviated from prior predictions. In contrast, the PKRNN-2CM model produced smoother curves by updating its state more gradually as new information became available. The curve behavior reflected the prospective plotting strategy, where earlier predictions were not retroactively adjusted. This emphasized real-time diagnostic interpretability over retrospective smoothing.

Table 4 compares the time-averaged, dose-normalized AUC/minimum inhibitory concentrations (MICs) calculated from different simulation options. The ground truth was derived from the PKRNN-2CM–derived model. The results indicate that the PKRNN-2CM–estimated model can provide a closer estimate of AUC/MIC compared with the PKRNN-1CM–estimated model. Additionally, for both models, predictions from combined peak and trough measurements yield smaller differences from the ground truth compared to peak or trough measurements. Thus, mixing peak and trough measurements can improve the accuracy of the predicted concentration-time curves for both estimated models. As before, given a fixed measurement budget, obtaining both peak and trough levels is superior to sampling measurements only at the peak or trough level across the doses.