Ron Ron - 3 months ago 12
R Question

How to structure ODEs in R based on multiple groups

I am trying to simulate cell uptake in R, having ported a model from Berkeley Madonna. The model is comprised of several constants and differential equations to calculate amounts and concentrations. A portion of the code is listed:

library(deSolve)

fb = 0.0510
Km = 23.5
Pdif = 0.429
Vmax = 270
Vol_cell = 9.33
Vol_media = 150
S = 10 #concentration of dosing media

yini = c(Amt_media=(S*Vol_media)-(S*fb*Vol_cell),
Amt_cell=S*fb*Vol_cell,
Amt_total=S*Vol_media,
Con_media=S-(S*fb),
Con_cell=S*fb)

Uptake = function(t, y, p){
dy1 = (- (Pdif * y[1]) + (Pdif * y[2]) - ((Vmax * y[4])/(Km + y[4])))
dy2 = (+ (Pdif * y[1]) - (Pdif * y[2]) + ((Vmax * y[4])/(Km + y[4])))
dy3 = dy1 + dy2
dy4 = dy1 / Vol_media
dy5 = dy2 / Vol_cell
list(c(dy1, dy2, dy3, dy4, dy5))}

times1 = seq(from=0, to=15, by=0.01)
out1 = ode(y=yini, times=times1, func=Uptake, parms=NULL, method="rk4")


The rest of the code is for output to dataframes and plotting. My question then is how to have the code structured to use "S" as a list of several concentrations such that each concentration can be applied to the differential equations (essentially giving me an out1 for S1, out2 for S2, etc, that can then be passed onto a dataframe)? In Berkeley Madonna this was achieved by writing over 35 differential equations, though I'd like to use a simplified approach in R if possible.

Answer

The only part where S is used is in the initialization of the yini values. Basically we just need to move that part and the part that runs ode with those values into a new function. Then you can call that function for what ever values you want. For example

#set up
library(deSolve)

fb <- 0.0510
Km <- 23.5
Pdif <- 0.429
Vmax <- 270
Vol_cell <- 9.33
Vol_media <- 150

Uptake <- function(t, y, p){
   dy1 = (- (Pdif * y[1]) + (Pdif * y[2]) - ((Vmax * y[4])/(Km + y[4]))) 
   dy2 = (+ (Pdif * y[1]) - (Pdif * y[2]) + ((Vmax * y[4])/(Km + y[4])))
   dy3 = dy1 + dy2
   dy4 = dy1 / Vol_media
   dy5 = dy2 / Vol_cell
   list(c(dy1, dy2, dy3, dy4, dy5))}

times1 <- seq(from=0, to=15, by=0.01)

# function with S as a parameter
runConc <- function(S) {
  yini <- c(Amt_media=(S*Vol_media)-(S*fb*Vol_cell), 
    Amt_cell=S*fb*Vol_cell, 
    Amt_total=S*Vol_media, 
    Con_media=S-(S*fb), 
    Con_cell=S*fb)

    ode(y=yini, times=times1, func=Uptake, parms=NULL, method="rk4")
}

#run for concentrations 10,20,30
out <- lapply(c(10,20,30), runConc)

This will result in a list object with the results for each concentration. So out[[1]] is the result for S=10, out[[2]] is S=20, etc. We can see the first few lines of each of the results with

lapply(out, head, 3)

# [[1]]
#      time Amt_media Amt_cell Amt_total Con_media Con_cell
# [1,] 0.00  1495.242  4.75830      1500  9.490000 0.510000
# [2,] 0.01  1488.103 11.89710      1500  9.442408 1.275145
# [3,] 0.02  1481.028 18.97216      1500  9.395241 2.033457
# 
# [[2]]
#      time Amt_media Amt_cell Amt_total Con_media Con_cell
# [1,] 0.00  2990.483  9.51660      3000  18.98000 1.020000
# [2,] 0.01  2976.550 23.44980      3000  18.88711 2.513377
# [3,] 0.02  2962.739 37.26072      3000  18.79504 3.993646
# 
# [[3]]
#      time Amt_media Amt_cell Amt_total Con_media Con_cell
# [1,] 0.00  4485.725 14.27490      4500  28.47000  1.53000
# [2,] 0.01  4465.153 34.84653      4500  28.33286  3.73489
# [3,] 0.02  4444.761 55.23920      4500  28.19690  5.92060