This was a retrospective study for PK and PD properties of vancomycin in 542 patients including 40 pediatric patients (<20 years), who received vancomycin intravenous IV) infusion at Severance Hospital in Seoul, Republic of Korea. The data were extracted from electronic medical records. The study encompassed patients who were treated with vancomycin and underwent Therapeutic Drug Monitoring (TDM) services. Patients lacking dosing history and concentration sampling time information were excluded from the analysis. Demographic information, laboratory test results such as creatinine and blood urea nitrogen (BUN), and comorbidity details were collected as potential covariates, as well as drug concentrations with 3 measurements per patient on average for PK analysis and CRP concentrations for PD analysis. PK and PD analyses were performed using NONMEM software version 7.5 (ICON Development Solutions, Dublin, Ireland) and exploratory data analysis was conducted using R software version 4.2.2.

Serum vancomycin levels were assessed using the kinetic interaction of microparticles in solution (KIMS) method on the Roche Cobas c702 analyzer (Roche, Basel, Switzerland). A competitive reaction takes place between the vancomycin-macromolecule conjugate and vancomycin in the serum to bind with the vancomycin antibody on the microparticles. The turbidity induced by the binding of the vancomycin conjugate to the antibody on the microparticles is measured photometrically, and this measurement is inhibited by the presence of vancomycin in the sample. The resulting turbidity is indirectly proportional to the amount of vancomycin present in the sample. The lower limit of quantitation for vancomycin is 4.0 μg/mL, determined as the lowest concentration that meets a total error goal of 20%.

CRP concentrations were determined through a particle-enhanced immunoturbidimetric assay on the Roche Cobas c702 as well. Serum creatinine levels were measured by a rate-blanked compensated kinetic Jaffe method on the Atellica CH 930 Analyzer (Siemens Healthcare Diagnostics, Marburg, Germany).

One- and two-compartment disposition models were tested for basic structural model since sparse samples near peak and trough concentrations were taken, where trough and peak concentrations were concentrations at the start and the end of IV infusion within one dosing interval, respectively. Systemic and inter-compartmental clearance and central and peripheral volume of distribution were parameterized using allometry scaling (N. H. Holford, 1996) as the following equations T V V = θ V × WT / 70 × COV V T V C L = θ CL × WT / 70 0.75 × F ren × F mat × COV CL where WT is body weight, TVV and TVCL denote population typical value of volume of distribution (V) and clearance (CL), and θ V and θ CL denote population median values of V and CL for the subject with WT = 70 kg. The same allometry scaling was also applied to the inter-compartmental clearance (Q) and peripheral volume of distribution (Vp). The exponential model was applied to inter-individual variation of these parameters, and for the residual error, additive, proportional and combined models were compared.

Assuming vancomycin clearance is affected by age-related renal function factor ( F ren ) and clearance maturation factor ( F mat ) associated with postmenstrual age (PMA) of up to 48 weeks (Anderson et al., 2007), covariate model building was conducted using the following model structure as a basis: T V V = θ V × WT / 70 × COV V T V C L = θ CL × WT / 70 0.75 × F ren × F mat × COV CL where COV V and COV CL denote functions of additional covariates to be searched for V and CL, respectively. Then, F ren was formulated as (Anderson et al., 2007): F ren = CLCr / CLCr TV λ C L C r = R Cr / Cr × e − k tox ∙ t R Cr = 64.2 × e k Cr ∙ a g e − 30 where for the F ren model, CLCr is creatinine clearance, CLCr TV is the typical value of CLCr corresponding to the subject with age of 30 years old, and λ is an exponent parameter, for the CLCr model, Cr is plasma creatinine concentration, R Cr is Cr production rate and k tox means a rate constant to describe the reduced creatinine clearance due to vancomycin-induced nephrotoxicity which could reduce clearance, and for the R Cr model, 64.2 (mg⁄h) is the value of R Cr for the subject with age of 30 years old, and k Cr is a rate constant related to age (Cockcroft and Gault, 1976) such that k Cr > 0 for age <30 and k Cr ≤ 0 for age ≥30.

On the other hand, F mat , which applied to patients under 4 years old, was formulated as (N. Holford et al., 2013): F mat = PMA γ / P M A 50 γ + PMA γ where PMA50 denotes PMA where maturation reaches 50% of the adult clearance and γ is a steepness factor of sigmoid function.

After integrating these into clearance, potential covariates such as age, gender, BUN, and the history of hypertension, diabetes, renal diseases, cardiovascular disease, hematological diseases, pleural effusion and edema, and sepsis were tested based on pharmacological and physiological plausibility, as denoted by COV V and COV CL . Stepwise covariate model building was conducted with a likelihood ratio test based on the criteria of p < 0.01 (ΔOFV = 6.63, df = 1) for forward selection and p < 0.001 (ΔOFV = 10.82, df = 1) for backward deletion, where OFV means objective function value of NONMEM.

To characterize the dynamics of CRP for the time delay between infection and biomarker level change in the plasma, a semi-mechanistic model with proliferation and transit compartments was attempted which was formulated as below: dProl dt = K in × 1 + S CRP ∙ D − K out × P r o l   ;   P r o l i f e r a t i o n d T r a n 1 dt = k tr × P r o l − T r a n 1   ;   T r a n s i t   1 d T r a n 2 dt = k tr × T r a n 1 − T r a n 2   ;   T r a n s i t   2 dCRP dt = k tr × T r a n 2 − k CRP × C R P   ;   C i r c u l a t i o n

Prol is the proliferation compartment, with K_in and K_out denoting production and degradation rate constant, respectively, which was assumed to be K_in = K_out = k tr for numerical simplicity, Tran1 and Tran2 are the transition compartments, with k tr denoting a transit rate constant, and CRP is the circulating compartment, with k CRP denoting an elimination rate constant. S CRP is a scaling factor and its effect is proportional to disease severity (D), which was estimated as the following equations where α is a scaling factor. dD dt = k D × D × 1 − E Drug E Drug = α × A U C

AUC obtained from the developed PK model, E Drug means drug effect and k D is a first-order progression rate constant. To find covariates for PD model parameters, a stepwise selection was done with a likelihood ratio test based on the criteria of p < 0.05 (ΔOFV = 3.84, df = 1) for forward selection and p < 0.01 for backward deletion.

To select the final model, the tested models were compared based on OFV or AIC values and the precision of parameter estimates. The selected model was then evaluated with goodness-of-fit plots such as observation versus model prediction (PRED), and conditional weighted residual (CWRES) versus PRED. Subsequently, a visual predicted check (VPC) was conducted for the validation of the selected model, ensuring that collected drug concentrations fall within the 95% confidence interval of the 2.5^th percentile, median, and 97.5^th percentile of predicted drug concentrations.

Using the final PK model, simulations were conducted to explore optimal dosage regimens for sub-populations stratified by selected covariates. To achieve this, for the m selected covariates, the range from the minimum to the maximum values for each covariate i was divided into equi-spaced intervals N i , resulting in a total of N 1 × N 2 × ∙ ∙ ∙ × N m scenarios of covariate pairs.

These scenarios were simulated using typical values for intravenous administration four times a day (QID). The objective was to determine the optimal dosage regimen that meets the criteria of achieving a trough concentration closest to the target value, which was set at 7 mg/L for children under the age of 4, 10 mg/L for those aged 4 to 19, and 15 mg/L for adults. Additionally, the peak concentrations were required to remain below the toxic level of 40 mg/L across all age groups. These criteria were based on findings from previous studies (Rybak M. et al., 2009; Rybak et al., 2020).

Using the developed PD model and the selected covariates, the time course of CRP plasma concentrations was simulated by varying the dose, infusion rate, and inter-dose interval in order to predict the optimal time to discontinue treatment and prevent unnecessary overtreatment.

A web-based application was developed using R Shiny to visualize the simulated time course of plasma concentrations of vancomycin and CRP, which is based on the developed PK-PD model and the selected dosage regimen. This application allows users to visualize the predicted individual PK-PD profile for a specific dosage regimen.

A total of 542 hospitalized patients were included in the PK model building, with 22 aged under 4, 18 aged 4 to 19, and 502 adults. Among them, 128 patients were eligible for the PD analysis. Vancomycin was administered at doses ranging from 500 to 1,500 mg, with dose intervals ranging from 6 to 24 h. The patient population spanned from neonates to the elderly, and 40 pediatric patients were included in both PK and PD analyses. The average length of hospital stay was 20 days, with a minimum of 2 days and a maximum of 113 days. The hospitalization of up to 113 days was due to a secondary infection caused by pneumonia, which necessitated the extension of vancomycin treatment. Detailed demographic information for the patients can be found in Tables 1, 2.

Various models were tested using the collected vancomycin concentration data, and a 2-compartment model incorporating allometry scaling was selected as the basic model based on the OFV value. For covariates, gender, BUN, and the history of diabetes and renal disease were selected for COV CL and age for COV V as follows: COV CL = e k BUN ∙ B U N − 15 × 1 + θ FEM ∙ F E M × 1 + θ DM ∙ D M × 1 + θ REN ∙ R E N COV V = e k V ∙ a g e − 40 where FEM = 1 for female and 0 for male, DM = 1 for diabetes and 0 for no diabetes, and REN = 1 for renal disease and 0 for no renal disease.

The incorporation of vancomycin-induced nephrotoxicity k tox ) into the model significantly improved model prediction (p < 0.0001). Vancomycin CL exhibited a gradual increase with BUN levels up to 15 mg/dL, followed by a subsequent decrease. In the final model, median values of CL and V were estimated to be 4.31 L/h and 38.6 L, respectively. The details of final parameter estimates are presented in Table 3, and these values were consistent with findings from other studies. Regarding the relative standard error (RSE), except for γ, most parameters, demonstrated reliable estimation. Goodness-of-fit plots of the PK model are shown in Figure 2, revealing the absence of noticeable trends. VPC plots in Figure 3 indicated that the majority of observations fell within the 95% confidence intervals of the predictions. In these figures, the vancomycin concentration data depicted represent observations during repeated dosing after hospitalization. The two concentration measurements at 2,500 h correspond to peak and trough samples obtained from the patient who was hospitalized for 113 days or 2,712 h.

A semi-mechanistic model with one proliferation compartment and two transit compartments was selected as the structural model. Among the covariates, the presence of pneumonia had a significant effect on the transit rate constant as in the following equation, where PNE = 1 indicates pneumonia and 0 indicates no pneumonia. k tr = θ k tr + θ PNE ∙ P N E

The final parameter estimates are presented in Table 4. Except for k D and CRP₀, the initial value of CRP, between-subject variability (BSV) could not be obtained due to numerical difficulties. Consequently, k CRP was fixed at 0.0365 h⁻¹ based on the prior knowledge that CRP’s half-life was 19 h (Vigushin et al., 1993). Mean transit times for pneumonia and non-pneumonia patients were 6.65 and 9.62 days, respectively. Figure 4 presented the goodness-of-fit plots where no obvious trends were observed, indicating, overall, the model adequately describes CRP concentrations.

Considering the range of covariates observed in the collected data, a total of 384 scenarios were simulated for adult patients, using 4 Cr levels, 4 BUN levels, 4 age levels, and 6 body weight (WT) levels. For pediatric patients, a total of 656 scenarios were simulated. Among these, 336 scenarios were designed for patients aged 4 years or younger, incorporating 7 PMA levels, 3 WT levels, 4 Cr levels, and 4 BUN levels. An additional 320 scenarios were created for patients aged over 4 years, employing 4 age levels, 5 WT levels, 4 Cr levels, and 4 BUN levels.

Daily doses were explored in increments of 0.2 g (0.05 g per dose; QID) for adults, and in increments of 5 mg/kg (1.25 mg/kg per dose; QID) for pediatrics. Detailed results are presented in Table 5 for adults, Table 6 for pediatric patients and Table 7 for those aged under 4. The MIPDs were provided for each combination of covariate values, and the shading increased with the dose. For precision dosing presented in the tables, the covariate model for COV CL indicated that dose reductions of 23%, 15%, and 20% were necessary for patients with renal diseases, diabetes and females, respectively. In simulation results, the optimal dose increased with weight, while it decreased with Cr and BUN levels. Regarding age, the dose exhibited an upward trend with age until 40, followed by a subsequent decrease. Regarding the attainment of target concentrations, it was confirmed that simulated concentrations given optimal doses were within the therapeutic range across all age groups. Figure 5 illustrates that when concentrations of the study patients were simulated using the optimal dose from a virtual patient with the most analogous covariates to the study patient, the majority of simulated concentrations lied within the therapeutic range. In contrast, a significant portion of observed concentrations fell outside this range.

The influence of dose and pneumonia on CRP concentration-time profile is illustrated in Figure 6. No dose-dependent differences were observed up to approximately 240 h or 10 days. However, after this point, CRP levels exhibited a faster decline with higher doses, returning to the normal range of 1–10 mg/L at around 1,100 h or 46 days with the recommended dose. This suggests that discontinuing vancomycin treatment at 46 days is viable to prevent unnecessary overtreatment. It is worth noting that the extended time required for CRP levels to normalize can be attributed to the prolonged hospitalization period experienced by the study patients.

The interface of the developed R Shiny tool for vancomycin dose optimization is showcased in Figure 7. This tool offers the capability to optimize the necessary dose for a patient to attain the desired drug concentration based on their individual covariates, achieved through the adjustment of the dose and/or dosing interval. By providing visualizations of drug concentration-time profiles, this tool serves as an alternative to the nomogram approach detailed in Tables 5–7. In addition to simulating drug concentrations, the application also facilitates the visualization of simulated CRP plasma concentrations that correspond to the employed dose, dosing rate, and covariates utilized in the PK simulation.