Applying MAP Bayes estimation for therapeutic drug monitoring (TDM) in R with mrgsolve

The development of population PK (and PD) models enable the use of individual Bayesian dose optimization. One could use the included covariates to derive the dose of an individual but one could also determine what an individual’s clearance and/or volume of distribution is by measuring one or multiple plasma concentrations after the first dose. Especially when a narrow therapeutic window is present and distribution/elimination kinetics are difficult, one could use this methodology to improve the treatment efficacy in a population.

This blog post builds further on the code and functions from the Maximum A Posteriori (MAP) Bayes Parameter Estimation example provided by mrgsolve. Currently available at the following link: https://mrgsolve.github.io/blog/map_bayes.html

In this blog post, we will compare the true pharmacokinetic profile with the estimated MAP Bayes Parameter estimates and show how the true versus the predicted multiple dose profile of an individual would look like when a sample is taken after the first dose.

First of all, MAP stands for Maximum A Posteriori estimate. In short, this is the value of a posterior distribution [random effects (e.g. IIV on clearance)] that best fits the available data of an individual. We can run an optimization algorithm across the posterior distribution for different clearances and identify which value would fit the data best, based on the distribution identified during model development. Even though this is an oversimplification of the process happening in the background of a Bayesian Therapeutic Drug Monitoring model, it is fit for purpose for this post.

For more information on obtaining the MAP estimate, have a look at the following publication: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/

---

Simulation model

The mrgsolve model code that is required for a MAP parameter estimation method is slightly different compared to a standard simulation model used in other posts. The main difference is that we need to simulate PK profiles at random, to simulate some realistic data, and simulate the PK during the optimization process. Therefore, we use both ETA(1), the normal variance parameter, and ETA1 as parameters in the model.

The following code is the population PK model saved in a separate file (Pop_PK_TDM.cpp).

$PARAM TVCL=0.5, TVVC=20, ETA1=0, ETA2=0

$CMT CENT

$PKMODEL ncmt=1

$OMEGA 0 0
$SIGMA 0

$MAIN
double CL = TVCL*exp(ETA1 + ETA(1));
double V =  TVVC*exp(ETA2 + ETA(2));

$TABLE 
capture DV = (CENT/V)*(1+EPS(1));
capture ET1 = ETA(1);
capture ET2 = ETA(2);

Compared to the mrgsolve example previously mentioned, we simulate a single dose of 750mg (intravenous) in a one compartment model, with a larger residual error component to mimic a more realistic scenario. We define all the omega’s and proportional error term and define where we take the plasma samples (0.5h, 12h, and 23.5h). Then, we run the mrgsolve simulation.

##################################################################################
############################ CODING BLOCK ########################################
##################################################################################
# Run-Time Environment: R version 3.6.1
# Author: MJ van Esdonk
# Short title: Pharmacokinetic TDM
# Version: V.1.0
##################################################################################
##################################################################################
###
rm(list = ls(all.names = TRUE))
###

###############################################
############ Load libraries
library(mrgsolve)
library(dplyr)
library(ggplot2)
library(minqa)
library(nloptr)


######################################################################
### Load in model for mrgsolve
mod <- mread_cache("Pop_PK_TDM")

### Set omega and sigma of the model
omega <- cmat(0.23,
              -0.78, 0.62)

sigma <- matrix(0.01)

dose <- ev(amt=750,cmt=1)

## Set the sampling times
sampl <- c(0.5,12,23.5)

# Set seed for random simulation
set.seed(1012) 

## Simulate mrgsolve model at selected timepoints
sim <- 
  mod %>%
  ev(dose) %>%
  omat(omega) %>%
  smat(sigma) %>%
  carry_out(amt,evid,cmt) %>%
  mrgsim(end=-1, add=sampl) %>% # only simulate at selected timepoints
  as.data.frame()

Once we simulated the PK profile of 1 individual and selected the time points at which we have our blood samples, we need to create a vector of the true simulated ETA1 and ETA2. We can use these later for comparison of the true PK profile of an individual compared with the MAP based prediction. As ETA1 and ETA2 are unknown when clinical data is obtained, we remove these from the dataset (ET1 and ET2 in the output object) during the optimization process.

## Create vector with the true estimates for both ETA's
true_ETA1_2 <- sim %>% slice(1) %>% select(ET1,ET2)


# Set DV to NA for the dosing record and round DV with 2 digits
data_for_optim <- sim %>%
  mutate(DV = ifelse(evid==1, NA_real_, round(DV,2)))%>% 
  select(-ET1, -ET2)

Now that we have obtained the PK concentrations of an individual, it is time to introduce the MAP function and optimize ETA 1 and 2. Have a look at the mrgsolve blog post for an introduction to this function. As good coding practice is that all parameters used in a function should also be forwarded as arguments to that particular function, the example MAP function has been slightly adapted and extensively annotated compared with the mrgsolve blog post.

### Annotated mapbayes function
mapbayes <- function(eta, # Eta vector to optimize
                     omega_matrix, # Matrix of omegas
                     sigmasq, # Matrix of sigma squares
                     init, # Initial eta's to retrieve the ETA labels
                     d, # Dataframe with observations and doses
                     ycol, # Column name with observations from d
                     m, # mrgsolve model (or other type of ODE solver model)
                     dvcol, # Column with DV from model 
                     pred=FALSE # Use parameters for predictions (TRUE) or optimization (FALSE)
                     ) {
  
  ## Make sure sigma squared is numeric
  sig2 <- as.numeric(sigmasq)
  
  ## solve omega matrix
  omega.inv <- solve(omega_matrix)
  
  ################
  # Make a list of the eta's to optimize
  eta <- as.list(eta)
  
  # Get ETA names from the init 
  names(eta) <- names(init)
  
  eta_m <- eta %>% 
    unlist %>% 
    matrix(nrow=1)
  
  ################
  # Update model parameters (ETA1 and ETA2) with current set of eta values
  m <-  param(m,eta)
  
  
  ################
  ## Run a simulation, with the new ETA1 and ETA2 values but without any other level of variability
  out <- m %>% 
    zero_re() %>% 
    mrgsim(data=d,output="df")
  
  ################
  ## If prediction is asked, return simulated values
  if(pred) return(out)
  
  
  ################
  ###### Objective function value to minimize
  # http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/
  sig2j <- out[[dvcol]]^2*sig2
  sqwres <- log(sig2j) + (1/sig2j)*(d[[ycol]] - out[[dvcol]])^2
  nOn <- diag(eta_m %*% omega.inv %*% t(eta_m))
  
  return(sum(sqwres,na.rm=TRUE) + nOn)
}

Then it is time to get back our estimates based on the simulated data. We define some random starting values for both ETA’s that we want to optimize and we start the newuoa function. This function changes the values defined in init in order to retrieve the lowest value (hence, best model fit). We expand the call to this function with all the different argument required.

### Initial values of the ETA's to optimize
init <- c(ETA1=-0.3, ETA2=0.2)

######## Optimizer function based on the data_for_optim
fit <- nloptr::newuoa(init,
                      mapbayes,
                      omega_matrix=omega,
                      sigmasq=sigma,
                      init=init,
                      d=data_for_optim,
                      ycol="DV",
                      m=mod,
                      dvcol="DV",
                      pred=F)

Running this function gives a fit object that we can use to print the optimized parameters. Congratulations, we have performed a MAP Bayesian optimization and estimated the individual clearance and distribution volume based on 3 plasma samples!

### Print the MAP Bayes estimates for the available dataset
fit$par

The beauty of this function is, that once the parameters have been optimized, we can use the same function for simulations! We only have to change the argument pred to TRUE and make sure that our model is simulating for a longer interval and not only for the selected timepoints. We save this in the object pmod (prediction mod) which is identical to the mod with an updated simulation end time and simulated delta.

##################################################
### Predict the PK profile
# Select dosing rows only
pdata <- data_for_optim %>% filter(evid==1)

# Set prediction model parameters
pmod <- mod %>% update(end=24, delta=0.1) 

# Simulate
pred_optimized <- mapbayes(fit$par,    
                           omega_matrix=omega,
                 sigmasq=sigma,
                 init=init,
                 d=pdata,
                 ycol="DV",
                 m=pmod,
                 dvcol="DV",
                 pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "MAP optimized individual PK")

We do this three times, once we use the fit$par as input to simulate the optimized parameters, once we use true_ETA1_2 to simulate the true PK profile of this individual, and once we use an empty vector to simulate the population profile (all ETA’s are 0). The three simulations are combined and plotted to judge the performance.

# Simulate true PK profile of this individuals
pred_true <- mapbayes(true_ETA1_2, # Replace with true ETA 1 and 2
                      omega_matrix=omega,
                           sigmasq=sigma,
                           init=init,
                           d=pdata,
                           ycol="DV",
                           m=pmod,
                           dvcol="DV",
                           pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "True individual PK")

# Simulate true PK profile of this individuals
pred_pop <- mapbayes(c(0,0), # Replace with population model (eta = 0)  
                     omega_matrix=omega,
                      sigmasq=sigma,
                      init=init,
                      d=pdata,
                      ycol="DV",
                      m=pmod,
                      dvcol="DV",
                      pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "Population PK profile")


## Combine dataframes together
pred <- pred_optimized %>%
  bind_rows(pred_true) %>%
  bind_rows(pred_pop)


######## Plot the pharmacokinetic profile 
plot_pk <- ggplot(pred, aes(x=time,y=DV,col=Simulation)) +
  ## Add line for each simulate profile
  geom_line()+

  # Set axis and theme
  scale_x_continuous("Time (h)",breaks=seq(0,24,4))+
  scale_y_continuous("Concentration (mg/L)")+
  
  theme_bw(base_size=16)+
  labs(col='Simulated profile')+
  
  ## Add observations to plot
  geom_point(data=sim %>% filter(evid==0), aes(time,DV), col="darkslateblue",size=3) 
  

## Save plot
ggsave(plot=plot_pk,filename='PMX_Solutions_TDM.png',width=9)

As can be observed from this figure, the true PK profile and the profile with the MAP optimized parameters are very similar, and deviations from the PK profile are mainly driven by the residual variability in the model in this scenario. One can play around with the parameters, the model structure, and sampling times to judge the performance.

Multiple dose simulation

Once this code is ready and the MAP estimates do indeed describe the observations, we use the estimated parameters for the simulation of some multiple doses. We only need to adapt the dosing dataframe and the simulation time and interval to simulate the PK over multiple doses. This simulation provides very useful information and can, for example, be used to show the improvement in the use of TDM compared with model-derived dosing regimens, which do not account for inter-individual variability on the parameters.

##################################################
################### Simulate multiple doses
### Modify dosing row
pdata_MD <- pdata %>% mutate(addl=5,ii=24) ## 5 additional doses with a 24h interval

# Set prediction model parameters
pmod <- mod %>% update(end=144, delta=1) 

# Simulate multiple dose
pred_optimized_MD <- mapbayes(fit$par,       
                              omega_matrix=omega,
                           sigmasq=sigma,
                           init=init,
                           d=pdata_MD,
                           ycol="DV",
                           m=pmod,
                           dvcol="DV",
                           pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "MAP optimized individual PK")

# Simulate true PK profile of this individuals multiple dose
pred_true_MD <- mapbayes(true_ETA1_2, # Replace with true ETA 1 and 2      
                         omega_matrix=omega,
                      sigmasq=sigma,
                      init=init,
                      d=pdata_MD,
                      ycol="DV",
                      m=pmod,
                      dvcol="DV",
                      pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "True individual PK")

# Simulate true PK profile of this individuals multiple dose
pred_pop_MD <- mapbayes(c(0,0), # Replace with population model (eta = 0)  
                        omega_matrix=omega,
                     sigmasq=sigma,
                     init=init,
                     d=pdata_MD,
                     ycol="DV",
                     m=pmod,
                     dvcol="DV",
                     pred=TRUE) %>% 
  filter(time > 0) %>%
  mutate(Simulation = "Population PK profile")


## Combine dataframes together
pred_MD <- pred_optimized_MD %>%
  bind_rows(pred_true_MD) %>%
  bind_rows(pred_pop_MD)



######## Plot the pharmacokinetic profile 
plot_pk <- ggplot(pred_MD, aes(x=time,y=DV,col=Simulation)) +
  ## Add line for each simulate profile
  geom_line()+
  
  # Set axis and theme
  scale_x_continuous("Time (h)",breaks=seq(0,144,24))+
  scale_y_continuous("Concentration (mg/L)")+
  
  theme_bw(base_size=16)+
  labs(col='Simulated profile')+
  
  ## Add observations to plot
  geom_point(data=sim %>% filter(evid==0), aes(time,DV), col="darkslateblue",size=3) 

## Save plot
ggsave(plot=plot_pk,filename='PMX_Solutions_TDM_MD.png',width=9)

Dosing adjustments

Based on the presented methodology, the simulated profiles of an individual can be analyzed and judged whether a certain target maximal concentration or area under the curve has been reached. A new dose can be given to this subject and the simulation repeated to determine whether an improved dose might be available.

COMMENT

Any suggestions or typo’s? Leave a comment or contact me at info@pmxsolutions.com!

---

Continue reading:

• 
  Exploring the Capabilities of ChatGPT: A Step-by-Step Guide to Creating a Pharmacokinetic Analysis Shiny App
• 
  Using R Shiny applications in scientific research – Can you spot a drug effect on blinded pharmacodynamic data…?
• 
  Introduction to pharmacometrics and its importance in drug development