[R] Using apply() with a function involving ode()

[Adam Zeilinger]

Hello,

I am trying to get the output from the numerical simulation of a system 
of ordinary differential equations for a range of values for three 
parameters.  I am using the ode() function (deSolve package) to run the 
numerical simulation and apply() to run the simulation function for each 
set of parameter values.  I am having trouble getting the apply() 
function to work.

Here is an example.  Consider an epidemiology model of West Nile Virus:

wnv.model <- function(Time, State, Pars){
                   with(as.list(c(State, Pars)), {
                     dSB <- -(alphaB*betaB*SB*IM)/(SB + IB)
                     dIB <- (alphaB*betaB*SB*IM)/(SB + IB) - deltaB*IB
                     dSM <- bM*(SM + EM + IM) - (alphaM*betaB*IB*SM)/(SB 
+ IB) - dM*SM
                     dEM <- (alphaM*betaB*SM*IB)/(SB + IB) - kappaM*EM - 
dM*EM
                     dIM <- kappaM*EM - dM*IM
                     return(list(c(dSB, dIB, dSM, dEM, dIM)))
                   })
}

# I would like to run a numerical simulation of this model for each 
combination of values for the parameters alphaB, betaB, # and kappaM:

av <- seq(0, 1, by = 0.2) # vector of values for alphaB
bv <- seq(0, 1, by = 0.2) # vector of values for betaB
kv <- seq(0, 1, by = 0.2) # vector of values for kappaM

# Here is my function with ode() for the numerical simulation.  The 
function returns the last value of the simulation for the # "IB" state 
variable ("Infected Birds").

library(deSolve)
wnv.sim <- function(x){
          Pars <- c(alphaB = x[1], betaB = x[2], deltaB = 0.2, bM = 
0.03, dM = 0.03, alphaM = 0.69, kappaM = x[3])
          State <- c(SB = 100, IB = 0, SM = 1500, EM = 0, IM = 0.0001)
          Time <- seq(0, 60, by = 1)
          model.out <- as.data.frame(ode(func = wnv.incubation,
                        y = State,
                        parms = Pars,
                        times = Time))
          model.out[nrow(model.out),]$IB
}

# Finally, here is my apply() function:
dat <- apply(expand.grid(av, bv, kv), 1, wnv.sim)

However, the apply() function returns an error:
Error in eval(expr, envir, enclos) : object 'alphaB' not found

I am not sure what is wrong.  Any help would be much appreciated.

Sincerely,
Adam

-- 
Adam Zeilinger
Post Doctoral Scholar
Department of Entomology
University of California Riverside
www.linkedin.com/in/adamzeilinger