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

Faelens Ruben ruben.faelens at altran.com
Tue Apr 23 18:38:06 CEST 2013


Hi Vincent,

As Thomas already noted, ode() expects a single state vector; you should collate your different states S, I, R together in a single vector (as you have done in your second example).

If I read the code correctly:
S -> I flow of beta*S*I
I -> R flow of gamma*R
+ migration due to 'movrate' between subpopulations

Q: "Do you think it is alright?"
A: I think it is.

Ruben FAELENS 
Consultant
Altran Energy and Life Sciences division

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of vincent laperriere
Sent: dinsdag 23 april 2013 16:46
To: Thomas Petzoldt; Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] metapopulation SIR model

I tried out to adapt the n-species LV model as you mention in your last message but did not succeed.
[...]

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

[modified version with adapted linebreaks]
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))
[end modified version]



More information about the R-sig-dynamic-models mailing list