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!

Mourad HAMIMEDDear Michiel,

Thank you very much for sharing this excellent and very interesting work.

I was curious, when we have a combined (proportional + additive) residual error model, how we should set the sigma (in the actual example : sigma <- matrix(0.01)) ? Should we calculate the sigma using the following equation:

W = SQRT(IPRED**2*SIGMA(1)**2 + SIGMA(2)**2)

_____________________________________

SIGMA(1) : is the proportional error

SIGMA(2) : is the additive error

Thank you in advance

MJvanEsdonkPost authorDear Mourad,

Thank you for your message.

Please find additional information on the following github thread discussing this question:

https://github.com/metrumresearchgroup/mrgsolve/issues/149

In short,

“You can modify sig2j <- sig2j <- out[[dvcol]]^2*sig2 to be sig2j <- out[[dvcol]]^2*sig2 + sig2add"

Mourad HAMIMEDThanks again Michiel, i am grateful for your help