Skip to content

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

9 messages · vincent laperriere, Thomas Petzoldt, Faelens Ruben

Original
Thomas Petzoldt
# 
Dear Vincent,

R is a vectorized language. Loops are often not necessary, and named 
vectors work differently. Try the example below. In addition, it would 
be also possible to write the model completely in matrix formulation.

Thomas



n <- 3
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, I, R)

sirmodel <- function(t, state, parameters) {
   S <- state[1:n]
   I <- state[(n+1):(2*n)]
   R <- state[(2*n+1):(3*n)]

   with (as.list(parameters),{
     dS <- -beta*S*I
     dI <- beta*S*I-gamma*I
     dR <- gamma*I
   # return the rate of change
   list(c(dS, dI, 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 = "-", mfrow=c(3,3))
Thomas Petzoldt
#
On 19.04.2013 16:42, vincent laperriere wrote:
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))
2 days later
Thomas Petzoldt
# 
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)
Thomas Petzoldt
# 
Dear Vincent,

I appreciate that you are fighting with the topic, but not all modelling
questions can be solved by volunteer support. As far as I see, your
fluxrate approach seems to go into a promising direction, but don't
have time to check if "it is all right".
On 23.04.2013 16:46, vincent laperriere wrote:
This is not possible because the deSolve help file ?ode tells us:

"If func is an R-function, it must be defined as:
func <- function(t, y, parms,...)."
No, it would not follow the syntax of R/deSolve's "ode" function. States
have to be bound together as a vector. But, you have the option to
reorganize your states internally in the model function, as long as
states at the beginning and derivatives at the end remain a vector and
occur in the same order.
See above, sorry. Nevertheless, I've put some code together yesterday
evening, just for recreation ;-) Maybe this gives you some idea:

http://hhbio.wasser.tu-dresden.de/projects/modlim/src/SIR-metapop.R


I'm sure this can be made even simpler, more realistic or whatever.

Maybe that someone else on this list can give you more insight ?!!!
Comments are welcome.

thpe
Faelens Ruben
# 
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]