This two part series will show you how to create a simple pharmacometric Shiny application using the mrg*solve* 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:

- Create the R code needed to run a single simulation
- 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 mrg*solve* 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!

### mrg*solve*

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

The homepage of mrg*solve*, 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 mrg*solve*. 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:

- Pharmacokinetic simulations of a simple 2-compartment model with oral or I.V. administration
- Manually specify the parameter estimates (theta’s, omega’s and sigma’s)
- Generate a figure of the simulation
- 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.

### mrg*solve* 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 mrg*solve* 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

#### mrg*solve* 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 mrg*solve* 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 mrg*solve* 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 mrg*solve* 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!

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

MJvanEsdonkPost authorDone! Thanks for your comment!

RyanDear MJVANESDONK,

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

MJvanEsdonkPost authorHello,

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?

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

MJvanEsdonkPost authorThere 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)

RyanThank you, That worked !!!

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

MJvanEsdonkPost authorDear,

Do you mean bioavailability?

You can set it with F_GUT as a parameter for example.

There is a lot of information on all mrgsolve options here, and on section 2.3.15 on the F:

https://mrgsolve.github.io/user_guide/model-specification.html

RyanThank you very much again for your help. I greatly appreciate it.

ThiDear MJVANESDONK,

I really appreciate your post! Your guidance is very much helpful for a junior like me.

Many thanks,

Thi

WilliamHi,

Great tutorial! However, I ran into problems regarding smat and omat. When I ran the code in your tutorial, after the mrgsim part it showed

‘Error in h(simpleError(msg, call)) :

error in evaluating the argument ‘.x’ in selecting a method for function ‘smat’: improper signature: omat’

I am not sure how to troubleshoot this.

Thank you in advance!

MJvanEsdonkPost authorHi,

Thank you for your comment.

Unfortunately I am unable to reproduce your error. The error does exist when the omega object is not a correct matrix or does not exist in the environment.

Could you check what response you get when you only try to execute omat(omega)? This should print a matrix in the console.

Please make sure you have the latest mrgsolve version and the latest version of R.

Michiel

drewQuestion: if you have a nonzero additive error component, then when you simulate, your initial concentration will be nonzero. Does that mean in simulations, you should ignore that residual error component?

MJvanEsdonkPost authorHello Drew, that really depends on the question you would like to answer if the simulation should be done with the residual error component.

Your question is valid. But if a concentration is 0, and it is being analysed, will it result in a non-zero concentration? If not, there might not be an additive error component present.

Another approach is to fix simulated values below the lower limit of quantification to 0, after the simulation is finished.

Michiel