for drugs that operate in a narrow therapeutic range , it can be desirable or even necessary to use pharmacodynamic measures to adjust the dose to achieve an appropriate level of response . for example , anticoagulants such as warfarin act by reducing the ability of the blood to form clots . thus , pharmacodynamic measurement of blood clotting ability is used to adjust the dose within a specified range.12 
  establishing a dose adjustment algorithm requires quite a few decisions . what endpoint will be used to make dose adjustment decisions ? what is the target exposure range to maintain this endpoint within ? what should be the starting dose ? will monitoring and adaptive dosing continue indefinitely , or only during an initial period after treatment initiation ? simulations of these scenarios can be quite helpful in understanding the impact of these decisions . in this tutorial for illustration purposes , we have chosen a system where the pk is described by a twocompartment model , and the pd effect is modeled with an indirect response13 ( figure 
 1 ) . in order to maintain efficacy while avoiding adverse side effects , it has been determined that inhibition of the target should be in the range of 4060% . a twocompartment pharmacokinetic model with an indirect response pharmacodynamic model . first , we will demonstrate how to set up and run a simulation using rxode . the rxode package can be installed from github at https://github.com/hallowkm/rxode ( see the supplementary material for installation instructions ) . rxode is made available as open source under the gnu general public license version 2 or later . figure 
 2 gives an overview of the workflow for performing a simple simulation in rxode with the model described above . c2  =  centr / v2;c3  =  peri / v3;d / dt(depot )  =  ka*depot;d / dt(centr )  =  ka*depot  cl*c2 kout*(1c2/(ec50+c2))*eff ;  
  d / dt(depot )  =  ka*depot ; d / dt(centr )  =  ka*depot q*c2  +  q*c3 ; d / dt(peri )  =  q*c2  q*c3 ; d / dt(eff )  =  kin  kout*(1c2/(ec50+c2))*eff names of parameters in the vector must be a superset of parameters in the ode model , and the order of parameters within the vector is not important . the number of ics must equal exactly the number of odes in the model , and the order must be the same as the order in which the odes are listed in the model . c(ka  =  0.3 , cl  =  7 , v2  =  40 , q  =  10 , v3  =  300 , kin  =  0.2 , kout  =  0.2 , ec50  =  8)inits <  c(0 , 0 , 0 , 1 ) 
  # define parameters and initial conditions params <  c(ka  =  0.3 , cl  =  7 , v2  =  40 , q  =  10 , v3  =  300 , kin  =  0.2 , kout  =  0.2 , ec50  =  8) inits < c(0 , 0 , 0 , 1 ) the creation of an eventtable ( ) object provides an extremely efficient and flexible way to specify dosing and sampling . the generated eventtable object has functions that allow easy addition of dosing and sampling events . the add.dosing ( ) function allows specification of the dose amount , number of doses , dosing interval , compartment to dose into , rate ( if an infusion ) , and dosing start time . more complex dosing schedules can be simulated by applying the add.dosing ( ) function multiple times . the add.sampling ( ) function allows specification of the time points to be included in the simulation output . the add.dosing ( ) function provides an efficient method to specify a variety of dosing schedules calling the rxode ( ) function , the model is translated into c code , compiled , and dynamically loaded into the running r process . a simulation can then be performed by calling the r model object 's function run ( ) , with the specified parameter vector , initial conditions vector , and event table as inputs . rxode(model  =  ode , modname  =   mod1  ) # run simulationx <  mod1$run(params , ev , inits ) 
  mod1 < rxode(model  =  ode , modname  =   mod1  ) x <  mod1$run(params , ev , inits ) all state variables as well as other variables computed in the model are returned in the output matrix , at the times specified in the eventtable . thus , the simulation results are readily available for performing calculations and generating plots in r using any of the existing r packages ( lattice , ggplot , etc ) . the user can also choose to specify the absolute and/or relative tolerance , as well as the type of solver to be used : 
 x <  m1$run(theta , ev , inits , stiff  =  f , atol  =  1e8 , rtol  =  1e6 ) 
  x <  m1$run(theta , ev , inits , stiff  =  f , atol  =  1e8 , rtol  =  1e6 ) rxode uses the lsoda and a rungekutta integrators for stiff and nonstiff equations , respectively . the lsoda ( livermore solver for ordinary differential equations ) fortran package is an automatic method switching for stiff and nonstiff problems throughout the integration interval . for purely nonstiff systems , rxode uses dop853 , an explicit rungekutta method of order 8(5,3 ) . a matrix of parameter values can be generated , where each row represents one set of parameter values . r sampling functions such as rnorm and mvrnorm can be utilized to build this matrix , but decisions on the level(s ) of uncertainty to include are left to the discretion of the modeler , since this depends on the specific questions a given simulation is designed to address . the following code generates parameters for 100 subjects with correlated interindividual variability on cl and v2 . mvrnorm(n = nsub , rep(0,2 ) , sigma ) # sample from covariance matrixcl < 40*exp(mv[,2])params.all <  cbind(ka=0.3 , cl = cl , v2=v2 , q=10 , v3=300 , kin=0.2 , kout=0.2 , ec50=8 ) 
  nsub matrix(c(0.09,0.08,0.08,0.25),2,2 ) # iiv covariance matrix mv <  mvrnorm(n = nsub , rep(0,2 ) , sigma ) # sample from covariance matrix params.all <  cbind(ka=0.3 , cl = cl , v2=v2 , q=10 , v3=300 , kin=0.2 , kout=0.2 , ec50=8 ) once this parameter matrix is generated , each subproblem can be simulated by looping through the parameter matrix , using each row as an input for the simulation , and collecting the output of each simulation in an output matrix . res <  null # loop through each row of parameter values and simulationfor ( i in 1:nsub ) { params < params.all[i,]x <  mod1$solve(params , ev , inits = inits)#store results for effect compartmentres <  cbind(res , x [ , " eff"])}the same result can be achieved more efficiently with the following code : res <  apply(theta.all , 1 , function(theta ) mod$run(theta , ev , inits ) [ ,  eff  ] ) 
  # loop through each row of parameter values and simulation params <  params.all[i , ] x <  mod1$solve(params , ev , inits = inits ) # store results for effect compartment res <  cbind(res , x [ , " eff " ] ) the same result can be achieved more efficiently with the following code : res <  apply(theta.all , 1 , function(theta ) mod$run(theta , ev , inits ) [ ,  eff  ] ) simulation results can be then be directly analyzed and visualized using any of the statistical and graphics tools available within r. figure 
 3 shows the results of the above simulation when the drug is given qd for 2 days . simulation of effect for 2 days of qd dosing , correlated interindividual variability on cl and v2 . since simulations are conducted in r , rxode simulation results can easily be summarized and visualized graphically . to return to the problem of simulating adaptive dosing , there are many factors that must be considered in a designing an adaptive dosing scheme . utilizing rxode within the r environment , simulations of the impact of various factors can be quickly performed and evaluated . in order to simulate an adaptive dosing scenario , we will first simulate the following decision rule : the drug is to be dosed once daily , and trough pd effect levels will be measured 24 hours after each dose . if the measured pd effect is less than 40% , the dose will be doubled . if the measured pd effect is greater than 60% , the dose will be cut in half . if the pd effect is between 40 and 60% , no change will be made . the following r code is used to perform this simulation over 25 days . in brief , parameters governing the simulation are defined , including the decision rule limits and dose adjustments , number of days , starting dose , and sampling frequency . then treatment is simulated one day at a time , and after each day the simulated trough level at the end of that day is used to determine the dose level for the next simulated day . this is repeated for the number of days specified , and the results are stored in a matrix . results contained in this matrix can then be plotted using any available r plotting tools . effect.limits  =  c(0 , 0.4 , 0.6 , 9 ) # decision rule limitsdose.multipliers  =  c(0.5 , 1 , 2 ) # decision rule effectsndays <  25 ; unit.dose <  10000 ; start.dose <  1 ; sampling.frequency <  1 # sample every day # simulate each day . at the end of each day , test the effect level , and adjust the dose level according to the decision rulefor ( i in seq(1 , ndays , by  =  sampling.frequency ) ) { if ( i==1 ) { # initialize on first dayinits <  c(0 , 0 , 0 , 1)last.multiplier <  start.dosethis.multiplier <  1 } else { # use end of previous day as initial conditionsfor next day , compare trough effect with#decision rule limits and determine dose multiplier accordinglyinits <  x[dim(x ) , vars]wh <  cut(inits[eff  ] , effect.limits)this.multiplier <  dose.multipliers[wh]}this.multiplier <  this.multiplier*last.multiplier # adjust doselast.multiplier <  this.multiplier # store new dose#specify dosing and samplingev <  eventtable()ev$add.dosing(dose <  mod1$run(params , ev , inits ) # run simulation # compile outputstime.total < length(time))x <  cbind(x , time.total , doses);res <  rbind(res , x ) } 
  effect.limits  =  c(0 , 0.4 , 0.6 , 9 ) # decision rule limits dose.multipliers  =  c(0.5 , 1 , 2 ) # decision rule effects ndays <  25 ; unit.dose <  10000 ; start.dose <  1 ; sampling.frequency <  1 # sample every day # simulate each day . at the end of each day , test the effect level , and adjust the dose level according to the decision rule for ( i in seq(1 , ndays , by  =  sampling.frequency ) ) { if ( i==1 ) { # initialize on first day inits < c(0 , 0 , 0 , 1 ) last.multiplier <  start.dose } else { # use end of previous day as initial conditions for next day , compare trough effect with # decision rule limits and determine dose multiplier accordingly inits <  x[dim(x ) , vars ] wh < cut(inits[eff  ] , effect.limits ) this.multiplier <  dose.multipliers[wh ] this.multiplier <  this.multiplier*last.multiplier # adjust dose last.multiplier <  this.multiplier # store new dose # specify dosing and sampling ev$add.dosing(dose  =  this.multiplier*unit.dose , dosing.interval = 24 , nbr.doses  =  sampling.frequency ) ev$add.sampling(0:(24*sampling.frequency ) ) x <  mod1$run(params , ev , inits ) # run simulation time.total <  ev$get.eventtable()[,time]+(i1)*(24 ) doses <  rep(last.multiplier , length(time ) ) x < cbind(x , time.total , doses ) ; 
 figure 
 4 shows the results of this simulation . after 10 days , a steady state dose that is 25% of the starting dose is reached , which maintains the trough levels within the desired range . however , during the first 5 days the level of inhibition is well beyond the desired target . the red dashed line shows the dose multiplier over time ( actual dose is the dose multiplier times the starting dose ) . a stable dose that maintains the trough effect within the target range is reached after 10 days . for instance , figure 
 5 shows the simulation results if both trough and peak levels are controlled within the 4060% range ( i.e. , biomarker measurements are taken at tmax of 12 hours and at trough , and dose is adjusted daily ) . this full script for performing this simulation can be found in the appendix , but it requires only the addition of the following line to the script above : simulated adaptive dosing to maintaining both trough and peak effect within the target range of 4060% . the red dashed line shows the dose multiplier over time ( actual dose is the dose multiplier times the starting dose ) . a pattern emerges , indicating that dosing 12.5% of the starting dose , but then giving a double dose every 5 days , would achieve this goal . # if effect at 12 hours is less than 0.4 , cut dose in halfif ( x[13 ,  eff  ]  <  effect.limits ) { this.multiplier <  dose.multipliers } 
  # if effect at 12 hours is less than 0.4 , cut dose in half if ( x[13 ,  eff  ]  <  effect.limits ) { this.multiplier <  dose.multipliers in this case , a stable dose is not reached . instead , a more complex pattern emerges , suggesting that dosing 12.5% of the starting dose , but then giving a double dose every 5 days , could achieve the desired result . it may also be of interest to explore the impact of pk and pd variability within the population on the resulting dose trajectory and final dose . we can simulate variability , as described above , by specifying a matrix of parameter sets , rather than a fixed set of parameters , and looping through this parameter matrix , performing the adaptive simulation as described above for each set of parameters . this is easily done by adapting the simulation of a single subject to contain an outer forloop that cycles through the parameter matrix . again , we assume correlated interindividual variability on cl and v2 , and simulate 100 subjects for 25 days . trough levels for all subjects are controlled within the specified range by the end of the time period . 62% of subjects end at 25% of the starting dose , 35% require only 12.5% , and a few subjects ( 2.5% ) require a higher dose ( 50% of starting dose ) . it took at least 20 days to reach a final dose in all subjects ( compared to 10 days for the typical subject in the previous simulation ) . simulation of adaptive dosing regimen in 1,000 subjects , with correlated interindividual variability on cl and v2 . ( a ) the adaptive dosing regimen successfully controls effect levels within the desired range after 5 days . ( b ) the dose trajectory and final dose varies among the population , and most subjects end on either 25% or 12.5% of the starting dose . the simulations performed thus far are static , in that fixed parameter inputs are defined , simulations are performed , and plots are generated . often , when the results of such simulations are reviewed , further questions arise  what happens if the range limits for our decision rule are relaxed or tightened ? what happens if we change the starting dose ? what happens if we sample every 2nd or 3rd day , rather than every day ? these new questions usually require the modeler to go off and perform new simulations , and the team must reconvene days or weeks later to discuss new results . this introduces significant lag time , especially since several iterations may be required to reach a final set of simulations . although it may be possible to compile a massive document containing all the permutations of scenarios of interest a priori , sifting through this document with team members in real time can be a daunting task , and is not advisable , as it can lead to confusion and/or loss of attention of key stakeholders . a better alternative is to encase the simulation procedure in an interactive application in which a wide range of scenarios can be evaluated in real time , as they arise . a major advantage of performing simulations within the r environment is the ability to take advantage of r 's shiny package for developing interactive web applications . these web applications can be made available to the team online , and used in meetings to explore simulation scenarios in real time . this is extremely advantageous in facilitating modeler  team interactions , because it can greatly reduce the number of iterations , separate meetings , and associated lag time . shiny apps are easy to write , requiring no web development skills and very limited programming skills . genshinyapp.template(appdir  =   shinyexample  , verbose  =  true)library(shiny ) # load the shiny packagerunapp(shinyexample  ) # run the example app 
  genshinyapp.template(appdir  =   shinyexample  , verbose  =  true ) library(shiny ) # load the shiny package runapp(shinyexample  ) # run the example app the function genshinyapp.template ( ) generates a folder that contains the r shiny ui.r and server.r files for the template app , as well as an rda file containing the saved example model , parameters , and initial conditions . to incorporate a different r model , the user needs can save the new model to an rda file and load this new file from the server.ui file . # save the model , parameters , init values , etc . in the file rx_shiny_data.rda to be loaded by the server.rsave(mod1 , params , inits , stiff  =  true , atol  =  1e8 , rtol  =  1e6 , file  =   rx_shiny_data.rda  ) 
  # save the model , parameters , init values , etc . in the file rx_shiny_data.rda to be loaded by the server.r save(mod1 , params , inits , stiff  =  true , atol  =  1e8 , rtol  =  1e6 , file  =   rx_shiny_data.rda  ) to tailor the app , the only two files that need to be altered by the user are the ui.r and server.r files . the ui.r file can be edited to add widgets for obtaining user input and to control the layout of the user interface . the server.r file can be edited to control how inputs are used and how outputs are displayed . there is excellent documentation on developing shiny app interfaces available online8 as well as a previous tutorial on this topic in this journal.14 
 
 figure 
 7 shows the shiny app interface for exploring different factors in an adaptive dosing regimen , including the starting dose , the number of days of simulation , the frequency at which the dose is adjusted , the lower and upper limits of the range for the decision rule , and whether trough , peak , or both levels are controlled by the decision rule . it displays the calculated % time above and below the decision rule limits , and it plots the effect and dose as a function of time . as the user moves the sliders or checks the check boxes , the outputs are adjusted interactively . the ui.r , server.r , rxode.run.r files and the compiled c model are available in the supplementary material . rxode provides a function for generating interactive shiny apps , which can then be customized . this app allows users to vary the dose , number of days , sampling frequency , and decision rules , and to view simulation results in real time . thus , starting within a baseline scenario , many different scenarios and combinations of scenarios can be quickly explored . while the modeler will likely facilitate evaluation of scenarios , these apps require no technical expertise or experience with r , and thus the apps are accessible to nontechnical partners . also , importantly , shiny apps can be easily deployed online and do not require local installation or r. this is in contrast to other software such as berkeley madonna . while berkeley madonna can facilitate interactive model exploration , it also requires purchase of a software license , local installation of the software , and some user familiarity with the software and modeling . rxode is an r package that provides tools for the efficient simulation of complex dosing regimens via pk / pd models described by odes . it provides great flexibility and speed in performing simulations with variability and uncertainty . as part of the r environment , rxode outputs can be combined with a multitude of r facilities to create advanced static and interactive visualization displays for effective communications with clinical team members and other consumers . furthermore , unlike many other simulation tools , simulation and preparation of graphics can be conducted completely within a single , freely available , and opensource software . no licenses are required , and it does not require linking with any external software . although currently focused on efficient simulations , rxode can also be used for parameter estimation through the many existing statistical estimation algorithms in r , including nonlinear mixed effects models,14 stochastic approximation expectationmaximization ( saem),16 and bayesian methods using gibbs sampling , e.g. , jags.17 future work includes developing functionality to aid users in linking rxode models with these estimation algorithms in a more efficient manner . supplementary material is linked to the online version of the article at wileyonlinelibrary.com.cpt 
  click here for additional data file .
this tutorial presents the application of an r package , rxode , that facilitates quick , efficient simulations of ordinary differential equation models completely within r. its application is illustrated through simulation of design decision effects on an adaptive dosing regimen . 
 the package provides an efficient , versatile way to specify dosing scenarios and to perform simulation with variability with minimal custom coding . 
 models can be directly translated to rshiny applications to facilitate interactive , realtime evaluation / iteration on simulation scenarios .