Applying a set of ordinary differential equations to each grid cell

Question

I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame landscape (see below) is updated at each time step. Here is an example at t = 0:

landscape <- data.frame(polygon_ID = seq(1, 10, by = 1), 
                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE), 
                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))

To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:

solve_sir_model <- function (times, parameters) { 

  sir_model <- function (times, states, parameters) { 

    with(as.list(c(states, parameters)), { 

      dSdt <- -beta*S*I 
      dIdt <- beta*S*I-gamma*I 
      dRdt <- gamma*I 
      dNdt <- dSdt + dIdt + dRdt 

      return(list(c(dSdt, dIdt, dRdt, dNdt))) 

    }) 
  } 

  states <- c(S = 99, I = 1, R = 0, N = 100) 
  return(ode(y = states, times = times, func = sir_model, parms = parameters)) 
} 

require(deSolve) 
output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))

At each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame landscape? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?

EDIT

It seems that the method iteration in the function ode can be useful in my case:

Method "iteration" is special in that here the function func should return the new value of the state variables rather than the rate of change. This can be used for individual based models, for difference equations, or in those cases where the integration is performed within func).

I have tested the method but I don't understand why it doesn't work with a single time step:

solve_sir_model <- function (times, parameters) { 

  sir_model <- function (times, states, parameters) { 

    with(as.list(c(states, parameters)), { 

      dSdt <- -beta*S*I 
      dIdt <- beta*S*I-gamma*I 
      dRdt <- gamma*I 
      dNdt <- dSdt + dIdt + dRdt 

      return(list(c(dSdt, dIdt, dRdt, dNdt))) 

    }) 
  } 

  states <- c(S = 99, I = 1, R = 0, N = 100) 
  return(ode(y = states, times = times, func = sir_model, parms = parameters, method = "iteration")) 
} 

require(deSolve) 
output <- as.data.frame(solve_sir_model(times = 1, parameters = c(beta = 400, gamma = 28)))

 Error in iteration(y, times, func, parms, ...) : 
   times should be equally spaced In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

You should check out the ReacTran package for models with a spatial component like this, as they have helper functions to set up your grid, etc. It sounds like you're talking about temporally variable parameters, which can also be accommodated by ReacTran, but beware that sharp (i.e. discrete) shifts in parameters cause problems with numerical solvers. I'd recommend approximating the parameters over time using something like approxfun. — Lyngbakr
Thank you for your comment. I don't know this package but it seems that it is used for partial differential equations. There are no diffusion and advection terms in my equations. In my agent-based model, the habitat polygons don't interact. Only host individuals (i.e., susceptible, infected and recovered individuals) can interact with their environment. — Marine
In that case, ReacTran is overkill. Yes, you can have multiple non-interacting polygons by simply having arrays of length n (where n is the number of polygons) for each of S, I, and R and N that are combined into one array. (Personally, I'd calculate N at the end to save time.) Pass the function your state variable array, which would be 4 * n long. Then, inside your model split it into four parts. So, for example, S <- states[1:n], I <- states[(n+1):(2*n)]. The rest stays pretty much the same, except you'd need to have an initial state array that is 4 * n in length. — Lyngbakr
Thank you for your comment. But how can I apply the model to each polygon at a given time step? When I run my model, I obtain results for S, I, and R over all time period (in my case, times = seq(0, 5, by = 1)). — Marine

Lyngbakr Lyngbakr · Accepted Answer · 2017-07-06T15:54:02

Something like this.

# Library *not* require
library(deSolve)

# Parameters: number of polygons, beta, & gamma
n.polys <- 10 
beta <- sample(c(100, 200, 400, 600), n.polys, replace = TRUE)
gamma <- sample(c(25, 26, 27, 28), n.polys, replace = TRUE)

# Model defintion
sir_model <- function (t, Y, pars) { 
  # Break up state variable into parts
  S <- Y[1:n.polys]
  I <- Y[(n.polys+1):(2 * n.polys)]
  R <- Y[(2 * n.polys+1):(3 * n.polys)]

  # Calculate rate of change
  dSdt <- -beta * S * I 
  dIdt <- beta * S * I - gamma * I 
  dRdt <- gamma * I 

  # Return rates of change after concatenating them 
  return(list(c(dSdt, dIdt, dRdt))) 
} 

# Initial conditions
Y.ini <- c(S = rep(99, n.polys), I = rep(1, n.polys), R = rep(0, n.polys)) 

# Solve the model
ode(y = Y.ini, times = 0:5, func = sir_model, parms = NULL)

Applying a set of ordinary differential equations to each grid cell

1 Answers