[R-sig-dyn-mod] Using Desolve to simulate drug concentrations

Dear group,

*I'm new to R and deSolve. Using an example I found in a publication, I
tried to develop a model that can simulate drug concentrations for various
scenarios (based on body weight, albumin, sex and other covariates).*

*The running model code is as below-*

```
#Example R script for simulating a population and concentrations as
described by:
#2-compartment model
#Dosing regimen:
# IV infusion (at 0 hours, duration = 10 min, 0.16 hr)
#
#Population parameter variability on and CL, V1, V2, Q
#Variance-covariance matrix for  and CL, V1, V2
#Covariate effects: CRCL on CL, ISM on CL
# ALB on V1
#         WT allometrically scaled to all V1, V2 and Q
#Proportional and additive residual error model

#---------------------------------------------------------------------------------------------------------------------------------------------
#Remove all current objects in the workspace
rm(list=ls(all=TRUE))

#Load package libraries
library(deSolve) #Differential equation solver
library(ggplot2) #Plotting
library(plyr)   #Split and rearrange data, ddply function
library(grid)   #Plotting
library(MASS)   #mvrnorm function
library(MBESS) #cor2cov function
library(compiler) #Compile repeatedly-called functions


#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 1.  Create a parameter dataframe with ID and parameter values for
each individual

#Define a population
n <- 1000
ID <- seq(from = 1, to = n, by = 1)
WT <- rlnorm(n, meanlog = log(58), sdlog = 0.205) #Log normal distribution
for body weight (kg)
SEX <- rbinom(n, size = 1, prob = 0.65)             #Binomial distribution
for sex  0 is Female, 1 is male  ; 65% males
CRCL<- rlnorm(n, meanlog = 4.28923 , sdlog = 0.2685873)
ALB<- rlnorm(n, meanlog = 1.28567, sdlog =0.1346275 )
BSA<-rlnorm(n, meanlog = 0.4507393, sdlog = 0.1184323)
#Check demographics
hist(WT)
hist(CRCL)
table(SEX)


#Define parameter values
#Thetas
CLPOP <- 3.14687  #Clearance L/hour
V1POP <- 5.46589 #Central volume L
V2POP <- 6.75549   #Peripheral volume L
QPOP <-  5.7451    #Inter-compartmental clearance L/h

COV1 <- 0.520245    #Covariate effect of CRCL on clearance (power model)
COV2 <- 1.22484    #Covariate effect of  gender on CL
COV3  <- -1.51075   #Covariate effect of ALB on V1

#Omegas (as SD)
ETA1SD <- 0.520860826
ETA2SD <- 0.507793265
ETA3SD <- 1.222178383
ETA4SD <- 0.881369956

#Specify a correlation matrix for ETA's
R12 <- 0.387845089 #Correlation coefficient for CL-V1
R13 <- 0           #Correlation coefficient for CL-V2
R23 <- 0.462686988 #Correlation coefficient for V1-V2
R14 <- 0            #Correlation coefficient for CL-Q
R24 <- 0            #Correlation coefficient for V1-Q
R34 <- 0            #Correlation coefficient for V2-Q

#Epsilons (as SD)
EPS1SD <- 0.287751281  #Proportional residual error
EPS2SD <- 0.643975931  #Additive residual error


#Calculate ETA values for each subject
CORR <- matrix(c(1,R12,R13,
R14,R12,1,R23,R24,R13,R23,1,R34,R14,R24,R34,1),4,4)



#Specify the between subject variability for CL, V1, V2
SDVAL <- c(ETA1SD,ETA2SD,ETA3SD,ETA4SD)

#Use this function to turn CORR and SDVAL into a covariance matrix
OMEGA <- cor2cov(CORR,SDVAL)

#Now use multivariate rnorm to turn the covariance matrix into ETA values
ETAmat <- mvrnorm(n = n, mu = c(0,0,0,0), OMEGA)
ETA1 <- ETAmat[,1]
ETA2 <- ETAmat[,2]
ETA3 <- ETAmat[,3]
ETA4 <- ETAmat[,4]

#Define individual parameter values
CL <- CLPOP*exp(ETA1)*(CRCL/75)^COV1*COV2^SEX
V1 <- V1POP*exp(ETA2)*(WT/58)*(ALB/3.6)^COV3
V2 <- V2POP*exp(ETA3)*(WT/58)
Q  <- QPOP*exp(ETA4)*(WT/58)^0.75


#Calculate rate-constants for differential equation solver
K12 <- Q/V1
K21 <- Q/V2
K10 <- CL/V1


#Collect the individual parameter values in a parameter dataframe
par.data <- data.frame(ID,CL,V1,V2,Q,K12,K21,K10,WT,SEX,CRCL,ALB,BSA)
head(par.data)


#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 2.  Make a time sequence and specify the dose information for a
system of differential equations
#There are a number of ways that doses can be coded for the deSolve package
- see help for deSolve

#Creating continuous infusion (for 24 hours) - this use the approxfun
function to make a "forcing function" for infusion rate in the differential
equations
CDOSE <- 700 #mg
CTinf <-  0.16
CRATE <- CDOSE/CTinf   #mg/h
      #hours

CTIMEinf <- c(0,CTinf,24) #100 - make sure the function works long after
the 24 hours
CRATEinf <- c(CRATE,0,0)

#Define an interpolation function that returns rate when given time -
"const"
#i.e. at TIME = 0, RATE = CRATE. At TIME = 0.16, RATE = 0 (infusion stopped)
Cstep.doseinf <- approxfun(CTIMEinf, CRATEinf, method = "const")
#Testing
Cstep.doseinf(0)  #Returns the infusion rate for time = 0 etc.
Cstep.doseinf(0.1111)
Cstep.doseinf(0.15)
Cstep.doseinf(2.1)
Cstep.doseinf(5)

#Make a TIME sequence (hours)
TIME <- seq(from = 0, to = 24, by = 0.1)

#The time sequence must include all "event" times for deSolve so add them
here and sort
#Do not repeat a time so use "unique" as well


#Function containing differential equations for amount in each compartment
DES <- function(T, A, THETA) {

  RateC <- Cstep.doseinf(T) #Infusion rate

  K12 <- THETA[1]
  K21 <- THETA[2]
  K10 <- THETA[3]

  dA <- vector(length = 2)
  dA[1] =  RateC -K12*A[1] +K21*A[2] -K10*A[1] #Central compartment
  dA[2] =  K12*A[1] -K21*A[2]         #Peripheral compartment

  list(dA)
}

#Compile DES function - it's called by lsoda for each individual in the
dataset
DES.cmpf <- cmpfun(DES)


#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 3.  Run the differential equation solver for each patient in par.data
#Function for simulating concentrations for the ith patient
simulate.conc <- function(par.data) {

  #List of parameters from input for the differential equation solver
  THETAlist <- c("K12" = par.data$K12,"K21 "= par.data$K21,"K10" =
par.data$K10)

  #Set initial compartment conditions
  A_0 <- c(A1 = 0, A2 = 0)

  #Run differential equation solver for simulated variability data
  var.data <- lsoda(A_0, TIME, DES.cmpf, THETAlist)
  var.data <- as.data.frame(var.data)
}

#Compile simulate.conc function - it's called by ddply for each individual
in the dataset
simulate.conc.cmpf <- cmpfun(simulate.conc)

#Apply simulate.conc.cmpf to each individual in par.data
#Whilst maintaining their individual values for V1, SEX and WT for later
calculations
sim.data <- ddply(par.data, .(ID,V1,SEX,WT,ALB,CRCL,BSA),
simulate.conc.cmpf)

#Calculate individual concentrations in the central compartment
sim.data$IPRE <- sim.data$A1/sim.data$V1


#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 4.  Add residual error

#Use random number generator to simulate residuals from a normal
distribution
nobs <- n*length(TIME) #number of observations = number of subjects *
number of time points per subject
EPS1 <- rnorm(nobs, mean = 0, sd = EPS1SD) #Proportional residual error
EPS2 <- rnorm(nobs, mean = 0, sd = EPS2SD) #Additive residual error
sim.data$DV <- sim.data$IPRE*(1 + EPS1) + EPS2


#---------------------------------------------------------------------------------------------------------------------------------------------
#Step 5.  Draw some plots of the simulated data

#Use custom ggplot2 theme
theme_bw2 <- theme_set(theme_bw(base_size = 16))
theme_bw2 <- theme_update(plot.margin = unit(c(1.1,1.1,3,1.1), "lines"),
                          axis.title.x=element_text(size = 16, vjust = 0),
                          axis.title.y=element_text(size = 16, vjust = 0,
angle = 90),
                          strip.text.x=element_text(size = 14),
                          strip.text.y=element_text(size = 14, angle = 90))

#Factor covariate values (for example, SEX)
sim.data$SEXf <- as.factor(sim.data$SEX)
levels(sim.data$SEXf) <- c("Female","Male")

#Function for calculating 5th and 95th percentiles for plotting
concentrations
CIlow <- function(x) quantile(x, probs = 0.05)
CIhi <- function(x) quantile(x, probs = 0.95)

#Generate a plot of the sim.data
plotobj <- NULL
plotobj <- ggplot(data = sim.data)
plotobj <- plotobj + stat_summary(aes(x = time, y = DV), fun.ymin = CIlow,
fun.ymax = CIhi, geom = "ribbon", fill = "blue", alpha = 0.2)
plotobj <- plotobj + stat_summary(aes(x = time, y = IPRE), fun.ymin =
CIlow, fun.ymax = CIhi, geom = "ribbon", fill = "red", alpha = 0.3)
plotobj <- plotobj + stat_summary(aes(x = time, y = IPRE), fun.y = median,
geom = "line", size = 1, colour = "red")
plotobj <- plotobj + scale_y_continuous("Concentration (mg/L) \n")
plotobj <- plotobj + scale_x_continuous("\nTime (hours)")
print(plotobj)
```

Now, I tried integrating it with Shiny with a UI to change parameters-
which changes simulations.

The codes are as below-

```
#ui.R

fixedPage(

#Application Title
fixedRow(
column(10,
h2("Drug Simulation tool"),
h6("Drug simulation tool"),
    offset = 1, align = "center")
), #Brackets closing "fixedRow"
hr(), #Add a break with a horizontal line

####################################################################################
#Sidebar panel with widgets
sidebarLayout(
sidebarPanel(
      #Slider input for BSA
sliderInput("BSA",
            "BSA (m2):",
            min = 1.08,
            max = 2.16,
            value = 1.08,
            step = 0.01),
textOutput("DOSE"),   #BSA Based dosing. This will show the output of the
dose to be
                      #administered on usual regimen 500 mg/ m2
#Slider input for Dose
sliderInput("CDOSE",
            "Dose administered as short infusion (mg):",
            min = 500,
            max = 1000,
            value = 1,
            step = 1),
#Numeric output for infusion rate
textOutput("RATE"),

#Radio button to choose sex of patient
radioButtons("ISM",
             "Sex:",
             choices = list("Female" = 0,
                            "Male" = 1),
             selected = 1),
#Slider input for weight
sliderInput("Weight",
"Weight (kg):",
min = 0,
max = 150,
value = 12,
step = 1),
#Slider input for CRCL
sliderInput("CRCL",
            "CRCL (ml/min):",
            min = 0,
            max = 150,
            value = 12,
            step = 5),
#Slider input for Albumin
sliderInput("ALB",
            "Albumin (g/dL):",
            min = 2.4,
            max = 4.8,
            value = 2.5,
            step = 0.1),
#Slider input for number of individuals to carry out simulation on
sliderInput("n",
"Number of Individuals:",
min = 10,
max = 1000,
value = 10,
step = 10),
br(),
#Button to initiate simulation
submitButton("Simulate"),
align = "left"), #Brackets closing "sidebarPanel"
mainPanel(
#Plot output for concentration-time profile
plotOutput("plotCONC", height = 650, width = 750),
align = "center") #Brackets closing "mainPanel"
) #Brackets closing "sidebarLayout"

) #Brackets closing "fixedPage"

````

*And server.R*

```
#Load package libraries
library(shiny)    #for application
library(ggplot2)  #plotting
library(deSolve)  #differential equation solver
library(plyr)     #Split and rearrange data, ddply function
library(reshape2)
library(compiler)
library(doParallel)
library(grid)   #Plotting
library(MASS)   #mvrnorm function
library(MBESS) #cor2cov function

#Code for functions and variables which are not reactive (not dependent on
"input$X")
#Setting up cores to run parallel processes, thus increasing speed
#Set up a cluster of cores to run the application over
# Creates a set of copies of R running in parallel and communicating over
sockets.
#package- parallel
cl <- makePSOCKcluster(detectCores() - 1)
#detectCores() searches for the number of cores that the local machine has
#Contents with makePSOCKcluster brackets can be changed to a whole number
if you
#want to assign an exact number
#List packages that are required to be sent to each core for the parallel
process
#The foreach package always needs to be included
#This example uses the .parallel argument in ddply which calls a function
that uses
#lsoda from the deSolve package
#ClusterevalQ- These functions provide several ways to parallelize
computations using a cluster.
clusterEvalQ(cl, list(library(deSolve), library(foreach)))
#Registers the parallel backend with the foreach package (automatically
loaded when
#doParallel is loaded)
registerDoParallel(cl)

#ggplot2 theme for plotting
theme_custom <- theme_set(theme_grey(18))
#TIME range - times where concentrations will be calculated
TIME <- seq(from = 0, to = 24, by = 0.25)
#----------------------------------------------------------------------------------------
#Define user-input dependent functions for output
shinyServer(function(input, output){
  #Reactive expression to generate a reactive data frame
  #This is called whenever the input changes
  all.data <- reactive({


#----------------------------------------------------------------------------------------
    #Step 1.  Create a parameter dataframe with ID and parameter values for
each individual

    #Define a population
    n <- input$n
    par.data <- seq(from = 1, to = n, by = 1)
    par.data <- data.frame(par.data)
    names(par.data)[1] <- "ID"

    #Input for weight from UI
    WT <- input$Weight

    #Input for sex from UI
    SEX <- input$ISM

    #Input for CRCL form UI
    CRCL<- input$CRCL

    #Input of Albumin level from UI
    ALB<- input$ALB

    #Input for BSA from UI
    BSA<- input$BSA

    #Define parameter values
    #Thetas
    CLPOP <- 3.14687  #Clearance L/hour
    V1POP <- 5.46589 #Central volume L
    V2POP <- 6.75549   #Peripheral volume L
    QPOP <-  5.7451    #Inter-compartmental clearance L/h

    COV1 <- 0.520245    #Covariate effect of CRCL on clearance (power model)
    COV2 <- 1.22484    #Covariate effect of  gender on CL
    COV3  <- -1.51075   #Covariate effect of ALB on V1

    #Epsilons (as SD)
    EPS1SD <- 0.287751281  #Proportional residual error
    EPS2SD <- 0.643975931  #Additive residual error

    OMEGA <-
matrix(c(0.271,0.102,0,0,0.102,0.257,0.291,0,0,0.291,1.5,0,0,0,0,0.774),4,4)

    #Now use multivariate rnorm to turn the covariance matrix into ETA
values
    #mvrnorm- package MASS- Produces one or more samples from the specified
    #multivariate normal distribution
    #n-number of samples required
    #mu- vector giving means of the variables
    #Sigma-omega-a positive-definite symmetric matrix
    #specifying the covariance matrix of the variables.

    ETAmat <- mvrnorm(n = n, mu = c(0,0,0,0), OMEGA)
    ETA1 <- ETAmat[,1]
    ETA2 <- ETAmat[,2]
    ETA3 <- ETAmat[,3]
    ETA4 <- ETAmat[,4]

    #Define individual parameter values
    #SEX=male
    #Clearance
     if(SEX==1){
       par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1*COV2
     }
    #SEX=Female
    #Clearance
    if(SEX==0){
      par.data$CL<- CLPOP*exp(ETA1)*(CRCL/75)^COV1
    }
    #Central volume
    par.data$V1 <- V1POP*exp(ETA2)*(WT/58)*(ALB/3.6)^COV3

    #Peripheral volume
    par.data$V2 <- V2POP*exp(ETA3)*(WT/58)

    #Inter-compartment clearance
    par.data$Q  <- QPOP*exp(ETA4)*(WT/58)^0.75

    #Calculate rate-constants for differential equation solver
    par.data$K12 <- par.data$Q/par.data$V1
    par.data$K21 <- par.data$Q/par.data$V2
    par.data$K10 <- par.data$CL/par.data$V1

    #  Add residual error

    #Use random number generator to simulate residuals from a normal
distribution
    nobs <- n*length(TIME) #number of observations = number of subjects *
number of time points per subject
    EPS1 <- rnorm(nobs, mean = 0, sd = EPS1SD) #Proportional residual error
    EPS2 <- rnorm(nobs, mean = 0, sd = EPS2SD) #Additive residual error


#------------------------------------------------------------------------------------------
    #Step 2.  Make a time sequence and specify the dose information for a
system of differential equations
    #There are a number of ways that doses can be coded for the deSolve
package - see help for deSolve

    #Creating continuous infusion (for 24 hours) - this use the approxfun
function to make a "forcing function" for infusion rate in the differential
equations
    CDOSE <- input$CDOSE #mg
    CTinf <-  0.16        #h
    CRATE <- CDOSE/CTinf   #mg/h


    CTIMEinf <- c(0,CTinf,24) #100 - make sure the function works long
after the 24 hours
    CRATEinf <- c(CRATE,0,0)

    #Define an interpolation function that returns rate when given time -
"const"
    #i.e. at TIME = 0, RATE = CRATE. At TIME = 0.16, RATE = 0 (infusion
stopped)
    #approxfun- package stats, Return a list of points which linearly
interpolate given data points, or a function
    #performing the linear (or constant) interpolation.
    Cstep.doseinf <- approxfun(CTIMEinf, CRATEinf, method = "const")

    #The time sequence must include all "event" times for deSolve so add
them here and sort
    #Do not repeat a time so use "unique" as well
    #specified above

    #Function containing differential equations for amount in each
compartment
     DES <- function(T, A, THETA) {

       RateC <- Cstep.doseinf(T) #Infusion rate

        K12 <- THETA[1]
        K21 <- THETA[2]
        K10 <- THETA[3]
       dA <- vector(length = 2)
       dA[1] =  RateC -K12*A[1] +K21*A[2] -K10*A[1]   #Central compartment
       dA[2] =  K12*A[1] -K21*A[2]         #Peripheral compartment

       list(dA)
     }

    #Compile DES function - it's called by lsoda for each individual in the
dataset
    #package- compiler-provide an interface to a byte code compiler for R
    DES.cmpf <- cmpfun(DES)



#---------------------------------------------------------------------------------------------------------------------------------------------
    #Step 3.  Run the differential equation solver for each patient in
par.data
    #Function for simulating concentrations for the ith patient
    simulate.conc <- function(par.data) {

      #List of parameters from input for the differential equation solver
      THETAlist <- c("K12" = par.data$K12,"K21"= par.data$K21,"K10" =
par.data$K10)

      #Set initial compartment conditions
      A_0 <- c(A1 = 0, A2 = 0)

      #Run differential equation solver for simulated variability data
      #lsoda- package desolve
      #Solving initial value problems for stiff or non-stiff systems of
      #first-order ordinary differential equations (ODEs).
      #The R function lsoda provides an interface to the FORTRAN ODE solver

      sim.data <- lsoda(A_0, TIME, DES.cmpf, THETAlist)
      sim.data <- as.data.frame(sim.data)
    }

    #Compile simulate.conc function - it's called by ddply for each
individual in the dataset
    simulate.conc.cmpf <- cmpfun(simulate.conc)

    #Apply simulate.conc.cmpf to each individual in par.data
    #Whilst maintaining their individual values for V1, SEX and WT for
later calculations

    sim.data <- ddply(par.data, .(ID,V1,V2), simulate.conc.cmpf,
.parallel=TRUE)

    #Calculate individual concentrations in the central compartment
    sim.data$IPRE <- sim.data$A1/sim.data$V1


#---------------------------------------------------------------------------------------------------------------------------------------------

    #simulating observations (Y)
    sim.data$DV <- sim.data$IPRE*(1 + EPS1) + EPS2
})

#---------------------------------------------------------------------------------------------------------------------------------------------
  #   #Step 5.  Draw a plot of the simulated data
   CIlow <- function(x) quantile(x, probs = 0.05)
   CIhi <- function(x) quantile(x, probs = 0.95)

   output$plotCONC <- renderPlot({
   plotobj <- ggplot(data=all.data())
   plotobj <- plotobj + geom_line(aes(x = TIME, y = DV), colour = "red",
size = 1)
   #plotobj <- plotobj + geom_ribbon(aes(x = time, ymin = CIlow, ymax =
CIhi), fill = "red", alpha = 0.3)
   plotobj <- plotobj + scale_y_continuous("Concentration (microg/mL) \n",
lim=c(0,600))
   plotobj <- plotobj + scale_x_continuous("\nTime (hours)",
breaks=c(0,8,16,24))
   print(plotobj)

    })

    output$DOSE <- renderText({
      paste("Dose to be administered =", signif(input$BSA*500, digits = 3)
,"mg")

    }) #Brackets closing "renderText" function

    output$RATE <- renderText({
      paste("Infusion rate =", signif(input$CDOSE*6, digits = 3) ,"mg/hr")

    }) #Brackets closing "renderText" function

  }) #bracket closing shiny server

```

*I'm getting the error which is as follows-*

```
Listening on http://127.0.0.1:6899
Warning: <anonymous>: ... may be used in an incorrect context: ?.fun(piece,
...)?

Warning: <anonymous>: ... may be used in an incorrect context: ?.fun(piece,
...)?

Warning: Error in : ggplot2 doesn't know how to deal with data of class
numeric
Stack trace (innermost first):
    105: fortify.default
    104: fortify
    103: structure
    102: ggplot.data.frame
    101: ggplot.default
    100: ggplot
     99: renderPlot [C:\Users\ms1225\Downloads\Shiny_code_2/server.R#222]
     89: <reactive:plotObj>
     78: plotObj
     77: origRenderFunc
     76: output$plotCONC
      1: runApp
```

*Can you please help me with what I might have missed?*

*Thank you so much for your help. Will be happy to clarify any doubts
regarding the code. *
*The server is essentially the final model code which is modified.*

Thanks!

Sincerely,