Interactive Pharmacometric Applications Using R and the Shiny Package

A population pharmacokinetic and/or pharmacodynamic model can be thought of as an elegant and concise description of the underlying mechanisms and variability in a data set. A robust simulation model is an important tool for predicting new data—the consequences of changes in the dosing regimen or covariate influences. Population models have found important roles in the development, regulation, and optimal use of pharmaceuticals.1, 2 However, the process of making predictions from population models has been largely time-consuming, and remained the province of dedicated pharmacometricians using specialized software limiting its wider applicability.3

Recent developments such as the ggplot24, 5 package for the R data analysis and statistical language6 have allowed sophisticated and flexible plotting of data and population model output. However, the plots are static and to investigate different values for model parameters the models need to be re-simulated and plots updated manually. The Berkeley Madonna software provides an alternative approach where models can be specified as differential equations and the simulated results shown in real-time for different parameter values by use of sliders.7 While Berkeley Madonna has found useful roles in presenting model predictions to nonpharmacometricians, it is proprietary software (albeit inexpensive) and is not readily adapted to simulations of variability.

Recently, developments in the R language, and in particular the Shiny package for R,8 have allowed R programmers to interactively show the output for R programs to Web-browsers. Given the general nature of the R language, it is possible to program interactive pharmacometric models with a Shiny Web-browser interface that can be viewed on the localhost (user's own computer) or on another computer accessed by means of the Internet. Applications such as a dosing-education tool aimed at high school students (https://acp-unisa.shinyapps.io/OpenDay/) and a population model simulation tool complete with simulated variability (https://acp-unisa.shinyapps.io/Ibuprofen/) have been developed using the Shiny package and R and can be viewed without an R installation or files containing R code. Unlike other Web-page design methods, only previous experience with the R programming language is required. R and Shiny maintain Berkeley Madonna's key feature of "reactively" updating output in response to changing input by means of widgets (such as sliders and radio buttons) but owing to the flexible R language and combination with extension packages, the pharmacometrician has control in coding all elements of a population model, the appearance of the application's user-interface, and generated output. We believe such Shiny Web applications have the potential to provide a convenient method of providing model simulations to a broad audience, which may be useful to pharmacometricians and nonpharmacometricians alike.

GETTING SHINY

This tutorial is primarily targeted toward intermediate/experienced R programmers such as those who use R for processing raw data and NONMEM^® output, statistical analyses, and have some experience in coding pharmacometric models in R (nonetheless, an annotated example script for coding a population model in R is available as Supplementary Pharmacometric Model Example). It is not recommended to learn R and Shiny concurrently. R is an open source data analysis language and can be downloaded from http://www.r-project.org. There are several online resources for learning R. Shiny is a package for R developed by RStudio, and needs to be installed in a given R installation. There are several ways to install packages in R depending on the operating system and R interface. One convenient method is to install the RStudio integrated development environment for R (http://www.rstudio.com/) and use the Tools/Install Packages menu. RStudio will automatically install any additional package dependencies. A list of software and R package versions used in the development of this tutorial's examples is available in Table 1. For those less familiar with R, it is important to note that R and its package libraries develop and update at a fast pace. Therefore, in the future the latest versions of the packages used in this tutorial may not be backward compatible with the code presented here.

USING SHINY

Shiny applications are built using two R scripts that communicate with each other: a user-interface script (ui.R), which controls layout and appearance; and a server script (server.R), incorporating instructions for user-input, processing data, and output by utilizing the R language and functions from user-installed packages. Tutorials and exercises for building Shiny applications are provided by RStudio on the Shiny website9 and are accompanied by articles that describe capabilities of Shiny beyond the scope of the tutorials, a reference page for Shiny functions, as well as a gallery of examples supported by code. The Shiny package also includes 11 built-in examples referred to by the Shiny by RStudio tutorials.

To run any of the applications locally in R or RStudio, an installation of the Shiny package (and any dependencies) is required, and ui.R and server.R scripts for the application saved in the same directory (i.e., a folder titled with the application's name). To launch the application from RStudio open each ui.R and server.R script in RStudio and click “RunApp,” which will appear in the top right hand corner of the source pane. Launching an application from R requires setting the working directory to where the application folder is located, and using the runApp function. When initiating, Shiny will open a Web-browser window for the application.

STRUCTURE OF A SHINY APPLICATION

In this section, we will discuss the core elements of the ui.R and server.R scripts required to create an application. Our first example for creating a user-interface in ui.R and application instructions in server.R is the Flip-Flop Kinetics application—Example 1 (supplementary code available online). This application implements a model for a hypothetical drug described by one-compartment first-order absorption kinetics. It illustrates the concept that for this type of model, there are two parameter sets that give identical results: one with absorption rate faster than elimination rate, and the other with elimination rate faster than absorption rate (flip-flop kinetics).10

Figure 1 features a screenshot of the application's browser window with a plot of the concentration–time profile for the drug over a 24-h period, with sliders that control the values for the absorption and elimination rate constants, k_a and k_e, respectively, and volume of distribution, V. Initially (Figure 1), the application's plot shows an example where elimination is rate limiting as the value of k_a is greater than k_e (k_a = 2, k_e = 0.2). When the sliders change, thus allowing the values of the rate constants to approach each other, the plot automatically updates in response. When k_a is less than k_e, the roles of the constants in describing the concentration–time profile swap and the absorption rate constant limits the decline of drug. For example, pairs of values such as k_a = 0.2 and k_e = 2 when flipped to correspond with the other rate constant (and V is also changed accordingly), do not change the outcome of the plot (Figure 1). Rather than a series of simulations and static plots to explore this system, Shiny has allowed the user to control and observe the roles of the rate constants simultaneously.

The ui.R script (supplementary code available online) encodes instructions for the application's layout, appearance, user-input widgets (interactive Web elements such as sliders, buttons, selection boxes, check boxes, etc.), and the output to be displayed. The main elements describing the user-interface of this application are shown below:

All code for the contents of the user-interface is required to be within the brackets of a layout function. Layout functions such as fluidPage (there are a number available) create a canvas for the interface and use fluidRow to position user-input widgets (sliderInput—a function for creating a slider) and output (plotOutput—function for plot object). Each layout function has its own framework for positioning elements where others, including fixedPage and navbarPage, can create different styled pages based on their own functions. However, they all follow the same hierarchical structure where element functions (such as widgets—sliderInput) are embraced within a positioning function (such as fluidRow), within a layout function (fluidPage). Functions of the same level are written as a string within their superior function, separated by “,”.

Error messages related to Shiny (and other loaded R packages) appear in the loaded Web browser upon application initiation, and once the application is terminated they also appear in the R Console. Examining if pairs of function brackets are closed before running the application can help limit a large number of error messages. In particular, functions placed earlier in the script are more difficult to track when ui.R becomes more detailed. Using free source code editor software, or directly writing scripts from RStudio, is encouraged when writing Shiny applications as they can aid in highlighting any unclosed brackets (among other navigational and editing benefits). If there is a major error in which the application is not functional a Web-browser page will still open, however, will be grayed-out. Other error messages generally provide a line number referring to the source code or the function name in question. When testing the coding of a pharmacometric model, it is recommended to write a generic R script ensuring that it runs successfully before incorporating the model into a Shiny application.

The structure of server.R code plays an important role in defining instructions for the application, while minimizing redundant computation and maximizing application speed (Figure 2). The function, shinyServer, requires an input object (values called from ui.R) and an output object (objects to be called by ui.R). Objects that change depending on input from widgets in the ui.R, such as input$KA in this example, are termed “reactive.” Every time input from a widget changes, the value for the reactive object updates to reflect this change.

Expressions for defining and processing reactive objects are required to be enclosed within a render* function (for returning reactive output to the user-interface, where * denotes an object description such as “plot,” “table,” or “text”) or a reactive expression (often for controlling a, or series of, reactive data frames that can be sent to the user-interface by render* functions). The renderPlot function (Figure 2) containing reactive input objects (KA, KE, and V from sliders), expressions for calculating concentration (CONC), and a ggplot25 object plotting resultant concentrations over time (plotobj), will re-execute with every reaction to a widget change and will store the updated plot object as plotCONC for the output object. Detailed code can significantly slow down the application; therefore, it is best to limit code within render* and shinyServer functions to just that is reactive or computationally simple. Code at the beginning of the script before shinyServer is only run once on application initiation. This area is the most appropriate place for loading libraries, source code, datasets, and constant expressions (such as defining time) that do not need to be called each occasion when input from a widget changes, thus limiting redundant computation time. These functions and libraries can also be stored in a third script called global.R.

BUILDING A USER-INTERFACE (ui.R)

Example 2 is a Shiny application that incorporates several layout functionalities and widgets in its ui.R script (supplementary code available online). Shiny's variety of customizable layouts and prebuilt widgets allows easy development of user-interfaces for applications. From the options available, developers can explore layouts that adapt to the different browser sizes of devices (i.e., phone, tablet, or computer), add a sidebar or tabs that organize and differentiate between input (widgets) and output (tables, plots, etc.), or change what is visible when specific input conditions are met. The example application implements a simple patient education tool developed for an Open Day at a university.

Figure 3 features a screenshot of the application's browser window with a plot of the 10-d concentration–time profile for a hypothetical orally administered drug with one-compartmental first-order absorption kinetics affected by patient age (on clearance) and weight (on volume of distribution). When the slider values for age, weight, and dose change, or when a dose frequency is selected, the plot of the concentration–time profile is updated. The application was intended to provide final year high school students and their parents with no background in pharmacy or pharmacokinetics an understanding of why some drugs require altered dosing regimens for different ages and weights to prevent subsets of patients experiencing toxicity or lack of efficacy. The students and their parents were able to also explore profiles with missed or doubled doses by means of integrated check boxes.

DEFINING APPLICATION INSTRUCTIONS (server.R)

The server controls the processing of user-input to display output to the user-interface. For those experienced with the R language, writing server.R consists of using the same packages, functions, and arguments to process information as previously in R with the only difference being considering whether specific tasks need to be placed in the reactive part of the script.

The server.R script (Example 2) is divided into two sections: an area for calling and processing reactive input and creating reactive output enclosed within shinyServer, and an area for nonreactive expressions and functions not dependent on widget-input (outside of shinyServer). The premise for the arrangement is based on separating expressions and functions that need to be re-evaluated upon input changes and those that never need to. All code before shinyServer (loading package libraries, defining a time sequence, ggplot25 themes and the function for a differential equation system) will be executed once on application initiation and will not be re-executed unless the Web-browser for the application is refreshed or the application is re-run from R. The results of these commands are not dependent on changes from input, thus their position outside of shinyServer limits the re-evaluation of nonreactive code. On the other hand, the reactive expression and the render* function within shinyServer of this example are re-executed on every change of a widget as they are dependent on widget input.

Reactive expressions

Some applications require more complex processing of multiple reactive objects and data, that when placed inside render* braces significantly slow the application's re-evaluation of the output. Furthermore, there may be instances where a reactive data frame may need to be called by multiple render* functions, which cannot be achieved when enveloped in a single one. In this application (Example 2), a data frame of time and concentrations at each time-point is constructed from reactive objects controlled by user-input widgets surrounded by reactive braces (just like a render* function). The content of the reactive data frame is used to produce a reactive plot of the concentration–time profile (renderPlot). A skeleton of the reactive expression of the server.R script for this application is provided below:

The reactive expression uses a series of R expressions that use widget input to return an object, or in this case a data frame—sim.data.df, that updates whenever widgets change. As discussed in Example 1, reactive objects are required to be enclosed within a render* function or a reactive expression. Therefore, in Example 2, all user-input widget values that are involved in the creation of sim.data.df are enclosed within the reactive function defining sim.data. In addition to widget values such as input$X, all expressions dependent on reactive input for defining sim.data are required to be within the reactive braces—this includes “if” statements, differential equations (as seen in Example 3), installed package functions, and standard R functions. However, reactive expressions cannot send their object to the user-interface; therefore, doing so requires a render* function, such as renderPlot.

Sending output to the user-interface

Objects appearing in the user-interface reactive to widget input or dependent on a reactive data frame (such as sim.data.df) are built into the output object by shinyServer. An element, such as plotCONC, defining a ggplot25 object is defined within the shinyServer braces (outside the reactive expression) as:

Each output object requires a render* function that corresponds to the type of reactive object needed to be made. In pharmacometrics, data are often in the form of a table, which is illustrated by a plot or further analyzed to obtain a statistical summary. Therefore, renderPlot, renderText, and renderTable are commonly required to display a reactive plot, text or table, respectively. Example 2 uses renderPlot to create plotCONC using the data frame—sim.data.df, where sim.data is called using its name followed by parentheses. As sim.data is reactive, it will check if the widget input has changed, update the data, and renderPlot will re-draw the plot accordingly.

Here, the ggplot25 package is used to create the plot object. The ggplot argument initializes a ggplot object and states the input data frame, from which defined variables can be used to create a plot. Although not shown as part of Example 2, creating a table object may require further processing of sim.data that can be achieved by an R function with sim.data as the input:

This demonstrates a scenario where sim.data can be simultaneously used to create multiple, or combinations of, reactive plots and tables using a single reactive expression, which could not be achieved if all code was embraced in a single render* function. For creating a reactive text object, the structure of the previous code can be followed using renderText in place of renderTable in server.R.

SIMULATING POPULATION MODELS WITH R AND SHINY

Shiny applications can use R packages such as deSolve12 (differential equation solver) to represent population models as differential equations. A population model in Shiny allows stochastic simulation of model predictions that embody, for example, between subject variability. Traditionally, this has been achieved by simulating, for example, concentration–time profiles for several individuals administered a dosing regimen (i.e., using NONMEM^®),13 plotting the outcome (often summarized by a median and confidence intervals) and saving the image, and repeating the process for a new dose. The plots created are static and the process of simulation can be time-consuming thus rationalization of the number and types of regimens tested is often required. Shiny applications allow for graphical output to be produced simultaneously with changing input (such as dosing regimen or covariate values); therefore, the user can rapidly view the results of candidate scenarios. Although other model visualization software such as Berkeley Madonna can also simultaneously update output, they are not readily adapted to simulate variability.

The third example Web application (Figure 4) demonstrates how stochastic simulation in Shiny can be used to explore a series of ibuprofen dosing regimens for patent ductus arteriosus (PDA) closure in preterm neonates less than 32 weeks gestation (Example 3; supplementary code available online). The application embodies a population pharmacokinetic model of ibuprofen enantiomers, R- and S-ibuprofen, in preterm infants14 described as a series of differential equations. The model includes one-compartment models for each of R- and S-ibuprofen with unidirectional bioconversion of R-ibuprofen into S-ibuprofen and the effect of increasing postnatal age increasing elimination of the R-enantiomer.14 Using the R programming language, the application simulates a patient population from random sampling of the parameter distributions of the population model to provide each individual with their own parameter set. Parameter sets and the differential equation system are provided as input to a differential equation solver to calculate the concentration–time profiles for R- and S-ibuprofen from 0 to 72 h at 15-min increments for each individual of the population. From the resultant data frame, the median and lower and upper percentiles for S-ibuprofen concentrations (red solid line and shaded ribbon, respectively) and median R-ibuprofen concentrations (blue solid line) are calculated, plotted and displayed in the user-interface (Figure 4).

The model embodies important information about the pharmacokinetics of each enantiomer in neonates, from which changes in the concentration–time profile can be examined and re-examined when input doses or regimens change. This application incorporates radio buttons to represent two dosing regimens (SELECT). The “loading-bolus-bolus” regimen reflects the most common pharmacotherapeutic management for PDA closure in preterm neonates consisting of three intravenous (IV) ibuprofen doses administered every 24 h; a 10 mg/kg IV loading dose on day one (LDOSE, often within first 12 h of life), and two subsequent 5 mg/kg IV doses (BDOSE) on days two and three (i.e., 10-5-5 mg/kg). A second course of a higher dosed regimen (20-10-10 mg/kg) is also commonly administered in individuals with PDA resistant to the first course,15 which can be simulated by adjusting the “Loading Dose” and “Bolus Dose” sliders. “Loading-continuous infusion” is a regimen which has not been commonly implemented in this practice. This Shiny application can explore the loading dose (LDOSE) and 72-h continuous infusion rate (CDOSE) required to achieve 90% of S-ibuprofen concentrations above the IC50 for cyclooxygenase-2 (COX-2; 16.5 μg/ml),16 and how it compares with current practice. Radio buttons take similar arguments to selection boxes by commonly sharing choices and selected:

Incorporating a large number of widgets allows a variety of scenarios to be simulated without altering the input dataset, model code or R processing code for output.

Comparing the Shiny application and NONMEM^® simulations

Performing an independent check of the differential equation system and simulated population characteristics of a Shiny application is important. For applications such as Example 3, R simulation output can be evaluated with the simulation output provided by NONMEM^®13 (or similar) when given an equivalent population (i.e., input dataset) and dosing regimen.

Figure 5 compares the simulated output for a loading-bolus-bolus regimen (20-10-10 mg/kg) and a loading-continuous infusion regimen (16–24 mg/kg) from Example 3 with the output provided if the same population model was coded in NONMEM^® and an equivalent population of 1,000 individuals (i.e., post-natal age of 12 h and 1 kg weight) administered the same dosing regimens and their resultant concentration–time profiles were simulated at 15-min intervals. The series of plots indicate that the application produces median concentrations for R- and S-ibuprofen (blue and red solid lines, respectively), and 80% empirical confidence intervals for S-ibuprofen (red shaded ribbon) comparable to those simulated by NONMEM^®.

Simulating the population (1,000 individuals) and updating the output plot using a Shiny application was considerably faster than that achieved by the combined NONMEM^® simulation and R processing times (NONMEM^® + R method). Table 2 shows the mean computation times to simulate and process output for the two dosing regimens shown in Figure 5 by three different methods (i.e., NONMEM^® + R, uni-core Shiny, multi-core Shiny with compiled functions) when run four times each. Note that using NONMEM^® also required considerably more user time. Whereas running Example 3's application with one core and no compiled functions largely improved the run time when compared with NONMEM^® + R (≈ four- to sixfold reduction), an over 60-s wait for plot results was required for both dosing regimens. Compiling DES and simulate.conc functions and the addition of parallel processing for concentration–time profiles further improved the Shiny run times (≈ two- to threefold reduction).

Table 2. Shiny versus NONMEM^® + R processing: comparisona of run times

Methodb	IV loading-bolus-bolus (20–10-10 mg/kg) (seconds)	IV loading-continuous infusion (16–24 mg/kg) (seconds)
NONMEM® + R Processingc	≈ 385	≈ 366
Shiny: one core	≈ 103	≈ 60
Shiny: four cores, compiled functions	≈ 35	≈ 27

• ^a Based on application Example 3 with n = 1,000 and evaluations at 15-min intervals over 72 h.
• ^b Simulations were performed using a Dell® Power Edge R910 server with 4 x 10 core Xeon 2.26 GHz processors and 256 GB of RAM running Windows 8 Server SP1 64-bit.
• ^c NONMEM^® Version VII Level 3.014 using the Wings for NONMEM (Version 734) interface (http://wfn.sourceforge.net/), with Intel® Visual Fortran Composer XE 2013 Update 1, and ADVAN8 subroutine.

Whereas only two dosing scenarios have been demonstrated using Example 3, the user-interface of this Shiny application provided an easy-to-use platform to explore other scenarios by changing widget values for both the pharmacometricians who developed the application and the nonpharmacometricians who raised the research question.

SHARING SHINY APPLICATIONS

Shiny applications can easily be shared by passing on the ui.R and server.R scripts to a colleague with an installation of R or RStudio and the installed Shiny package. However, sharing is not limited to this mode. Shiny applications can also be hosted on Web servers making them accessible to those unfamiliar with, or do not have an installation of, R. There are three ways that RStudio offers to host a Shiny Web application: ShinyApps.io, Shiny Server, and Shiny Server Pro. Users will find this a slower option than running Shiny on their local machine—predominantly dependent on the speed and latency of their Internet connection. When a Shiny application is updating, the output is grayed out indicating to the user that processing is taking place.

CONCLUSIONS

The Shiny package of R is part of a larger movement amongst data scientists (e.g., Google charts, https://developers.google.com/chart) exploring interactive data visualization by means of Web-browsers. Like data science, the problem of communicating and sharing the information embodied in pharmacometric models is challenging. Representing pharmacometric models in R requires a degree of competency, but it is a skill that is increasingly common in the pharmacometric work force. We note that models coded in other software such as NONMEM^® are relatively easily converted to R, and believe that the Shiny package of R is an exciting development that allows pharmacometric models coded in R to be made accessible to a wider audience (e.g., drug development teams, clinicians). Elegant and flexible interfaces to models can be produced relatively quickly once the basics of Shiny are mastered. We also believe that pharmacometricians themselves may find Shiny applications a quick and informative way of simulating from their models. We hope that the information and examples provided here is sufficient to allow interested readers to start creating their own pharmacometric Shiny Web applications.