Thank you Thomas for your reply.

Actually, my aim is to model connected, but distinct subpopulations. Do you have an idea of how to use a spatially implicit SIR metapopulation model, using ODE and a specified matrix of migration rates between subpopulations?


Best regards,
Vincent.






________________________________
 De : Thomas Petzoldt <Thomas.Petzoldt@tu-dresden.de>

Cc : Special Interest Group for Dynamic Simulation Models in R <r-sig-dynamic-models@r-project.org> 
Envoyé le : Vendredi 19 avril 2013 17h30
Objet : Re: [R-sig-dyn-mod] metapopulation SIR model


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))
	[[alternative HTML version deleted]]

