Model building was based on the data from two datasets. The first dataset (patient dataset) contained retrospective PK data and patient characteristics extracted from electronic health records of St. Antonius Hospital in the Netherlands. Patients were included in this analysis when they (1) were overweight or obese (BMI ≥ 25 kg/m2), (2) received vancomycin as a continuous infusion or intermittent dosing and had at least one vancomycin plasma concentration measurement within 168 h after the last dose, and (3) had regular serum creatinine measurements across vancomycin treatment. Patients could be admitted to the ward or ICU and were excluded if they received renal replacement therapy during vancomycin treatment. Vancomycin was administered according to protocols for routine clinical practice, including TDM-based initial dosage and dose adjustments. The second dataset (healthy obese dataset) was obtained from the AMIGO trial, a prospective PK study that has been published previously [7]. Data from 20 (morbidly) obese (BMI ≥ 35 kg/m2) subjects undergoing bariatric surgery were included in the analysis. All individuals had normal renal function defined as an estimated glomerular filtration rate (eGFR) ≥ 60 mL/min/1.73 m2. Subjects from this study were all dosed with 12.5 mg/kg vancomycin (maximum 2500 mg) by single intravenous infusion. Rich sampling was performed, yielding a median number of 11 samples per individual within 24 h after infusion.

In both datasets, vancomycin concentrations were measured with a validated and commercially available immune assay (VANC3) using Cobas platforms (Roche Diagnostics GmbH, Mannheim, Germany). The lower limit of quantification was 4 mg/L and the upper limit of quantification was 80 mg/L. The coefficient of variation was ≤ 4.4.

For the first TDM dataset that was obtained during routine clinical care in the hospital, the Institutional Review Board (IRB) approved the protocol and waived the requirement of informed consent. The second dataset was part of the AMIGO trial, which was approved by the local human research and ethics committee (registered in the Dutch Trial Registry [NTR6058]) and all participants provided informed consent. All studies were conducted in accordance with the Declaration of Helsinki and the International Conference on Harmonization Guidelines in Good Clinical Practice.

A population PK analysis was performed using NONMEM 7.4 (ICON Development Solutions, Hanover, MD, USA). Pearl-speaks-NONMEM (PSN) 4.9.0 [11], Pirana 2.9.7 (Certara USA, Inc, Princeton, USA) [12], and R 4.0.3 (Xpose4 package 4.7.1) were used to organize, evaluate, and visualize data. ADVAN11 TRANS4 subroutine with the first-order conditional estimation method with interaction was used.

The structural model was selected from one-, two-, and three-compartment candidate models. Inter-individual variability on the parameters was assumed to follow a log-normal distribution and tested for statistical significance on each parameter. Inter-occasion variability (IOV) was investigated. Additive, proportional, and combined error models were tested to describe residual unexplained variability. For model selection, a drop in objective function value (OFV, − 2 log-likelihood) of more than 3.84 was considered statistically significant (p < 0.05, assuming a χ2 distribution).

Covariates that were considered in the analysis were age, sex, bodyweight-related descriptors [including total bodyweight (TBW), lean bodyweight (LBW), adjusted bodyweight (ABW), ideal bodyweight (IBW), and BMI], and renal function estimates (quantified using Cockcroft-Gault calculated with LBW [CG-LBW] or TBW [CG-TBW], modification of diet in renal disease [MDRD], or chronic kidney disease epidemiology collaboration [CKD‐EPI]). MDRD and CKD‐EPI de-indexed for body surface area (BSA) were also tested. The equations for calculating bodyweight descriptors and renal function can be found in Supplementary Table S1. ICU admission status, indicating whether a patient was admitted to ICU, was also explored as a covariate. Potential covariates were implemented into the model separately using a linear function or a power function with an estimated exponent. A covariate relationship was considered statistically significant when OFV dropped more than 10.83 (p < 0.001). The model fits were assessed using goodness-of-fit (GOF) plots stratified for the BMI group [overweight (BMI 25–30 kg/m2) and obese (BMI > 30 kg/m2)], stratified for the estimated renal function group (CKD-EPI < 30, 30–50, 50–90, > 90 mL/min/1.73 m2) and ICU admission status. Plots of individual influence on change in OFV (also known as shark plots) in Xpose 4 were used to ascertain that statistical significance was driven by a sufficient number of individuals and it was checked that the inclusion of a covariate reduced the random inter-individual variability on a parameter [13]. In addition, covariate relationships were visually inspected in plots of eta value versus individual covariate value for time-independent covariates, conditional weighted residuals (CWRES) versus individual covariate value for time-dependent covariates, and individual post hoc parameter estimates versus covariate value overlain with lines representing the obtained covariate relationships, with the stratifications previously mentioned. Finally, relative standard errors of estimates (RSE) were considered to indicate adequate precision if they were below 50%.

In addition to the diagnostic plots used during model development, the robustness of the parameter estimates of the final model was validated by a bootstrap resampling analysis using PsN [12], with 1000 replicates stratified on the two datasets.

The predictive performance of the final model was assessed and visualized by a normalized prediction distribution error (NPDE) analysis using the NPDE package in R [14]. The analysis was based on 1000 simulations of the dataset.

On the basis of the obtained covariate relationships in the final model and clinical feasibility, AUC-targeted dosing recommendations for continuous infusion and intermittent dosing regimen were developed, which were assessed in stochastic simulations using the final vancomycin PK model. For representative virtual individuals with assigned values of indexed CKD-EPI (20, 40, 60, 80, and 100 mL/min/1.73 m2), bodyweight (80, 120, 160, and 200 kg), and ICU admission status (non-ICU versus ICU), 1000 simulations with inter-individual variability were performed. The target exposure for efficacy was defined as AUC24/MIC ≥ 400 mg·h/L. An AUC24/MIC threshold of 600 mg·h/L was selected as the safety limit to minimize the probability of nephrotoxicity according to the literature [15,16,17]. Since the MIC was assumed to be 1 mg/L in this study, AUC24/MIC is abbreviated as AUC24 for simplicity. Therefore, optimal dose regimens should achieve the target AUC24 from the first day of dosing onward. To explore the concentrations that are obtained upon these AUC24 400–600 mg·h/L targeted dosing recommendations, concentration–time profiles for the different bodyweight and renal function groups were also explored.