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:

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!: 

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")
PK simulation in mrgsolve

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!

9 thoughts on “Creating a simple pharmacometric Shiny application with mrgsolve in R – Part 1

  1. David Busse Reply

    Great and well structured tutorial! Thanks a lot for the efforts. Would be great to have a link to part 2 on this page.

  2. Ryan Reply

    Dear MJVANESDONK,

    When i changed the operator to F (for Oral dose) the system is hanging up. can you help?

    • MJvanEsdonk Post authorReply

      Hello,

      The object ‘oral’ can be set to TRUE or FALSE. If you set to true:
      oral <- TRUE It should work. If not, from what section onwards is it hanging up?

  3. Ryan Reply

    Thank you for a quick reply. One clarifications is i have modified the model to 3 compt. The code is hanging up at perform simulation. I am pasting my mrgsolve code for 3cmpt and R code for running. Apologies for a long post. Your help is greatly appreciated.

    $PARAM

    TVKA=1, TVCL=1, TVVC=20,TVQ1=2,TVVP1=10,TVQ2=2, TVVP2=1

    $CMT

    EV CENT P1 P2

    $MAIN

    double KA = TVKA*exp(ETA(1));

    double CL = TVCL*exp(ETA(2));

    double VC = TVVC*exp(ETA(3));

    double Q1 = TVQ1*exp(ETA(4));

    double VP1 = TVVP1*exp(ETA(5));

    double Q2 = TVQ2*exp(ETA(6));

    double VP2= TVVP2*exp(ETA(7));

    double k10=CL/VC;

    double k12=Q1/VC;

    double k13=Q2/VC;
    double k21=Q1/VP1;
    double k31=Q2/VP2;

    $OMEGA 0 0 0 0 0 0 0

    $SIGMA @labels PROP ADD
    0 0

    $ODE

    dxdt_EV = -KA*EV;

    dxdt_CENT = KA*EV+k21*P1+P2*k31-(k12+k13+k10)*CENT;
    dxdt_P1 = EV*k12 – P1*k21;

    dxdt_P2 = EV*k13 – P2*k31;

    $TABLE

    double IPRED = CENT/VC;

    double DV = IPRED*(1+PROP)+ADD;

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

    }

    $CAPTURE IPRED DV

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

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

    IM <- F

    dose <- 2 # mg

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

    etaka <- 0.09
    etacl <- 0.09
    etavd <- 0.09
    etavd2 <- 0.09
    etavd3 <- 0.09
    etaq1 <- 0.09
    etaq2 <- 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

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

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

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

    ######################################################################
    ### Load in model for mrgsolve
    mod <- mread("3cmpt")

    ## Specify the omegas and sigmas matrices
    omega <- cmat(etaka,
    0, etacl,
    0,0,etavd,
    0,0,0,etavd2,
    0,0,0,0,etavd3,
    0,0,0,0,0,etaq1,
    0,0,0,0,0,0,etaq2)

    sigma <- cmat(sigmaprop,
    0,sigmaadd)

    ## Set parameters in dataset
    data$TVKA <- ka
    data$TVCL <- cl
    data$TVVC <- vd
    data$TVVP1 <- vd2
    data$TTVP2<- vd3
    data$TVQ1 <- q1
    data$TVQ2 <- q2

    #################################################################
    ###### Perform simulation
    out %
    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)

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

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

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

    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”)

    • MJvanEsdonk Post authorReply

      There are 2 reasons why it is not working. The most important part is that the differential equations in your mrgsolve model are incorrect, you should change them to:

      $ODE

      dxdt_EV = -KA*EV;

      dxdt_CENT = KA*EV+k21*P1+P2*k31-(k12+k13+k10)*CENT;
      dxdt_P1 = CENT*k12 – P1*k21;

      dxdt_P2 = CENT*k13 – P2*k31;

      Second, it looks like a part of the dplyr piping function went missing. There is no ‘mod %>%’ in your code. When I use the code below it works:

      #################################################################
      ###### Perform simulation
      out <- mod %>%
      data_set(data) %>%
      omat(omega) %>%
      smat(sigma) %>%
      mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE)

  4. Ryan Reply

    Dear MJVANESDONK, can you help me with one curious question. How do we incorporate fraction absorbed in the simulation?

Leave a Reply

Your email address will not be published. Required fields are marked *