Thank you Thomas,

I tried out to adapt the n-species LV model as you mention in your last message but did not succeed.
Is it possible to specify a model as:

                model <- function(t, S, I, R, parms) {......}

and then run it as:

                out <- ode(S0, I0, R0, times, model, parms)

?



Nevertheless, I tried another way to develop a SIR metapopulation model, like this:
do you think it is all right with the commands below?


N <- 10
S <- rep(1, times=N)
I <- rep(0, times=N)
R <- rep(0, times=N)
S[1] <- 0.99
I[1] <- 0.001
R[1] <- 0.009
state <- c(S, I, R)
movrate <- matrix(nr = N, nc = N, 0.0001)
movrate[1,2] <- 0.002
for (i in 1:N){movrate[i,i]=0}
FluxS <- matrix(nr = N, nc = N, 0.00)
for (i in 1:N){FluxS[i,i]=0}


SIRMETAPOP <- function (time, state, parms, N, movrate) {
  S <- state[1:N]
  I <- state[(N+1):(2*N)]
  R <- state[(2*N+1):(3*N)]
with (as.list(parms), {
## without migration terms
dS <- -beta*S*I
dI <- beta*S*I-gamma*I
dR <- gamma*I
## Migration terms
FluxS <- (S*movrate) # described as a matrix[i,j], flux of susceptibles from subpop i to subpop j
FluxI <- (I*movrate) # described as a matrix[i,j], flux of infectives from subpop i to subpop j
FluxR <- (R*movrate) # described as a matrix[i,j], flux of recovered from subpop i to subpop j
## with migration terms
dS <- dS - (rowSums(FluxS) - colSums(FluxS))
dI <- dI - (rowSums(FluxI) - colSums(FluxI))
dR <- dR - (rowSums(FluxR) - colSums(FluxR))
list(c(dS,dI,dR))
})}
pars <- c(beta=1000,gamma=365/13)
yini <- c(S=S,I=I,R=R)
times <-seq(from=0,to=60/365,by=1/365/4)
print(system.time(
out <- ode(y = yini, times = times, func = SIRMETAPOP, parms = pars, N=N, movrate=movrate)
))
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-", mfrow=c(3,3))


Best regards, Vincent Laperrière








________________________________
 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 : Lundi 22 avril 2013 21h28
Objet : Re: [R-sig-dyn-mod] metapopulation SIR model


Hi,

I would suggest to use an interaction matrix approach
similar to the multi-species Lotka-Volterra example, cf.:

http://desolve.r-forge.r-project.org/slides/tutorial.pdf#91

or the example below.

thpe


##----------------------------------------------------------
library(deSolve)

model <- function(t, n, parms) {
   with(parms, {
     dn <- r * n  + n * (A %*% n)
     list(dn)
   })
}

parms <- list(
   r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
   A = matrix(c(
     0.0, 0.0, -0.2, 0.0, # prey 1
     0.0, 0.0, 0.0, -0.1, # prey 2
     0.2, 0.0, 0.0, 0.0,  # predator 1; eats prey 1
     0.0, 0.1, 0.0, 0.0), # predator 2; eats prey 2
     nrow = 4, ncol = 4, byrow = TRUE)
)

times = seq(0, 500, 0.1)
n0  = c(n1 = 1, n2 = 1, n3 = 2, n4 = 2)

out <- ode(n0, times, model, parms)
plot(out)
	[[alternative HTML version deleted]]

