Introduction

Therapeutic drug monitoring (TDM) concerns the measurement of drug concentrations in patients to optimize dosing schedules in individual patients and is commonly used in treatment optimization of patients of intensive care unit (ICU). Population pharmacokinetic (PK) models are regularly used to derive optimized dosing regimens based on TDM data. These model-based dosing regimens are guided by Bayesian forecasting through maximum a posteriori (MAP) estimation. Using MAP estimation, individual PK parameters are estimated using a previously developed population PK model based on collected historical drug concentration TDM data of a patient who has received the drug. These individual PK parameters in turn can be used to perform Bayesian forecasting to predict the prospective concentrations to further derive a dosing schedule that meets therapeutic concentration targets associated with efficacy or toxicity.

There is currently a lack of consensus or guidelines regarding the use of historically collected TDM PK data for Bayesian forecasting. Given that TDM is typically associated with sparse PK samples, inclusion of all available TDM data may arguably lead to optimal use of data and as a consequence individual PK parameter estimates closer to the true value in a patient. However, in ICU patients, rapid alterations in organ functions associated with PK might also actually lead to increased bias in the current PK parameter estimates if too much historical TDM data is included (1). In ICU patients it is therefore unclear if and how such historical TDM data should be included for Bayesian forecasting, and, if current MAP-based approaches optimally make use of available historical data in estimating individual PK parameters. In the current analysis we aimed to address these questions using a representative TDM data set for vancomycin in ICU patients. Vancomycin forms a cornerstone antibiotic agent for treatment of sepsis associated with gram-positive infections in the ICU. In the ICU population, vancomycin usually shows large inter-individual variability as well as large intra-individual variability in PK (2,3). TDM is routinely performed for vancomycin in the ICU due to its narrow therapeutic window and rapid alterations in organ function that could lead to changes in PK (3,4,5,6). As such, vancomycin represents a relevant paradigm drug to study the optimization of model-based TDM approaches.

In this analysis we propose and evaluate two methods of MAP estimation for model-based dose optimization using TDM data. We will refer to the newly proposed methods as adaptive MAP and weighted MAP. The proposed methods are compared to standard MAP estimation.

Materials and Methods

Data

A retrospective vancomycin TDM data set collected in ICU patients was used for this analysis, containing vancomycin concentration-time data, vancomycin dosing histories, and patient demographics, collected from two hospitals in The Netherlands. The patients from one hospital received vancomycin with a loading dose of between 1000 mg and 2000 mg followed by doses of 1000 mg twice a day until doses were adjusted at the discretion of the treating clinician. The patients from the other hospital received vancomycin with a loading dose, followed by continuous infusion. Both loading doses and continuous infusion doses were individualized based on the advice from the in-house-developed model-based TDM software AutoKinetics (7). The infusion duration time ranged from one to two hours for non-continuous infusion and was 24 h for continuous infusion. Blood samples were obtained on average twice or three times every week mainly at trough level, i.e. prior to the next dose. A serum creatinine measurement was also measured for each sample. The creatinine clearance of the patients were calculated by default based on the MDRD formula, provided that the commonly seen formulae for calculating glomerular filtration rate were expected to perform similarly in ICU patients (8,9). We excluded the patients if plasma concentration samples were only available for a single day of treatment since model predictions cannot be validated for such patients. As a result, a data set consisting of 2435 of concentration data points from 408 patients was used (Table I).

Table I Characteristics of the Patients in this Study
Full size table

Population Pharmacokinetic Model

A previously published one-compartmental population PK model of vancomycin by Roberts et al. 2011 (Table II) was used to perform both the MAP estimation and Bayesian forecasting (10). The model has been validated in our own ICU population, with the same data as used in this study (11).

Table II The Vancomycin PopPK Model (10)
Full size table

MAP Estimation

Data splitting was needed in order to evaluate different MAP methods. The data set was split into segments where each segment contained the concentration data of a single day. Thus, one segment corresponded to one day which contained one trough sample. For each patient, we estimated individual PK parameters based on one to up to seven consecutive segments of data to predict the concentration data of the succeeding segment. Hence, a segment of data was regarded as either historical data if used for MAP estimation, or as prospective data if used to validate the model predictions. We compared three estimation methods: standard MAP, adaptive MAP and weighted MAP. Figure 1 depicts the conceptual differences between these methods, and are defined further in below.

Fig. 1
figure 1

Schematic diagram of methods of MAP estimation used in the study. Standard MAP executes estimation once using all historical TDM data (a); Adaptive MAP executes estimation iteratively using each segment of historical TDM data with updated prior mean by posterior mode from its preceding iteration and repeats until the last iteration, i.e. 0 to m1, …, mn-1 to mn (b). Weighted MAP executes estimation once using all historical TDM data with weighted importance (likelihood) of each segment of data during the estimation (c).

Full size image

Standard MAP

Given a model with random-effect (of clearance (CL) and volume of distribution (V)) parameter vectors η, the standard MAP estimation constructs a posterior density distribution as follows (12):

$$ z\left(\boldsymbol{\upeta} \right)=f\left(\boldsymbol{\upeta} |\mathbf{Y}\right)g\left(\boldsymbol{\upeta} \right) $$
(1)

where f(η|Y) is the likelihood of η given concentration observation vector of Y. g(η) denotes the prior density distribution of η, which followed a multi-normal distribution informed by the vancomycin population PK model (Table II). z(η) is the posterior density distribution that MAP estimation aims to maximize so that the mode (maximum) value of z(η) is the final estimate of η.

Adaptive MAP

The adaptive MAP shared a similar idea of model predictive control theory (13). The estimation was executed iteratively, in which each iteration used one segment of the historical data. For the 1st iteration, we executed MAP estimation using the earliest one segment of data to obtain the random-effect parameters (posterior mode) which were denoted as m1. For the 2nd iteration, the means of the variance of the random-effect parameters (prior mean), which are 0 under usual assumptions, were substituted by m1. Then the MAP estimation was executed using the next segment of data. By this way, the prior means of each iteration (except for the first one) were updated by the posterior modes of the last iteration. Such process was repeated until only one segment of data was left, which was used to validate model predictions (Fig. 1b).

Weighted MAP

Like the standard MAP estimation, the weighted MAP estimation was executed once using all historical data (Fig. 1c). However, the contribution of each segment of historical data to the total likelihood was weighted according to the following functions:

$$ z\left(\boldsymbol{\upeta} \right)=f{\left(\boldsymbol{\upeta} |\mathbf{Y}\right)}^{w\left({\mathbf{Y}}_s\right)}g\left(\boldsymbol{\upeta} \right) $$
(2)
$$ w\left({\mathbf{Y}}_s\right)={\left(\frac{\Delta {T}_{Ref}}{{\Delta T}_{{\mathbf{Y}}_s}}\right)}^{\alpha } $$
(3)

where w(Ys) is a weighting function for the likelihood of the parameter to be estimated given observation Ys, which is a subset of Y e.g. the concentration data at a particular day prior to the data to be predicted. \( {\Delta T}_{{\mathbf{Y}}_s} \) stands for the time distance in days between the historical data Ys and the data to be predicted. ΔTRef and α are both weighting factors. ΔTRef is the reference day, defined as the cutoff value of the time distance where w(Ys) is 1, i.e. where ΔTRef equals \( {\Delta T}_{{\mathbf{Y}}_s} \). α is the unitless effect size of the weighting function w(Ys) on the likelihood f(η|Y). For the sake of simplicity, we only tested a number of integer combinations of ΔTRef and α including integers from 1 to 7 for ΔTRef, and integers 1 to 5 for α. By this method, the individual likelihood of each data point to the total likelihood was exponentiated by w(Ys) so that its importance for the MAP estimation was weighted.

Evaluation of MAP Methods

The evaluation of the MAP methods was based on the use of different amounts of historical TDM data for each patient to estimate the PK parameter and predict a prospective drug concentration. To this end, a segment of data may be used as historical data for MAP estimation or as prospective data to validate model predictions in different scenarios. The prediction of prospective data was done using historical data from preceding one up to seven days (Fig. 2).

Fig. 2
figure 2

Schematic diagram of including historical data for MAP estimation and Bayesian forecasting of one patient.

Full size image

To evaluate predictive performance of each method, we calculated the difference between the observed and predicted vancomycin concentration, by calculating the percentage error (PE) as follows:

$$ {PE}_i=\frac{1}{n_i}\sum \limits_{j=1}^{n_i}\frac{{\hat{Y}}_{ij}-{Y}_{ij}}{Y_{ij}}\times 100\% $$
(4)

Here, PEi denotes percentage error for ith subject; ni equals to the number of forecasted concentration of ith subject; \( {\hat{Y}}_{ij} \) and Yij represent the predicted and observed value of the jth prospective concentration of ith subject, respectively. We normalized the PE to the number of data points of the patient, to equalize the weight of each patient in the results. The accuracy and the precision were measured as the median and the interquartile range (IQR) of the PE, respectively. The results of PE were also visualized using a bar plot for standard MAP and adaptive MAP methods, and a heat map for weighted MAP method.

Implementation

We executed the MAP estimation with all three methods using historical data to obtain individual PK parameters. We also estimated PK parameters using each segment of data set to explore the trend of the parameters’ change over time. Nonlinear mixed-effects modeling software (NONMEM, version 7.4.4; ICON Development Solutions, MD, USA) was used for both MAP estimation and Bayesian forecasting. The NONMEM code of the implementation of all methods is available as supplementary material. Data organization and visualization were carried out with R (version 3.6.0; R-project.org).

Results

Impact of the Methods of MAP Estimation

The median PEs of including preceding one up to seven days of historical data for the standard MAP method ranged from −3.4% to −11.5% and for the adaptive MAP method, −3.4% to −5.4% (Fig. 3a, b). For the weighted MAP method, the median PEs ranged from −2.3% to −8.5% (using optimal values for weighting factors ΔTRef=4, α=2 based on the mean and median values of median PEs that were closest to 0) (Fig. 3c). The adaptive MAP method thus outperformed the standard MAP method, whilst the effect of weighted MAP method was limited. We did not find a clear effect of the weighting factors ΔTRef and α (Fig. 4). The IQR of the PE for the standard MAP method and the weighted MAP method were both visibly wider than that for the adaptive MAP method, indicating the adaptive MAP method was most precise. An alternative error bar plot based on mean values and 95% confidence intervals is available in the supplementary material.

Fig. 3
figure 3

The percentage error of Bayesian forecasting using the standard MAP method (a), the adaptative MAP method (b), and the weighted MAP method using optimal weighting factors ΔTRef=4 and α=2 (c). The squares are the median values, and the bars represent 25% and 75% quantile lines respectively. The labeled text is the number of patients that were included for the calculation.

Full size image
Fig. 4
figure 4

The percentage error of Bayesian forecasting using the weighted MAP method for all combinations of weighting factors ΔTRef and α. ΔTRef and α are both weighting factors. ΔTRef is the reference day, defined as the cutoff value of the time distance where w(Ys) is 1, i.e. where ΔTRef equals \( {\Delta T}_{{\mathbf{Y}}_s} \). α is the unitless effect size of the weighting function w(Ys) on the likelihood f(η| Y).

Full size image

Impact of the Number of Days of Included Historical Data

Including only one preceding day of historical data led to the most accurate Bayesian forecasting (Figs. 3 and 4). The PE worsened noticeably when including more historical data and using the standard MAP method and the weighted MAP method, while was consistent when using the adaptive MAP method (Figs. 3 and 4). The adaptive MAP method was therefore most robust to the included amount of historical data. The IQR of PE did not show clear relationship with the number of days of included historical data, regardless of the MAP method. However, the adaptive method resulted in the narrowest IQR which decreased markedly when using more than one preceding day of historical data (Fig. 3b). Thus, including the preceding two days of historical data for the adaptive MAP method guaranteed both accuracy and precision. Given the current setting of data splitting, including previous e.g. three segments of historical data can be also seen as including one segment of historical data which contains three days of concentration data. Hence, the similar results would be expected.

Time Course of the Random-Effects Parameters

The random-effects of both PK parameters, CL and V, were estimated for each day (Fig. 5). The random-effect of CL showed a decreasing trend over time while that of V was rather stable. Such a result is in agreement with the observation that the PE values are all negative, since using a retrospective CL of a higher value would result in underprediction of the prospective data.

Fig. 5
figure 5

Time course of random-effects on PK parameters over time. CL, clearance; V, volume of distribution. The squares are the mean values, and the bars represent mean plus standard deviation and mean minus standard deviation lines respectively.

Full size image

Discussion

Adequate dosing is crucial for successful treatment of critically ill patients. Bayesian forecasting is probably the most commonly used approach in model-based TDM for ICU patients. There is unfortunately no evidence how Bayesian forecasting should be made properly. Our study investigated to optimize the use of historical TDM data in this special population.

Although potentially counter-intuitive, we found that including more historical TDM data can result in reduced forecasting accuracy. From the clinical perspective, this is not surprising since PK properties of ICU patients can change rapidly with disease progression. This means that older historical data contains outdated information that leads to bias in Bayesian forecasting. Thus, TDM samples that are drawn more closely to the ones to be forecasted in terms of time, are more likely to have better predictive abilities. Similar findings were reported in a previous study that if only a single sample was utilized, Bayesian predicted concentrations were less accurate when obtained using the first (‘oldest’) observed concentration compared with the most recent observed vancomycin concentration (14). By contrast, in the general patient population, it was reported that taking more TDM data into account did improve the performance of Bayesian forecasting for vancomycin (15).

The median PE values across scenario’s (Figs. 3 and 4) are all negative indicating concentrations are typically underpredicted. We indeed observe a decreasing trend of the random-effects on the CL over time (Fig. 5). Considering that the time-varying covariate creatinine clearance, which is the only potential time-varying covariate identified in most published vancomycin PK models (11,16), was taken into account, a likely explanation is that there could still be other time-varying factors that are difficult to be captured by covariates, e.g. health condition. Patients staying long at the ICU usually require more clinical care indicating a deteriorated health condition and possibly worsening organ function. Such results demonstrated that changes in health condition and related PK in ICU patients may occur and could limit the utility of previously published models for model-based TDM. Our analysis provided a possible solution to best deal with this situation. Clinical practitioners should be well aware of the disease progress of the patient and should primarily take the most recent history into account for Bayesian forecasting. Further improvement of the model is desirable to account time-varying PK, but is not the purpose of this study.

The adaptive MAP method turned out to be most accurate and precise for Bayesian forecasting. A similar approach has been recently proposed for a tacrolimus PK model and it showed that such an approach can be more robust with regard to model misspecification (17). Compared to the standard MAP method, the adaptive MAP method led to a smaller median PE with a smaller accompanying IQR except for including only the preceding one day of historical data. Only one (the first) MAP iteration needs to be executed in such a case, so that the adaptive MAP is identical to the standard MAP (Fig. 3a, b). When including historical data from further past, the adaptive MAP method extracts the information of a segment of historical data in one iteration and carries it forward into the next iteration. Hence, the prior distribution is always updated by the last posterior which is tailored by the historical data from the preceding iteration. After the last iteration, the accumulated information of the historical data was optimally “resembled” through the iterative procedure rather than averaged at once as is done in the standard MAP method. Consequently, the prediction is generated with smarter use of data and thus the accuracy and the precision are improved. Although the most accurate Bayesian forecasting was observed when only including the preceding one day of historical data for the adaptive MAP method, including historical TDM data from the preceding two days is favorable as precision improved drastically with negligible loss in accuracy, which may cause clinically relevant consequences. Provided an actual AUC0–24 is 480 mg∙h/L, the IQR of its prediction would be from 382 to 537 mg∙h/L using the preceding one day of TDM data (IQR of PE from −15.2% to 11.9%), while that would be from 443 to 517 mg∙h/L using the preceding two days of TDM data (IQR of PE from −7.7% to 7.7%). As an AUC0–24between 400 and 600 mg∙h/L is now recommended the target range for effective and safe vancomycin treatment, the former would unrightfully indicate a dose adjustment while the latter would not (3).

Regarding the weighted MAP method, the underlying hypothesis was that with increasing “age” of the historical data their carried information on the PK of the drug would be of less importance for predicting the data in the future. Mathematically, it equates with weighting the likelihood contribution of each data point to the total likelihood. The time distance between the historical data and the prospective data was chosen to construct the weighting function given our hypothesis. The main reason we additionally added two parameters ΔTRef and α to the weighting function was to incorporate flexibility in the exploration on how much effect time distance could have on the likelihood as well as on the Bayesian forecasting. Unfortunately, the weighted MAP method showed limited improvement in Bayesian forecasting compared to the standard MAP method. The reason is probably that the weighted MAP method alters the likelihood by relocating the relative importance of data points based on the time distance but does not affect the prior distribution. Given that TDM data are usually sparse, the prior distribution is likely to have a greater impact on the MAP estimation than the data does. Hence, the weighted MAP method is not able to optimize the MAP estimation as much as the adaptive estimation does. We evaluated the values of ΔTRef and α using only a limited number of integer combinations, but this should not have a harmful impact on the results. For ΔTRef, there is an obvious interpretability that it refers to the time distance. For α, its value is not anticipated to be considerably large or small, which indicates to either diminish the contribution of prior distribution or the likelihood to the MAP estimation. There might be some precision loss not evaluating non-integers, but the result should be covered in the investigated range and the precision may not be much off.

Despite that this study focuses on vancomycin only, the findings are potentially generalizable for other TDM drugs as well. First, the ICU data set used in this study is large. This ensures the ICU patient population is well represented which guarantees the validity of finding. Second, the evidence of the findings in this study is not drug specific. It is ICU patient population rather than vancomycin that results in the observed altering PK properties in the data. It is thus physiologically plausible to assume that similar behaviors may be observed for other drugs in ICU patient. Nonetheless, the applicability of the finding to other TDM drugs in ICU patients needs to be verified and cautions should be taken when making decisions in clinical practice.

Since Bayesian forecasting in routine practice is usually based on the MAP estimation, it inherits the properties as well as the drawbacks of the point estimation. As the PK models used for Bayesian forecasting are usually nonlinear models, the point estimate of MAP does not necessarily relate to the most probable clinical outcome. MAP estimation does not address the uncertainty of the posterior distribution and hence is not able to quantify the potential risk of the relevant clinical outcome. Recently, MAP-based approaches have been questioned in certain circumstances, such as when trough concentration is of interest because even a subtle discrepancy could be damaging for concentrations at a low level (18). Although a full Bayesian approach is superior to MAP-based approaches in handling mentioned issues, the benefit from such an approach could be limited or perhaps clinically irrelevant, due to the relatively large noise in the routine TDM data from ICU patients which influences the results as well. Considering the requirement of extensive domain knowledge in non-clinical fields, a full Bayesian approach is conceivably difficult to be implemented for routine practice and has poor translational ability to communicate to the clinical professionals.

Nonparametric estimation has been proposed and has shown to outperform parametric estimation in the cases where the assumption of normal distribution of the data is violated, while the nonparametric estimation does not make any assumptions as such (19). Indeed, nonparametric approach by theory could be better suited for a particular type of tasks. Yet, it has not been widely adopted in clinical practice and thus the experience is limited in the community. Further investigation and understanding of the differences between parametric and nonparametric approaches and their applicability are certainly required but are out of the scope of this study. Although the parametric MAP-based approaches may have limitations, we believe they are likely to retain the mainstay across the community for a considerable time.

There are limitations worth mentioning regarding this study. First, the creatinine clearance of patients was calculated using the MDRD equation normalized to body surface area. We were not able to calculate the absolute MDRD (which we should ideally have as the original population PK model used measured creatinine clearance) since height (required to estimate body surface area) was not available at the time the data set was obtained. This may affect the descriptive ability of the model for our data set. However, considering that the used model was externally validated using the same data set (11) and MAP estimation can correct individual PK parameters, the consequence of using normalized MDRD should be modest. Second, since vancomycin concentration is usually measured sparsely in routine clinical practice, the impact of the intensity of samples within one dosing interval was not investigated with TDM data in this study. This however may be of interest for other TDM drugs for which samples can be measured multiple times in a dosing interval. Third, we identified that the most recent historical data have best predictive value. Yet, we didn’t investigate what the time distance to the last historical data could be to still have acceptable predictive performance of Bayesian forecasting. This was because the sampling scheme of vancomycin followed the local policy, where trough samples are collected prior to the next dose which are usually scheduled regularly. The time distance to the last historical data is thus usually fixed. Yet, it is our expectation that a greater time distance to the last historical data comes with higher PE, as was also shown by Broeker et al. in a mixed patient population receiving vancomycin (14). Last, although the time-varying covariate creatinine clearance was taken into account in the analysis, further inclusion of time-varying characteristics is expected to improve the performance of Bayesian forecasting using the standard MAP method. This can however be challenging due to the need for sufficient data abundancy over the course of therapy. Besides, the inclusion of inter-occasion variability can probably improve the model and also predictive performance of the standard MAP method, but it is often difficult to define occasions in the local data set in the same way as was done for the development of the population PK model.

In conclusion, the adaptive MAP method seems to be promising for model-based TDM of vancomycin in ICU patients with better performance in Bayesian forecasting than the standard MAP method. One should include historical TDM data of the preceding one day when using the standard MAP method or the weighted MAP method, or of the preceding two days when using the adaptive MAP method for the most reliable Bayesian forecasting of vancomycin for ICU patients.