[R-sig-dyn-mod] Using Desolve to simulate drug concentrations
Meenakshi Srinivasan
m.srinivasan29193 at gmail.com
Mon Mar 26 20:49:53 CEST 2018
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,
[[alternative HTML version deleted]]
More information about the R-sig-dynamic-models
mailing list