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

Shiny applications are a great way to show your complicated models in an interactive way. In recent years, many different examples have been published online showcasing a wide range of functionality. Shiny applications are also being increasingly used in pharmacometrics, to quickly see the effect of a different dosing regimen or even develop a Bayesian TDM application.

---

Getting started

If you are completely new to Shiny, I can recommend visiting a few pages to get the basics:

• https://shiny.rstudio.com/articles/build.html
• https://shiny.rstudio.com/tutorial/

As a basis, a shiny application consists of two parts/scripts. The user interface (UI) and the server script. 

UI script: Consists of all the graphical components needed to interact with the model. This can be buttons, checkboxes, textboxes, labels, graphical output, tables, etc.

Server script: Main script where all the calculations are being performed, being updated by the input from the user. This scripts consists of input components (retrieved from the UI script), reactive expressions (see later), and outputs (such as graphs and tables). We will use the script from Part 1 as the basis for our Shiny application. 

Reactive expressions

In order to change our previous script for PK simulations with mrgsolve to work with Shiny, we need to make use of reactive expressions. Reactive expressions are expressions, ´objects´, that can change based on the user input. For example, the dose, the parameter estimates, the simulation time are all examples of objects that need be included in a reactive expression.

R scripts – ui.R

If you use RStudio, you can easily create a new Shiny Web Application which would automatically generate the ui.R and the server.R script. When you click the Run App button, this default Shiny app is being generated: 

Even though the Old Faithful Geyser Data is interesting, we have different problems to solve… However, this app provides a basis which we can further adapt. Like in part 1 of this blog, we need to define which options we need to have to fulfill all of our criteria. You can check with the R script of part 1 to see how I decided what I need.

Select from a list:

• Oral or intravenous administration

Button:

• Download the simulated dataset

Numerical input:

• Population parameters
• ETA’s and Sigma’s
• Dose
• Number of samples for simulation
• Simulation time
• Min and max probabilities for figure

Graph

• The PK simulation

We can now start updating the sidebarlayout in our ui.R script to incorporate all the components. Previously, it contained a sliderInput to modify the number of bins used in the histogram.

  # Sidebar with a slider input for number of bins 
  sidebarLayout(
    sidebarPanel(
       sliderInput("bins",
                   "Number of bins:",
                   min = 1,
                   max = 50,
                   value = 30)
    ),

Now, we want to change this to include headers, buttons, labels, checkboxes, and numerical input. Let’s start with the dosing details. We need to select the route of administration and be able to set the dose:

  # Sidebar 
  sidebarLayout(
    sidebarPanel(
      h4("Dosing details"),
      selectInput('admin', label='Route of administration', c("Depot","I.V. bolus")),
    
      numericInput("dose", label='Dose', value = 10,min=0)

    ), 

Then, we add the simulation details and a number of numeric inputs for all the parameters. Disclaimer: I fully agree that the design of the interface can be much better but it does the job.

# Sidebar 
  sidebarLayout(
    sidebarPanel(
      h4("Dosing details"),
      hr(),
      selectInput('admin', label='Route of administration', c("Depot","I.V. bolus")),
      numericInput("dose", label='Dose', value = 100,min=0),
      br(),
      h4("Simulation details"),
      hr(),
      numericInput("nsim", label='Number of samples for simulation', value = 1000,min=0),
      numericInput("simtime", label='Simulation time', value = 12,min=0),
      numericInput("minprob", label='Minimal probability', value = 0.25,min=0,max=1),
      numericInput("maxprob", label='Maximal probability', value = 0.75,min=0,max=1),
      
      br(),
      h4("Population parameters"),
      hr(),
      numericInput("ka", label='Absorption rate constant', value = 10,min=0),
      numericInput("v1", label='Central volume', value = 10,min=0),
      numericInput("v2", label='Peripheral volume', value = 10,min=0),
      numericInput("q1", label='Inter-compartmental clearance', value = 10,min=0),
      numericInput("cl", label='Clearance', value = 10,min=0),
      
      br(),
      h4("IIV"),
      hr(),
      numericInput("etaka", label='ETA - Absorption rate constant', value = 10,min=0),
      numericInput("etav1", label='ETA - Central volume', value = 10,min=0),
      numericInput("etav2", label='ETA - Peripheral volume', value = 10,min=0),
      numericInput("etaq1", label='ETA - Inter-compartmental clearance', value = 10,min=0),
      numericInput("etacl", label='ETA - Clearance', value = 10,min=0),
      
      br(),
      h4("Residual variability"),
      hr(),
      numericInput("sigmaprop", label='Proportional', value = 0.01,min=0),
      numericInput("sigmaadd", label='Additive', value = 0.1,min=0)
      
    ), 

Now, in the main panel, we want to have the graph as an output and add the button here where we can download our simulated dataset. Luckily Shiny knows that this is a common feature which is why we can use the downloadButton function.

 mainPanel(
       plotOutput("simPlot"),
       downloadButton("downloadData", "Download")
    ) 

This code means that we need to have two output objects in our server.R script, the simPlot, and the downloadData.

Remember that you need to provide a unique name to all the components that you are adding to your application. We can then call these components from the server script to see the user input.

---

R scripts – server.R

The RStudio template contains just a simple graphical output to start with. No reactive expressions. Nothing complicated.

library(shiny)
# Define server logic required to draw a histogram
shinyServer(function(input, output) {
  
  output$distPlot <- renderPlot({
    
    # generate bins based on input$bins from ui.R
    x    <- faithful[, 2] 
    bins <- seq(min(x), max(x), length.out = input$bins + 1)
    
    # draw the histogram with the specified number of bins
    hist(x, breaks = bins, col = 'darkgray', border = 'white')
    
  })
  
})
 

Now we want to start implementing our previous script into the server side. The way I approach this is by identifying which parts of the code need to be updated when a user changes the input. Unfortunately for us, almost the complete script is dependent on the user input. From the dose that is given to the model that is used.

As a guideline, every reactive expression results in a single dataframe or object that will be updated. In our script, we however need to have two outputs: The dataframe with the simulated values, which we need to be able to download, and the dataframe that will be used for the generation of the graph.

For starters, we create a reactive expression in our script in which we copy the R script until we have a dataframe with the simulations:

  sim_dataframe <- reactive({
    
    ###############################################
    ### Dosing
    
    oral <- T
    
    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
    
    ###################################################
    ############# 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),]
    ######################################################################
    ### 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
      
#################################################################
    ###### 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)
    
    
  }) 

If there is anything unclear about this script, go back to part 1.

We can then change the objects that are currently fixed based on the user input. We can do this by retrieving the information with input$'component name'. For example, input$ka will give the value of the numericinput of the absorption rate constant.

The same piece of code with this modifications can be found below:

  sim_dataframe <- reactive({
    
    ###############################################
    ### Dosing
    if(input$admin =='Depot'){
      oral <- T
    }else{
      oral <- F
    }
    
    dose <- input$dose
    ###############################################
    ## Parameter estimates 
    ka  <- input$ka # Absorption rate constant
    cl  <- input$cl # Clearance
    vd  <- input$v1 # Central distribution volume
    vd2 <- input$v2 # Peripheral distribution volume
    q1  <- input$q1 # Inter-compartmental clearance (set to 0 for 1-cmt model)
    
    ###############################################
    ## Insert etas and sigmas
    
    etaka  <- input$etaka
    etacl  <- input$etacl
    etavd  <- input$etav1
    etavd2 <- input$etav2
    etaq1  <- input$etaq1
    
    #################################################
    ## Insert residual error
    
    sigmaprop <- input$sigmaprop # Proportional error
    sigmaadd  <- input$sigmaadd # Additive error
    
    ###############################3
    ## Simulation info
    
    nsamples <- input$nsim ### Number of simulated individuals
    sim_time <- input$simtime ## Time of simulation
    
    # set probabilities of ribbon in figure
    minprobs=input$minprob
    maxprobs=input$maxprob
    
    ###################################################
    ############# 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),]
    
    
    
    ######################################################################
    ### 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
    
#################################################################
    ###### 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)
    
    return(df)
  }) 

This piece of code will result in the reactive dataframe called sim_dataframe which will be updated every time a user changes the input. If we want to use this dataframe later in our script, we need to use the following syntax: sim_dataframe()

Since objects cannot be transferred between reactive expressions, we need to define the min and maxprobs again here. We use the simulated dataset for the summary calculations as follows:

sum_stat <- reactive({
    
  #################################################################  
  ## Calculate summary statistics
  
# set probabilities of ribbon in figure
  minprobs=input$minprob
  maxprobs=input$maxprob
  sum_stat <- sim_dataframe() %>%
    group_by(time) %>%
    summarise(Median_C=median(DV),
              Low_percentile=quantile(DV,probs=minprobs),
              High_percentile=quantile(DV,probs=maxprobs)
    )  
  return(sum_stat)
  
})
  
   

Notice how I use the sim_dataframe() as the first chunk of the piping function.

Finally, our output, the graph that we want is similar to the code as we developed in part 1, with the only difference being the brackets behind sum_stat(). In this piece of code, we just say that the output$simPlot in the UI, needs to contain the generated ggplot object.

  output$simPlot <- renderPlot({
    
    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")
    
  }) 

So what about our second output, the download button? As I mentioned previously, Shiny is capable of handling these kind of things with limited lines of code. We can specify a default filename and we specify the contents (in our case the sim_dataframe()).

# Downloadable csv of selected dataset ----
  output$downloadData <- downloadHandler(
    filename = function() {
      "Simulated_dataset.csv"
    },
    content = function(file) {
      write.csv(sim_dataframe(), file, row.names = FALSE)
    }
  ) 

Congratulations! You just created your pharmacometric Shiny application! As you may have noticed, starting with the development of a basis script, without any Shiny components, saves a lot of time in the end. Development of a Shiny application can result in some strange error messages and you do not want to waste a lot of time debugging your scripts. 

---

Fancy user interface

In order to get away from the very standard Shiny template, a very interesting package is shinythemes. This allows you to change the look of your application with just the click of a button.

Deploying your Shiny application for free

Information on how to make your app publicly available can be found on the following website:

https://shiny.rstudio.com/deploy/

The free package on shinyapps.io allows 25 hours of use per month.

Conclusion

Lets take a look of what the functionality of the Shiny app should be, which we specified in Part 1:

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

As you now can see from your own application, we fulfilled all the points! Of course, there are many remaining items that can be improved. For example, now, when a single number is being changed, the whole simulation is being repeated. You can add an action button to your application which holds all actions until you press it (especially for large simulations this makes sense). Also, I load the model code in the reactive expression, which means that it is being loaded on every change (whereas the model remains the same...). Furthermore, you can write a code to generate pdf reports from your app, by using RMarkdown. The options are limitless.

If you know the basics, it is just a small step from the app that we have now created to the one that is located on my website. Good luck!

You can download the created ui.R and server.R scripts below

COMMENT

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