Dear all,

I am trying to apply the classir SIR epidemic model in a metapopulation context (i.e. in several subpopulations of a metapopulation), using the package DeSolve.

At the moment, migration is not included between subpopulations (it will be in a next step).
I am just trying to replicate the same numerical simulation results in each subpopulation as those for the basic model as reference, but it does not work.

As I am not familiar with ODE modelling, could somebody tell we what's wrong with the following commands ?

thank you for your help.
Vincent


library(deSolve)

############ first version SIR 
parameters <- c(beta=1000,gamma=365/13)
state <- c(S=0.99,I=0.001,R=0.009) # initial conditions
sirmodel <- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
dS <- -beta*S*I
dI <- beta*S*I-gamma*I
dR <- gamma*I
# return the rate of change
list(c(dS, dI, dR))
}) # end with(as.list ...
}
times <- seq(from=0,to=60/365,by=1/365/4) # returns a sequence
out <- ode(y = state, times = times, func = sirmodel, parms = parameters)
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")


############ first version SIR METAPOP
n <- 10
parameters <- c(beta=1000,gamma=365/13)
S <- rep(0.99, times=n)
I <- rep(0.001, times=n)
R <- rep(0.009, times=n)
state <- c(S=S, I=I, R=R)
sirmodel <- function(t, state, parameters) {
with (as.list(c(state, parameters)),{
# rate of change in subpopulation i
dS <- rep(0,n)
dI <- rep(0,n)
dR <- rep(0,n)
for (i in 1:n)
{
dS[i] <- -beta*S[i]*I[i]
dI[i] <- beta*S[i]*I[i]-gamma*I[i]
dR[i] <- gamma*I[i]
}
# return the rate of change
return(list(c(as.vector(dS), as.vector(dI), as.vector(dR))))
})
}
times <-seq(from=0,to=60/365,by=1/365/4)
print(system.time(
 out <- ode(y = state, times = times, func = sirmodel, parms = parameters)
 ))
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
	[[alternative HTML version deleted]]

