[R-sig-dyn-mod] metapopulation SIR model

[Thomas Petzoldt]

On 19.04.2013 16:42, vincent laperriere wrote:
> dear list,
> sorry to come again with my metapopulation simulation problem, using
> DeSolve package.
>
> I seek to add migration terms to an ODE system describing the dynamics
> of a SIR epidemic into n subpopulations, how should I proceed?
> Migration terms involve migrations rates as defined by a n*n matrix.
>
>
> Frow your answer to my previous post below are the commands of the model
> without migration terms.
>
> Should I rewrite it completely in a matrix formulation to integrate
> migration terms?
>
> Should it look something like that or am I on the wrong way?
> This is inspired from LVmod2D (Karline Soetaert, Thomas Petzoldt, R.
> Woodrow Setzer,Package deSolve: Solving Dierential Equations in R:
> Package deSolve, Journal of Statistical Software, February 2010, Volume
> 33, Issue 9).
>

Yes, this example can be a good starting point, because L&V and SIR are 
quite similar, but you can also use the example below that we used as an 
example in a talk about deSolve last year. The example uses an 
additional package ReacTran for transport, but you can do this also 
without this package (just ask).

Please note: these examples below describe populations on a grid, but 
for modelling connected but distinct sub-populations, I would use a 
spatially implicit approach.

Hope it helps

Thomas P.


### A simple SIR 2D model ####################################
library(deSolve)
library(ReacTran)


## The model equations
SIR2D <- function (t, y, parms)  {
   S  <- matrix(y[1:N^2], N, N)
   I  <- matrix(y[N^2 + 1:N^2], N, N)
   R  <- matrix(y[2*(N^2) + 1:N^2], N, N)

   infect <- beta * I * S
   recovr <- nu * I

   dS <- -infect
   dI <- infect - recovr  + tran.2D(I, dx = dx, dy = dy, D.x = D, D.y = 
D)$dC
   dR <- recovr
   list(c(dS, dI, dR))
}

## grid
N     <- 81
zero <- rep(0, N)
ndx <- 1:N^2
dx <- 10/N
dy <- 10/N

## model parameters
D     <- 1e-3
S0    <- 100
I0    <- 1
R0    <- 0.0
beta  <- 0.001
nu    <- 0.02


# initial conditions
S  <- matrix(nrow = N, ncol = N, data = S0)
I  <- matrix(nrow = N, ncol = N, data = 0)
R  <- matrix(nrow = N, ncol = N, data = R0)

set.seed(123)
I[sample(1:length(I), 10)] <- I0

y <- c(S = S, I = I, R = R)

times <- seq(0, 200, 4)

system.time(
   out <- ode.2D (y = y, func = SIR2D, t = times, nspec=3, parms = NULL,
                 dim = c(N, N), method = "adams")
)


windows(width=15, height=5)

image(out, select=1:3, axes=FALSE, legend = TRUE, ask=FALSE,
   main=c("Susceptible", "Infectious", "Recovered"), mfrow=c(1,3))