Creating a simple pharmacometric Shiny application with mrgsolve in R – Part 1

This two part series will show you how to create a simple pharmacometric Shiny application using the mrgsolve package in R. I think that Shiny applications are the most promising way to communicate our complicated mathematical models to others and I would therefore like to show how to create such an application from scratch, including all the codes used with an explanation.

There have been already some publications regarding Shiny applications for pharmacometrics online which I can strongly advice to take a look at:

• https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4394611/
• http://webpopix.org/shiny/ShinyExamples.html

In order to create a Shiny app, I commonly apply this two step approach:

1. Create the R code needed to run a single simulation
2. Modify this script to implement it in a Shiny app

By dividing the development process in two steps, we can quickly examine our script and identify possible bugs. This first part will show the development of a relatively simple R code in mrgsolve to perform PK simulations. If you already have a working R script that you would like to convert to a Shiny app, feel free to skip to Part 2!

---

mrgsolve

mrgsolve is my go-to ordinary differential equation (ODE) solver due to its speed, ease of use, and similarities with NONMEM.

The homepage of mrgsolve, with many resources, can be found here and is definitely a page to add to your bookmarks!: 

• http://mrgsolve.github.io

What I noticed in the last year of working with this package is that there are many different ways to create your dataset to be used for simulation. Therefore, the way I suggest here is not the only solution and there are (probably) shorter codes or more efficient ways possible to perform your simulations with mrgsolve. Send me an e-mail if you have any suggestions!

Before we embark on this project, lets define what kind of simulations we want our shiny app to perform in the end:

1. Pharmacokinetic simulations of a simple 2-compartment model with oral or I.V. administration
2. Manually specify the parameter estimates (theta’s, omega’s and sigma’s)
3. Generate a figure of the simulation
4. Have the option to save the simulated data to a .csv file

Defining a list with the requirements at the start helps us to keep track of what functionality our Shiny app needs to have and avoids spending time on meaningless functions.

mrgsolve model code

As said previously, there is a resemblance between mrgsolve and NONMEM. Mainly in the way the model code is written as you can see below.

This is an example code of a two-compartment model that I commonly use for my simulations. Of course you can also have a look at the pre-specified models by using the following command in R:

modlib(list=TRUE)

Two-compartment mrgsolve model code (popPK.cpp)

$PARAM
TVKA=1, TVCL = 1, TVVC=1, TVVP1=1, TVQ1=0

$CMT GUT CENT P1

$MAIN
double KA = TVKA*exp(ETA(1));
double CL = TVCL*exp(ETA(2));
double VC = TVVC*exp(ETA(3));
double VP1 = TVVP1*exp(ETA(4));
double Q1 = TVQ1*exp(ETA(5));

double k10 = CL/VC;

double k12 = Q1/VC;
double k21 = Q1/VP1;

$OMEGA 0 0 0 0 0

$SIGMA @labels PROP ADD
0 0

$ODE
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT -k12*CENT+k21*P1 - k10*CENT;
dxdt_P1 = k12*CENT-k21*P1;

$TABLE
double IPRED = CENT/VC;
double DV = IPRED*(1+PROP)+ADD;

while(DV < 0) {
  simeps();
  DV = IPRED*(1+PROP)+ADD;
}

$CAPTURE IPRED DV

So let's go through the model step by step. 

The $PARAM section at the top gives the parameters used in the model. These are now given the values 1, 1, 1, 1 and 0 respectively but we will change this in our simulation (remember, we needed to be able to manually change these estimates).

The $CMT defines the names of the compartments and the $MAIN defines the THETA's, ETA's and the calculation of the rate constants (can you spot the similarities with a NONMEM code).

The $OMEGA and $SIGMA initialize the different sources of variability (inter-individual and residual) at 0. This will also be updated by the script.

The $ODE are the 3 differential equations needed to simulate a 2-cmt model with a depot. 

The $TABLE and $CAPTURE is the place in the model where we are defining what variables we want to obtain from our simulation. We can calculate the individual prediction (without residual variability) by dividing the amount in the central compartment by its respective volume of distribution (CENT/VC). If you are simulating a PD response, dependent on a simulated concentration, this can also be implemented at this location. Then, we can simulate the DV (with residual variability) by adding the additive and/or proportional sigma's as: 
double DV = IPRED*(1+PROP)+ADD;
However, if we have an additive residual error, it can result in negative concentrations. This is of course not physiologically possible (but it is mathematically) and this is where the while loops comes into play. This loop is being executed when there are negative simulated observations and it keeps sampling from the residual distribution until a positive observation is obtained.

The R code

I have a personal preference of defining all the objects that can/may change in a script at the top. In our scenario, this would be the
dosing information, parameter estimates, and the simulation characteristics.

We define the oral or I.V. dosing as a boolean (true or false) and the other objects with their corresponding values.

##################################################################################
############################ CODING BLOCK ########################################
##################################################################################
# Run-Time Environment:               R version 3.5.1
# Author:			      M.J. van Esdonk
# Short title:			      Pharmacokinetic simulations in R 
# Version:			      V.1.0
##################################################################################
##################################################################################
###
rm(list = ls(all.names = TRUE))
###

version <- 1.0

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

###############################################
### Dosing

oral <- F
dose <- 100 # mg

###############################################
## Parameter estimates 
ka  <- 1 # Absorption rate constant
cl  <- 1 # Clearance
vd  <- 1 # Central distribution volume
vd2 <- 1 # Peripheral distribution volume
q1  <- 1 # Inter-compartmental clearance (set to 0 for 1-cmt model)


###############################################
## Insert etas and sigmas

etaka  <- 0.09
etacl  <- 0.09
etavd  <- 0.09
etavd2 <- 0.09
etaq1  <- 0.09

#################################################
## Insert residual error

sigmaprop <- 0.01 # Proportional error
sigmaadd  <- 0 # Additive error

###############################3
## Simulation info

nsamples <- 1000 ### Number of simulated individuals
sim_time <- 12 ## Time of simulation

# set probabilities of ribbon in figure
minprobs=0.1
maxprobs=0.9

Then, we need to create the dosing events. This can either be oral (dosing in cmt 1) or I.V. (dosing in cmt 2), based on a simple if statement. We also need to fill in the pre-specified dose. We assume here that there is a bolus injection but this can always be altered by specifying the rate. The same columns as in NONMEM can be used to specify different dosing regimens with ADDL, ii, rate, etc. You can also rbind different administration objects to generate more complex dosing regimens.

###################################################
############# Set Dosing objects

if(oral){
  ## Oral dose
  Administration <-  as.data.frame(ev(ID=1:nsamples,ii=24, cmt=1, addl=9999, amt=dose, rate = 0,time=0)) 
}else{
  ## IV BOLUS
  Administration <-  as.data.frame(ev(ID=1:nsamples,ii=24, cmt=2, addl=9999, amt=dose, rate = 0,time=0)) 
}

## Sort by ID and time
data <- Administration[order(Administration$ID,Administration$time),]

We can then load the model code (saved in the working directory as popPK.cpp), define the omega and sigma matrix and add the user specified parameters as columns to the dosing dataset. Adding the columns with the same names as in $PARAM will overwrite these values which enables us to easily update our model code instead of changing the .cpp file every time.

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

## Specify the omegas and sigmas matrices
omega <- cmat(etaka,
                0, etacl,
                0,0,etavd,
                0,0,0,etavd2,
                0,0,0,0,etaq1)
  
  
sigma <- cmat(sigmaprop,
                0,sigmaadd)
  
## Set parameters in dataset
    data$TVKA <- ka
    data$TVCL <- cl
    data$TVVC <- vd
    data$TVVP1 <- vd2
    data$TVQ1 <- q1

---

mrgsolve simulation

We now have all the information needed to perform a simulation. The model code, the dosing dataset with the parameters, and the two matrices with the variability. This can now all be implemented in a single piping function which results in a single dataframe as output to be used later (remember, this is point 4 of our requirements list).

#################################################################
###### Perform simulation
out <- mod %>%
      data_set(data) %>%
      omat(omega) %>%
      smat(sigma) %>%
      mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE) # Simulate 100 observations, independent of the total simulated time
    
    
### Output of simulation to dataframe
df <- as.data.frame(out)

With this simulated dataset df, we can do some summary calculations and plot it in ggplot which will result in the output of our first mrgsolve simulation!

    
#################################################################  
## Calculate summary statistics

sum_stat <- df %>%
    group_by(time) %>%
    summarise(Median_C=median(DV),
              Low_percentile=quantile(DV,probs=minprobs),
              High_percentile=quantile(DV,probs=maxprobs)
    )  
  

################################### Graphical output

plot_pk <- ggplot(sum_stat, aes(x=time,y=Median_C)) +
  
      ## Add ribbon for variability
      geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+
      
      ## Add median line
      geom_line(size=2) +

      # scale_y_log10()+

      # Set axis and theme
      ylab(paste("Concentration",sep=""))+
      xlab("Time after dose (h)")+
      theme_bw()+
      
      # Remove legend
      theme(legend.position="none")

Conclusion

As a conclusion of this first part, we have now developed a script that enables us to perform basic simulations of a 2-cmt PK model in R using the mrgsolve package. We are able to manually adjust the parameter estimates, the dose, and the route of administration. In the end, we have an object holding all the raw simulated values and a graph which shows the median prediction including the variability.

Up next: how to change this script into a Shiny application

Downloads

You can download the R script and the mrgsolve model file below. Make sure you rename the files correctly (remove the .txt extension).

COMMENT

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