FormalPara Key Points

Current dosing recommendations fail to achieve consistent vancomycin exposure throughout life and across patient subgroups.

This pharmacokinetic model for vancomycin adequately characterises vancomycin pharmacokinetics across a broad range of patient populations, including those at the extremes of age and weight.

1 Introduction

The glycopeptide antibiotic vancomycin plays an important role in the treatment of Gram-positive bacterial infections. It is currently considered a key therapeutic option in the context of the treatment of serious infections (e.g. community- and hospital-acquired pneumonia, complicated skin and soft-tissue infections, infective endocarditis) and for the empiric treatment of patients at high risk of infections caused by methicillin-resistant Staphylococcus aureus or multidrug-resistant Staphylococcus epidermidis [1, 2].

Although vancomycin has been in clinical use since its first registration almost six decades ago, questions with respect to the optimal dosing regimen remain. Consequently, a plethora of reports on subpopulation- and context-specific pharmacokinetics and derived dosing regimens have surfaced over the last 25 years. In 2012, these numerous population-pharmacokinetic (PK) studies led Marsot and co-workers [3] to conclude that initial dosing regimens for vancomycin should depend on creatinine clearance and body weight in adults, and creatinine clearance, body weight and age in children. The assessment report on vancomycin drug products issued in 2017 by the European Medicines Agency [1] is in line with the recommendations by Marsot and co-workers. Nevertheless, it also highlights some known unknowns. For example, the report comments that specific dosing regimens might be required for neonatal and elderly populations and for patients with renal impairment.

Besides getting the initial dose right, optimisation of vancomycin exposure entails the use of therapeutic drug monitoring (TDM), especially in populations where PK variability is large (e.g. preterm neonates, critically ill patients) [1]. Individualised dosing based on TDM data is most efficiently achieved through Bayesian forecasting of PK parameters based on an a priori population-PK model. This technology is currently offered by a range of TDM software packages [4], such as DoseMe® (DoseMe Pty Ltd., Brisbane, QLD, Australia), InsightRx® (Insight Rx Inc., San Francisco, CA, USA), MwPharm ++® (Mediware a.s., Prague, Czech Republic) and DosOpt (University of Tartu, Tartu, Estonia) [5].

At present, different vancomycin population-PK models are, however, used across these software packages; and within a package often different models are used for different patient subpopulations. This poses a challenge to clinicians who have to consider the limitations of the models they utilise. Moreover, the clinician might have to switch models when treating different patient populations.

It has previously been shown that when applied to external datasets, there are wide differences in the performance of various currently used PK models [6, 7]. In anaesthesia, this has led to data-sharing initiatives with the intent of replacing different subgroup-specific models with a single pooled population-PK model [8, 9]. Such a model is expected to be more generalisable than other models, making it simpler and more useful for everyday clinical use [10].

In this study, we aimed to develop a single PK model for vancomycin, derived from a broad dataset covering the extremes of patient characteristics. Moreover, as a benchmark for current dosing recommendations, we performed simulations using the final population-PK model and evaluated the expected vancomycin exposure throughout life and for specific patient subgroups. Finally, based on these simulations optimised doses were recommended.

2 Materials and Methods

2.1 Component Datasets

Studies on vancomycin pharmacokinetics were identified though a PubMed search (until 27 September, 2017) using search terms: “vancomycin AND population AND pharmacokinetics [Title/Abstract]” OR “vancomycin AND pharmacokinetics [Title/Abstract]”. We excluded studies of patients receiving continuous renal replacement therapy, haemodialysis, haemodiafiltration, extra-corporeal membrane oxygenation and non-intravenous administration of vancomycin. Corresponding or senior authors from these studies were invited to contribute their anonymised data and participate in this modelling study. All studies obtained necessary institutional review board approval, as declared in the original papers or as declared by the corresponding or senior author (personal communication).

Postmenstrual age (PMA), postnatal age (PNA), weight, height and serum creatinine (SCR) were extracted from the datasets. Postmenstrual age for patients other than neonates was assumed to be 40 weeks longer than the recorded postnatal age (years). Missing values for SCR were imputed as the population median SCR from that study. Missing values for height were imputed as the population median height from that study or from a similar study (for neonatal studies with missing height) or with height data from a national health survey in that specific population.

Questionable covariate data were corrected whenever possible. These primarily concern overlapping drug infusion records and unrealistic combinations of age, weight and height within a subject. To avoid computational difficulties during model building due to excessively long follow-up times, the dataset was restricted to observations from the first 31 days of therapy.

2.2 Population-Pharmacokinetic Modelling

The vancomycin concentration vs. time data were fitted using the FOCE-I estimation algorithm in NONMEM® (Version 7.3; GloboMax LLC, Hanover, MD, USA). The “tidyverse” package (Version 1.1.1.; Wickham H. 2017) in R® (R Foundation for Statistical Computing, Vienna, Austria) was used to graphically assess the goodness of fit and for simulations.

As a starting point, one-, two- and three-compartmental PK models were fitted to the data. Inter-individual variability on the typical population parameter estimates was assumed to be log-normally distributed. Residual unexplained variability was modelled using a combined proportional and additive error model. Inter-occasion variability was not tested in the model owing to difficulties in defining an occasion in the context of this dataset where dosing and sampling times/frequency vary considerably.

Modifications to the structural and/or covariate model were accepted only if they resulted in a decrease in the objective function value (OFV) and an increased in-sample predictive performance. A decrease in OFV was judged statistically significant if inclusion of an additional parameter decreased the OFV by more than 3.84 points. As a measure of in-sample prediction performance, the absolute relative prediction errors (APEs) were calculated (Eq. 1) by comparing the measured vancomycin concentrations for each individual i at time point j (VANij) with the population-predicted concentrations (PRED). The median of the distribution of APEs (MdAPE) was used to compare the imprecision of candidate models during model building and to compare the final model against earlier published models.

$${\text{APE}} \left( \% \right) = \left| {\left. {\frac{{{\text{VAN}}_{ij} - {\text{PRED}}_{ij} }}{{{\text{PRED}}_{ij} }}} \right|} \right. \times 100\% .$$
(1)

To ascertain that model building was not driven by the most populated subgroup, we stratified the calculation of the prediction error. The following subgroups were created: pre-term and term newborns (PMA < 0.87 years), children and adolescents (age < 18 years), adults (age < 65 years), elderly (age < 80 years), very elderly (age ≥ 80 years), underweight adults (age > 18 years and body mass index < 18.5 kg m−2) and obese adults (age > 18 years and BMI > 30 kg m−2). During model building, the average MdAPE across the aforementioned subgroups was taken as an overall measure of predictive performance. Covariates tested for inclusion in the model were: weight (kilograms), PMA (years), PNA (days), SCR (mg dL−1), sex, critical illness and presence of severe burn injuries. Height and related covariates such as body mass index and fat-free mass were not tested as height was missing in most neonates and some adults.

2.3 A Priori Included Covariate Models

In line with earlier work by Holford et al. [11] and Germovsek et al. [12, 13] prior to inclusion of additional covariates in the model we corrected for size and maturational changes. For this, PK parameters were scaled to weight according to allometric theory [14], with an exponent of 1 for volume terms (V1, V2, V3) and an exponent of 0.75 for clearance terms (CL, Q2, Q3). A sigmoidal maturation function was used to scale clearance with PMA according to Eq. 2 with γ defining the steepness of the non-linear relationship and PMA50 being the PMA when clearance reaches 50% of maximal values.

$${\text{Maturation function}} = \frac{{{\text{PMA}}^{\gamma } }}{{{\text{PMA}}^{\gamma } + {\text{PMA}}_{50}^{\gamma } }}.$$
(2)

2.4 Testing Serum Creatinine as a Covariate in the Model

The influence of SCR on clearance was tested using an exponential function as shown in Eq. 3. This function predicts decreasing clearance with increasing SCR values with θSCR defining the steepness of the relationship. Serum creatinine was included in the dataset as a time-varying covariate with backward constant interpolation between observations.

$$F_{\text{Scr}} = e^{{\left( { - \theta_{\text{SCR}} \times ({\text{SCR}} - {\text{SCR}}_{\text{std}} )} \right)}} .$$
(3)

Parameter values were standardised by centering SCR observations according to a standardised SCR (SCRstd) value. Different approaches for defining SCRstd were compared. A first approach consisted of using the median observed SCR from the dataset. Other approaches were explored in an effort to derive age-, weight- and sex-adjusted SCRstd values. For this, empiric functions were fitted to the covariate data using fractional polynomial regression (R package “mfp: Multivariable Fractional Polynomials”; Version 1.5.2) in line with the work of Ceriotti et al. [15], or custom non-linear models in R®.

2.5 Evaluation and Optimisation of Current Dosing Recommendations

As examples of contemporary dosing guidelines for vancomycin approved by the European Medicines Agency and the US Food and Drug Administration (FDA), we used the summary of product characteristics (SmPC) for “Vancomycin 500 mg Powder for Solution for Infusion” (Consilient Health Ltd; available from www.medicines.org.uk; consulted on 23 May, 2018) and the label for “Vancomycin Hydrochloride for Injection USP” (ANI Pharmaceuticals, Inc.; available from www.fda.gov; consulted on 23 May, 2018). The posologies outlined by the SmPC and the FDA label are summarised in Tables S1 and S2 of the Electronic Supplementary Material (ESM), respectively.

Steady-state vancomycin exposure was simulated using the post hoc-predicted clearance for all patients in our dataset according to Eq. 4:

$${\text{AUC}}_{{24\;{\text{h}}}} = \frac{{{\text{DOSE}} \times \frac{24}{\tau }}}{\text{CL}},$$
(4)

where AUC24 h is the area under the concentration–time curve for a 24-h period in steady state and DOSE and τ are the amount of vancomycin administered and the dosing interval, respectively. The predicted AUC24 h was then used to calculate the proportion of patients who were under- and over-dosed. Under-dosing was assumed when the AUC24h was below 400 mg L−1 h, a frequently cited threshold for efficacy [16, 17] for the treatment of pathogens with a minimum inhibitory concentration of 1 mg L−1. Over-dosing was defined as an AUC24h exceeding 700 mg L−1 h, a recently advocated safety threshold for renal toxicity [18].

In the final part of this study, we attempted to optimise the daily doses in the SmPC and the FDA label. For this, in line with the above, post hoc-predicted clearance for all patients was used to predict the proportion of patients achieving an AUC24h above 400 mg L−1 h and below 700 mg L−1 h (fTarget). A non-linear optimisation routine based on the Nelder-Mead Simplex [19] as implemented in R® (R package “nloptr”; Version 1.0.4) was used to optimise fTarget using the doses in the SmPC and the label as control variables. For practical purposes, optimised daily doses were rounded to the nearest 50 mg and 0.5 mg kg−1.

3 Results

3.1 Data

In total, 39 publications were identified. All senior and/or corresponding authors for these publications were contacted via e-mail on multiple occasions. We obtained 14 previously published data sets [20,21,22,23,24,25,26,27,28,29,30,31,32,33]. Our pooled dataset contains information across a broad range of patient subgroups ranging from premature neonates [20, 24, 30, 32] to extremely obese patients [26] and from adult Japanese healthy volunteers [27] to critically ill and trauma patients [25, 31, 33]. A summary of the included studies is shown in Table 1. The distributions of patient characteristics and vancomycin concentrations are shown in Fig. S1 of the ESM. In total, 8300 vancomycin drug concentrations from 2554 individuals were included in the dataset. Individuals were adults (720), newborns (559), elderly patients (512), obese adults (274), very elderly patients (213), underweight adults, (158) and children and adolescents (118).

Table 1 Details of the component datasets. The distribution of patient covariates is summarised by the median and the range; the percentage of missing values is shown in square brackets

Initial analysis of the dataset identified outliers (n = 10). The data for these individuals were primarily biased because of missing dosing information prior to the first observation, potential sampling errors (e.g. reported vancomycin concentrations > 100 mg L−1) or dosing regimens that were significantly different from other comparable individuals from the same study. Because of the small number of outliers identified (0.4% of individuals), we decided to remove these patient records from the dataset.

3.2 Population-Pharmacokinetic Modelling

We found that a two-compartment model better described the data compared with a one-compartment model (ΔOFV = − 1155). A three-compartment model had a lower OFV (ΔOFV = − 85). However, goodness-of-fit plots for both models were indistinguishable and likelihood profiling [34] revealed high uncertainty on the estimate for the inter-compartmental clearance to the slow peripheral compartment (Q3). Therefore, we decided to retain the two-compartment model as the final model.

We found significant maturation in vancomycin clearance, with 50% of maximal clearance being reached by 46.4 weeks PMA (PMA50). After correcting for maturation, it became apparent that vancomycin clearance deteriorates with ageing. The sigmoidal function fitted to describe the declining clearance with age revealed that 50% of maximal clearance is lost by 61.6 years PMA (AGE50). Figure 1 shows the typical-for-PMA standardised clearance (L h−1 70 kg−1) for all patients in our dataset. Figure 1 also shows the maturation–decline function as estimated by Lonsdale et al. [35] in a recent pooled analysis of three commonly used beta-lactam antibiotics.

Fig. 1
figure 1

Standardised clearance [CLstd] (L h−1 70 kg−1) for vancomycin throughout life. Typical vancomycin clearance according to our final model is shown with a solid black line. Post-hoc CLstd values for all patients in our dataset are shown with solid grey circles. The grey shaded area denotes the region between the 10 and 90% percentile of all observations. The dashed black line is the maturation–decline function for beta-lactam antibiotics according to Lonsdale et al. [35] PMA postmenstrual age

After accounting for size- and age-related changes in vancomycin pharmacokinetics, we tested the influence of SCR on clearance according to Eq. 3. Table 2 provides an overview of the different model building steps and shows the different approaches taken to define SCRstd. When the median observed SCR (0.80 mg dL−1 or 70.7 µmol L−1) is used as SCRstd, the lowest OFV is obtained (ΔOFV = − 1810). However, this approach led to an increase in the estimates for PMA50 and AGE50. A fractional polynomial model relating SCRstd to PNA, gestational age, weight and sex performed less well (ΔOFV = − 1692). However, standardisation of SCR values using this approach had no influence on the estimates for PMA50 and AGE50. Ultimately, we found that the empiric function producing the lowest OFV with the least influence on the estimates for PMA50 and AGE50 was a non-linear function relating SCRstd to PMA (Eq. 5).

Table 2 Model building hierarchy
$${\text{SCR}}_{\text{std}} = e^{{\left( { - 1.228 + \log_{10} \left( {{\text{PMA}}\left( {\text{yr}} \right)} \right) \times 0.672 + 6.27 \times e^{{\left( { - 3.11 \times {\text{PMA}}({\text{yr}}} \right)}} } \right)}} .$$
(5)

A plot of the predicted SCRstd as a function of PMA is shown in Fig. S2 of the ESM. Figure S2 also shows the SCRstd used in the work of Johansson et al. [36] and Hennig et al. [37], which is based on the reference intervals reported by Ceriotti et al. [15] and Junge et al. [38] When we implemented the latter approach, we found a similar decrease in OFV compared to our best performing function (ΔOFV = − 1790 vs. −1805, respectively). However, as seen in Table 2, under this approach the estimates for PMA50 and AGE50 slightly increase (+1.8 weeks and +4.1 years, respectively). Furthermore, as seen in Fig. S2 of the ESM, the predicted SCRstd does not align well with the median of our observed SCR values. Therefore, we decided to include the empiric function described in Eq. 5 in our final model. The inclusion of SCR decreased in-sample MdAPEs for all subgroups ranging from − 0.7% for children and adolescents to − 10.2% in elderly patients. Inclusion of SCR as a baseline covariate rather than a time-varying covariate was significantly worse (ΔOFV = + 876).

At this step in model development, we tested compartmental allometry for clearance and Q2, scaling them allometrically to the estimated size of the corresponding compartment (V1 and V2). In line with earlier work on propofol [8], remifentanil [9] and dexmedetomidine [39], compartmental allometry improved the goodness of fit (ΔOFV = − 212) without disturbing the interpretation of the fixed-effect parameters in the models.

Altered pharmacokinetics was observed in two composite datasets. On the one hand, in the study by Buelga and co-workers on patients with haematological malignancies [29], clearance was significantly higher (+29.4%). On the other hand, the study by Lo and co-workers [32], who used heel-prick sampling as opposed to arterial/venous sampling, we found a lower V1 (− 31.2%) and Q2 (− 59.7%). Accounting for these subgroup-specific effects resulted in a better fit (ΔOFV = − 184 and − 115, respectively).

To achieve adequate numeric stability of the model, i.e. smooth gradient minimisation and acceptable likelihood profiles, we removed the estimate for the population variability on Q2 (ΔOFV = + 5.0).

The final population-PK model is shown in Eqs. 613 and Table 3, the NONMEM output for the final model is shown in the ESM.

Table 3 Parameter estimates and associated relative standard errors (RSEs) for the final population-pharmacokinetic model. Inter-individual variability (IIV) associated with the typical parameters is expressed as coefficient of variation %
$${\text{CL}} ({\text{L}}\;{\text{h}}^{ - 1} ) = \theta_{\text{CL}} \times \left( {\frac{V1}{{\theta_{V1} }}} \right)^{0.75} \times F_{\text{Mat}} \times F_{\text{Decline}} \times F_{\text{SCR}} \times \left( {1 + \theta_{{{\text{STDY}}10}} } \right) \times e^{{\eta_{1} }}$$
(6)
$$V1 (L) = \theta_{V1} \times \left( {F_{\text{Size}} } \right)^{1} \times \left( {1 - \theta_{{{\text{STDY13}}\_V1}} } \right) \times e^{{\eta_{2} }}$$
(7)
$$V2 (L) = \theta_{V2} \times \left( {F_{\text{Size}} } \right)^{1} \times e^{{\eta_{3} }}$$
(8)
$$Q2 ({\text{L h}}^{ - 1} ) = \theta_{Q2} \times \left( {\frac{V2}{{\theta_{V2} }}} \right)^{0.75} \times \left( {1 - \theta_{{{\text{STDY}}13\_Q2}} } \right)$$
(9)
$$F_{\text{Size}} = \frac{{{\text{WGT}} ({\text{kg}})}}{70}$$
(10)
$$F_{\text{Mat}} = \frac{{{\text{PMA}} ({\text{wk}})^{{\gamma_{1} }} }}{{{\text{PMA}} ({\text{wk}})^{{\gamma_{1} }} + {\text{PMA}}_{50}^{{\gamma_{1} }} }}$$
(11)
$$F_{\text{Decline}} = \frac{{{\text{PMA}} ({\text{yr}})^{{ - \gamma_{2} }} }}{{{\text{PMA }}({\text{yr}})^{{ - \gamma_{2} }} + {\text{AGE}}_{50}^{{\gamma_{2} }} }}$$
(12)
$$F_{\text{SCR}} = e^{{\left( { - \theta_{\text{SCR}} \times \left( {SCR \left( {{\text{mg}} {\text{dL}}^{ - 1} } \right) - {\text{SCR}}_{\text{std}} } \right)} \right)}} ,$$
(13)

In these equations, CL (L h−1) is vancomycin clearance, V1 (L) and V2 (L) are the central and peripheral volume of distribution and Q2 (L h−1) is the inter-compartmental clearance. FSize, FMat, FDecline and FSCR describe size-related changes, maturational changes, age-induced deterioration and SCR-related changes in vancomycin pharmacokinetics, respectively. θSTDY10, θSTDY13_V1 and θSTDY13_Q2 denote the increased clearance in patients with haematological malignancies and the effect of heel-prick sampling on V1 and Q2. Symbols η1η3 denote inter-individual variability (with variances ω1ω3) of the typical PK parameters.

Backwards elimination of FSize, FMat, FDecline or FSCR from the model consistently resulted in an increased OFV and increased unknown inter-individual variability in clearance (data not shown). Goodness-of-fit plots for the final model are shown in Fig. 2. Figures S3, S4 and S5 of the ESM show the goodness-of-fit plots stratified by study and the prediction-variance-corrected visual predictive check [40] for the pooled data and data stratified by patient category, respectively. Overall, these diagnostics show that our final model is adequately developed.

Fig. 2
figure 2

Goodness-of-fit plots for the final population-pharmacokinetic model. Scatterplots show the distribution of the observed vancomycin concentrations (Observed Cplasma) vs. population and individual predictions and the (absolute) conditionally weighted residuals (|CWRES| and CWRES) vs. individual predictions and time after the end of the dose. Negative times denote observations taken when drug was infused, whereas positive times are observations after stopping the infusion. A dashed line denotes the line of unity or the zero line, whilst a red solid line shows a non-parametric smoother

3.3 Evaluation and Optimisation of Current Dosing Recommendations

The distribution of predicted steady-state vancomycin AUC24h resulting from the dosing recommendations in the SmPC and the FDA label are shown in Fig. 3. The proportion of patients achieving an AUC24h below 400 mg L−1 h (fAUC<400), between 400 and 700 mg L−1 h (fTarget), and above 700 mg L−1 h (fAUC>700) are shown in Tables S3 and S4 of the ESM.

Fig. 3
figure 3

Simulated area under the concentration–time curve at steady state for a 24-h period (AUC24h) according to the different patient groups in the summary of product characteristics (SmPC) and the US Food and Drug Administration (FDA) label. White boxplots show the AUC24h distributions resulting from the original SmPC and FDA label whereas grey boxplots show the distribution of AUC24h for the optimised doses for the SmPC and the label. The grey shaded area denotes the target exposure, i.e. an AUC24h between 400 mg L−1 h and 700 mg L−1 h. EMA European Medicines Agency

Based on the dosing regimen in the SmPC, 19% of patients attain an AUC24h below 400 mg L−1 h, whilst 46% of patients attain an AUC24h above 700 mg L−1 h. Elderly (76%), very elderly (70%) and obese patients (81%) are mainly at risk for attaining a steady-state AUC24h above 700 mg L−1 h, whereas for underweight adults (25%), newborns (43%), and children and adolescents (65%) a considerable proportion of patients is expected to attain an AUC24h below 400 mg L−1 h.

The dosing regimens in the FDA-approved label result in a more consistent steady-state exposure with a very low risk of attaining an AUC24h above 700 mg L−1 h across patient subgroups. However, except for newborns, all patient subgroups are significantly at risk for a steady-state AUC24h below 400 mg L−1 h, with probabilities ranging from 52% in very elderly patients to 70% in children and adolescents.

The optimised daily doses for the SmPC and FDA label are shown in Table 4. In general, optimised daily doses are higher. Exceptions are adults and children with an estimated glomerular filtration rate < 50 mL.min−1, where the optimised daily doses for the SmPC are 14.5 and 28.5 mg kg−1 instead of 15 and 30 mg kg−1. For both the optimised SmPC and the FDA label, fTarget increases consistently across subgroups with the overall fTarget increasing from 35 to 46% and 37 to 60% for the SmPC and the FDA label, respectively. For the optimised label, fTarget is > 50% in all subgroups and the risk for over- or under-dosing is consistent (± 20%) across subgroups. However, for the optimised SmPC, children and adolescents (30%) and underweight adults (46%) remain at risk for under-dosing, whereas elderly patients (46%), very elderly patients (53%) and obese patients (58%) are at risk for over-dosing.

Table 4 Posology from the US Food and Drug Administration (FDA) label for “Vancomycin Hydrochloride for Injection USP” from the FDA website (www.fda.gov; consulted on 23 May, 2018) and the summary of product characteristics for “Vancomycin 500 mg Powder for Solution for Infusion” available from www.medicines.org.uk (consulted on 23 May, 2018)

4 Discussion

In this population-PK modelling study, we showed that vancomycin pharmacokinetics undergoes significant changes throughout human life. Besides size-related changes in clearance, maturation is the main driver for clearance during early childhood with 90% of adult clearance being reached by 2 years PMA. Between 3 and 16 years PMA, size-corrected clearance is highest (> 5.05 L h−1 70 kg−1). Afterwards, clearance deteriorates significantly, reaching 2.66 h−1 70 kg−1 (i.e. 50% of maximum clearance) by 61.6 years PMA. Not surprisingly, we found that SCR is a significant covariate on clearance for both between- and within-subject variability. When SCR decreases or increases 0.20 mg dL−1 (17.7 µmol L−1) from the typical-for-PMA SCRstd, clearance changes with + 13.8 and − 12.2%, respectively.

According to our model, typical vancomycin clearance in a 60-year-old, 65-kg patient with a SCR of 0.97 mg dL−1 (85.7 µmol.L−1) is 2.55 L h−1 (0.039 L h−1 kg−1). Vancomycin clearance in a 32-week, 1.5-kg neonate with a SCR of 0.64 mg dL−1 (56.6 µmol L−1) is 0.0756 L h−1 (0.0504 L.h−1.kg−1). Both values are in line with the range of clearance estimates reported by Marsot et al. 3 for adults (0.031–0.086 L h−1 kg−1) and paediatric patients (0.020–0.112 L h−1 kg−1). In contrast to Marsot et al. [3] who reported different estimates for vancomycin volume of distribution for adults (0.864 L kg−1) and children (0.565 L kg−1), we found that the size-standardised volume of the distribution is constant throughout life at 0.61 L.kg−1 for V1 and 0.59 L kg−1 for V2.

Our estimate for PMA50 of 46.4 weeks is in good agreement with earlier work by Rhodin et al. [41] who found that glomerular filtration function reaches 50% of adult values by 47.7 weeks PMA. This suggests, in line with previous reports [42], that vancomycin elimination depends predominantly on glomerular filtration. Moreover, our PMA50 estimate is not significantly different from the findings of Anderson et al. [20] who reported a PMA50 of 33.3 weeks (95% confidence interval 14.8–51.8) in a cohort of preterm neonates. It also inspires confidence that our conclusion of maturation being nearly complete at 2 years PMA is supported by a recent physiologically based PK simulation study by Calvier et al. [43].

Model building started out with an a priori model consisting of allometric scaling and a maturation function as advocated by Holford et al. [11] and more recently by Germovsek et al. [12]. By extending this standardised model with an additional sigmoidal decline function, we were able to describe the deterioration in vancomycin clearance with age. Recently, Lonsdale et al. [35] have used the same methodology to describe beta-lactam antimicrobial pharmacokinetics from early life to old age. Notwithstanding that tubular secretion and/or re-absorption are likely involved in beta-lactam elimination, these authors found a PMA50 estimate of 49.7 weeks, which is in line with our findings. Their estimated AGE50 of 86.8 years, however, is significantly higher than what we found for vancomycin (61.6 years). In addition to differences in elimination processes, the different handling of renal function in the study by Lonsdale et al. [35] compared with our analysis might explain the different AGE50 estimates. Indeed, as seen from Table 2, during model building, AGE50 ranged from 67.4 to 79.8 years depending on the method used to standardise SCR.

This is by far the largest population-PK modelling study on vancomycin pharmacokinetics. On the one hand, our findings confirm earlier fragmented findings on the influence of weight, kidney function and age on vancomycin pharmacokinetics by other groups [3] as well as the increased clearance in patients with haematological diseases [29]. On the other hand, we were able to negate some prior hypotheses. We found that after correcting for age, body size and kidney function, PK parameters were not different between patients with or without burns [23], between critically ill or non-critically ill patients [25], or between patients treated with intermittent or continuous dosing [44]. In addition, we found no differences in PK parameters with respect to the different vancomycin assay methods used (i.e. assays based on fluorescence polarisation or turbidimetric inhibition principles) [45]. Supportive evidence for these findings is provided in Figs. S6–S9 of the ESM.

Patients from two studies differed significantly from the rest of the population, thereby forcing us to include two study-specific parameters. First, patients in the Buelga et al. study [29] had a significantly higher clearance. This finding is in line with earlier work by Jarkowski et al. [46] and Zhao et al. [47] who found higher clearance in adult and paediatric patients with malignant haematological diseases. Second, the concentration–time profile for patients in the Lo et al. study [32] was different from other studies. This is likely owing to the heel prick sampling that was used in the Lo et al. study as opposed to venous or arterial blood sampling in the other studies. This reasoning is in line with earlier work by Chiou who showed that drug concentrations depend on the sampling site [48].

Using the post hoc PK parameters derived from our model, we evaluated two currently used dosing regimens for vancomycin and showed that these dosing regimens do not consistently result in efficacious and safe steady-state vancomycin exposure across patient populations. The optimised daily doses increase the proportion of patients attaining an AUC24h between 400 and 700 mg L−1 h for both the European Medicines Agency-approved SmPC and the FDA-approved label. From these simulations, we see that weight-based dosing for neonates is acceptable, yet the dose-per-kilogram should be increased. For adults, weight-based dosing is inappropriate and results in under-dosing of underweight adult patients whilst obese adult patients are generally over-dosed. Moreover, age- or kidney function-adjusted dosing is necessary to avoid over-dosing in elderly patients and very elderly patients. The optimised FDA label, which defines kidney function-adjusted doses for patients aged older than 18 years and weight-based dosing otherwise, performs well with consistent target attainment and minimal (± 20%) risk of under- and over-dosing across patient subgroups.

These recommendations should be interpreted by taking into account that: (1) further optimisation of these dosing regimens could be achieved by further stratifying patient subgroups based on age and/or weight; (2) the thresholds used for defining adequate and or excessive AUC24h are still under debate and regional differences exist in minimum inhibitory concentration distributions; and (3) TDM is still necessary to further optimise the treatment of individual patients. Furthermore, our study excluded specific patient categories (patients receiving renal replacement therapy, haemodialysis, haemodiafiltration and extra-corporeal membrane oxygenation) and in the pooled dataset some patient categories were only sparsely populated (e.g. the group of children and adolescents). As such, the generalisability of our results and the predictive performance of the final model should be validated in a prospective study prior to implementation of the model or derived recommendations in clinical practice.

5 Conclusion

In this article, we show that vancomycin pharmacokinetics changes dramatically throughout human life and a single population-PK model is able to capture these changes. Through simulations, we showed that current dosing regimens do not succeed in providing similar steady-state exposure across patient subgroups and as such the probability for therapeutic success or vancomycin-induced toxicity is not constant throughout life. As illustrated in this work, this new model has the potential to overcome these limitations and could be used to further develop age- and renal function-specific dosing regimens for vancomycin. Furthermore, this model could facilitate individualised vancomycin dosing through Bayesian forecasting based on therapeutic drug monitoring.