Note: sometimes, the embedded shiny apps don’t show up on Github. Please use this link as a temporary work around in such cases.

Introduction

The following examples show how to control dynamic models implemented in the R language interactively from a web interface. We use a technically simple ordinary differential equation (ODE) model as an example, then we will see how to add GUI elements and finally how we can make it faster by using compiled code.

This tutorial was originally presented at the useR!2017 conference in Brussels.

The HTML version embeds the shiny apps from a remote server and needs no local installation. The source code of this HTML file is in Rmarkdown (.Rmd) format that can be run locally from within RStudio when the shiny apps locally installed.

Note: this is an interactive document. Please use the Code buttons to unfold the code of individual code cunks or for the complete document. The content is best viewed with Chrome or Firefox.

Ingredients needed

  1. A recent installation of the R language
  2. A dynamic deSolve model
  3. The shiny package and an R script app.R
  4. A shiny server

The Brusselator

The Brusselator (Brussels Oscillator) a mathematical model of an autocatalytic reaction of two substances A and B, two intermediate products X and Y and two final products C and D. The latter are continuesly removed.

  1. A → X
  2. B + X → Y + C
  3. 2X + Y → 3X
  4. X → D
  • developed by Ilya Prigogine and R. Lefever 1967,68 at the Université Libre de Bruxelles
  • important theoretical example for studying oscillations and bifurcation
  • see Wikipedia for details
\[ \begin{aligned} \frac{dX}{dt} &= k_1 A - k_2 B X + k_3 X^2 Y - k_4 X \\ \frac{dY}{dt} &= k_2 B X - k_3 X^2 Y \end{aligned} \]

The Brusselator in R

Implementation of the model in R is straightforward. We need the package deSolve, translate the differential equations of the Brusselator into an R function, define initial values (\(y_0\)), parameters and time steps and then use ode for the numerical solution. Finally we plot the results.

library("deSolve")

brusselator <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dX <- k1*A   - k2*B*X    + k3*X^2*Y - k4*X
    dY <- k2*B*X - k3*X^2*Y
    list(c(X=dX, Y=dY))
  })
}

y0 <- c(X=1, Y=1)
parms <- c(A=1, B=3, k1=1, k2=1, k3=1, k4=1)
out <- ode(y = y0, times=seq(0, 100, .1), brusselator, parms)
plot(out)

A minimalistic Shiny app

library("shiny")

server <- function(input, output) {
  output$distPlot <- renderPlot({
    hist(rnorm(input$obs))
  })
}

ui <- fluidPage(
    sliderInput("obs", "Number of observations:", min = 10, max = 500, value = 100),
    plotOutput("distPlot")
)

shinyApp(ui = ui, server = server)

The second part creates the web interface. The app consists of a server element doing the computations and plotting, and a user interface (ui), that is responsible for handling user inputs and for sending the outputs back to the user. Both parts, the server and the ui can be saved in a single file ‘app.R’. It can then be run directly from RStudio or placed on a shiny-enabled web server.

Exercise: Copy the app to Rstudio, save it to the disc (fille app.R) and try it out

More about this: https://shiny.rstudio.com

Combine the parts

Now let’s combine the shiny ui/server environment with our Brusselator model. We start with a very simple application that allows to change two parameters A and B (essentially the start concentrations of substances A and B) and then plots the time series of the model.

A first approach

The code is stunningly simple, consisting of three parts:

  1. The model from above (without any modifications).
  2. The server component. It consists of a single function getting input and output. The two parameters A and B are taken from the input and the plot is then directed to an element output$brusselsof the output. Note also that we use deSolve’s function matplot.0d instead of plot as before, so that both state variables appear in a single figure.
  3. The ui calls then a shiny element fluidpage with two numericInput calls and one call to plotOutput.
library("deSolve")

brusselator <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dX <- k1*A   - k2*B*X    + k3*X^2*Y - k4*X
    dY <- k2*B*X - k3*X^2*Y
    list(c(X=dX, Y=dY))
  })
}

server <- function(input, output) {
  output$brussels <- renderPlot({
    parms <- c(A=input$A, B=input$B, k1=1, k2=1, k3=1, k4=1)
    out <- ode(y = c(X=1, Y=1), times=seq(0, 100, .1), brusselator, parms)
    matplot.0D(out)
  })
}

ui <- fluidPage(
  numericInput("A", label = "A", value = 1),
  numericInput("B", label = "B", value = 3),
  plotOutput("brussels")
)

shinyApp(ui=ui, server=server)

Improvments

Now, we can start to adapt the app to our needs: change the layout, unlock more parameters of the model, add a second figure and so on. An especially convenient layout is the sidebarLayout that separates the control elements (sidebarPane left or top) and output of the app (mainPanel right or bottom).

HTML tags like h2, h3, or p can be used to annotate the app, see help pages for more.

library("deSolve")

brusselator <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dX <- k1*A - k2*B*X + k3*X^2*Y - k4*X
    dY <- k2*B*X - k3*X^2*Y
    list(c(X=dX, Y=dY))
  })
}


server <- function(input, output) {
  output$brussels <- renderPlot({
    y0 <- c(X=input$X, Y=input$Y)
    parms <- c(A=input$A, B=input$B,
               k1=input$k1, k2=input$k2, k3=input$k3, k4=input$k4)
    times <- seq(0, 100, .1)
    out <- ode(y0, times, brusselator, parms)
    par(mfrow=c(1, 2))
    matplot.0D(out, main="time series")
    plot(out[,2:3], type="l", xlab="A", ylab="B", main="state diagram")
  })
}

ui <- fluidPage(
  headerPanel("Brusselator"),
  sidebarLayout(
    sidebarPanel(
      h3("Init values"),
      numericInput("X", label = "X",
                   min = 0.0, max = 5,  value = 1, step = 0.2, width=100),
      numericInput("Y", label = "Y",
                   min = 0.0, max = 5,  value = 1, step = 0.2, width=100),

      h3("Parameters"),
      numericInput("A", label = "A",
                   min = 0.0, max = 5,  value = 1, step = 0.2, width=100),
      numericInput("B", label = "B",
                   min = 0.0, max = 5,  value = 3, step = 0.2, width=100),
      numericInput("k1", label = "k1",
                   min = 0.0, max = 10, value = 1, step = 0.2, width=100),
      numericInput("k2", label = "k2",
                   min = 0.0, max = 10, value = 1, step = 0.2, width=100),
      numericInput("k3", label = "k3",
                   min = 0.0, max = 10, value = 1, step = 0.2, width=100),
      numericInput("k4", label = "k4",
                   min = 0.0, max = 10, value = 1, step = 0.2, width=100)
    ),
    mainPanel(
      h3("Simulation results"),
      plotOutput("brussels")
    )
  )
)

shinyApp(ui = ui, server = server)

The code above should work as intended, but there are of course further options of improvement.

  1. Improve layout and design of the app: this can be quite time consuming and wil not be handled here.
  2. Generalize the app, so that new models can instantly be implemented.
  3. Make the model faster.

Here, we want to focus on options 2 and 3.

Generalization

R is a programming language, so why not creating the user interface directly from information found already in the initial values and parameters of the model? This would allow us to save stupid typing (and copy and paste!) and to adapt the user interface ``automagically’’ if we change something in the model.

The solution is, to construct lists of lists beforehand, in our case one for the parameters L_parms and one for the initial values L_y0. Each element of the top level list contains all the parameters needed by a call to numericInput, including some constant optional values like the witdh of the input field. The templates stored in these lists are then expanded with lapply and do.call. The result is exacly the same as before, compare: https://weblab.hydro.tu-dresden.de/apps/bruss03

library("deSolve")

brusselator <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dX <- k1*A - k2*B*X + k3*X^2*Y - k4*X
    dY <- k2*B*X - k3*X^2*Y
    list(c(X=dX, Y=dY))
  })
}

## prepare data structures to create UI programmatically
y0    <- c(X=1, Y=1)
parms <- c(A=1, B=3, k1=1, k2=1, k3=1, k4=1)

makelist <- function(i, obj, min=NA, max=NA, step=NA, width=NULL) {
  list(inputId=names(obj[i]), label=names(obj[i]),
       value=unname(obj[i]), min=min, max=max, step=step,
       width=width)
}

## two lists of lists
L_parms <- lapply(1:length(parms), makelist, obj=parms, min=0, max=10, step=0.2, width=100)
L_y0 <- lapply(1:length(y0), makelist, obj=y0, min=0, max=10, step=0.2, width=100)

server <- function(input, output) {
  output$brussels <- renderPlot({
    L_input <- reactiveValuesToList(input) # to enable width
    y0    <- with(L_input, c(X=X, Y=Y))
    parms <- with(L_input, c(A=A, B=B, k1=k1, k2=k2, k3=k3, k4=k4))
    times <- seq(0, 100, .1)
    out <- ode(y0, times, brusselator, parms)
    par(mfrow=c(1, 2))
    matplot.0D(out, main="time series")
    plot(out[,2:3], type="l", xlab="A", ylab="B", main="state diagram")
  })
}

ui <- fluidPage(
  headerPanel("Brusselator: generic UI creation"),
  sidebarLayout(
    sidebarPanel(
      ## generic creation of UI elements
      h3("Init values"),
      lapply(L_y0, function(x) do.call("numericInput", x)),   # <--------

      h3("Parameters"),
      lapply(L_parms, function(x) do.call("numericInput", x)) # <--------
    ),
    mainPanel(
      h3("Simulation results"),
      plotOutput("brussels")
    )
  )
)

shinyApp(ui = ui, server = server)

Make it faster

As a second improvement of our app, we want to make our model faster. Quick response of a server based app is crucial when models become bigger, because our users may think that the connection is broken otherwise.

This is fortunately possible, given that some development tools (e.g. the Rtools on Windows) are installed.

Package deSolve supports models in compiled code provided by a .dll / .so directly. It allows the solver to communicate directly with the model, and without any overhead of a back-call to R. Our experience shows, that this can then be 2 … 200 times faster, depending on the type of model.

Several approaches were developed, depending on model complexity and your skills:

  • write your model ``by hand’’ in Fortran or C
  • use a code generator:
    • package ccSolve from Karline Soetaert (from Github)
    • package cOde from David Kaschek for smaller models
    • package rodeo from David Kneis for complex models

We use cOde here, because it is especially convenient for small models, while for more complex models, the table based notation of rodeo can be more advantageous.

The Brusselator in C

Package cOde allows to compile our (technically!) simple Brusselator without knowledge of Fortran or C. The only prerequisite is, that we have the development tools (the compiler etc. …) installed, that the system PATH correctly set, and that write access to the working directory is granted.

Then it is sufficient, to write the differential equations in a character vector and then to compile it with funC. This creates a C file and compiles it into a shared library (.dll or .so), which then can be numerically solved with odeC.

library("deSolve")
library("cOde")    # 3rd party package by Daniel Kaschek

y0 <- c(X=1.1, Y=1)
times <- seq(0, 100, .1)
parms <- c(A=1, B=3, k1=1, k2=1, k3=1, k4=1)

brussC <- funC(c(
    X = "k1*A - k2*B*X + k3*X^2*Y - k4*X",
    Y = "k2*B*X - k3*X^2*Y"
), modelname="brusselator", compile=TRUE)

out <- odeC(y0, times, brussC, parms)
matplot(out[,1], out[,2:3], type="l")

Shiny app with compiled code

As a final exercise, we modify our app from before, to use the compiled C version of the model instead of the R version.

library("deSolve")
library("cOde") # 3rd party package by Daniel Kaschek

## helper function, adapted from package cOde v.2.2 (GPL >= 2.0)
loadDLL <- function(func, cfunction="derivs") {
  .so <- .Platform$dynlib.ext
  checkDLL <- try(getNativeSymbolInfo(cfunction), silent=TRUE)
  if(inherits(checkDLL, "try-error")) {
    dyn.load(paste0(func, .so))
    #cat("Shared object is loaded and ready to use\n")
  } else if ((is.null(checkDLL$package))) {
    # We are on Windows (always overload)
    dyn.load(paste0(func, .so))
  } else if((checkDLL$package)[[1]] != func) {
    # We are on Unix
    dyn.load(paste0(func, .so))
  }
}

compile <- !file.exists(paste0("brusselator", .Platform$dynlib.ext))

brussC <- funC(c(
  X = "k1*A - k2*B*X + k3*X^2*Y - k4*X",
  Y = "k2*B*X - k3*X^2*Y"
), modelname="brusselator", compile=compile)

if (!compile) loadDLL("brusselator")

y0 <- c(X=1, Y=1)
parms <- c(A=1, B=3, k1=1, k2=1, k3=1, k4=1)

makelist <- function(i, obj, min=NA, max=NA, step=NA, width=NULL) {
  list(inputId=names(obj[i]), label=names(obj[i]),
       value=unname(obj[i]), min=min, max=max, step=step,
       width=width)
}

L_y0    <- lapply(1:length(y0), makelist, obj=y0, min=0, max=10, step=0.2, width=100)
L_parms <- lapply(1:length(parms), makelist, obj=parms, min=0, max=10, step=0.2, width=100)

server <- function(input, output) {
  output$brussels <- renderPlot({
    L_input <- reactiveValuesToList(input)
    y0    <- with(L_input, c(X=X, Y=Y))
    parms <- with(L_input, c(A=A, B=B, k1=k1, k2=k2, k3=k3, k4=k4))
    times <- seq(0, 100, .1)
    #out <- ode(y0, times, brusselator, parms) # R version
    out <- odeC(y0, times, brussC, parms)      # C version
    par(mfrow=c(1, 2))
    matplot.0D(out, main="time series")
    plot(out[,2:3], type="l", xlab="A", ylab="B", main="state diagram")
  })
}

ui <- fluidPage(
  headerPanel("Brusselator"),
  sidebarLayout(
    sidebarPanel(
      h3("Init values"),
      lapply(L_y0, function(x) do.call("numericInput", x)),
      h3("Parameters"),
      lapply(L_parms, function(x) do.call("numericInput", x))
    ),
    mainPanel(
      h3("Simulation results"),
      plotOutput("brussels")
    )
  )
)

shinyApp(ui = ui, server = server)

The result looks identical as before, just a little bit faster.

If you want to see the C code: here is it

Further reading

Original author: tpetzoldt