I am trying to create an epidemiological SEIR model and simulate the model compartments over time. More specifically, the code below code defines a function called Metapopulation_SEIR. This function represents a compartmental model known as the Metapopulation SEIR model, which is commonly used in epidemiology to simulate the spread of infectious diseases in multiple interconnected populations.See the code:
library(deSolve)
library(ggplot2)
# Metapopulation SEIR model function
Metapopulation_SEIR <- function(time, current_state, params, connectivity_matrix){
with(as.list(c(current_state, params)),{
# Calculate total population size in each subpopulation
N <- apply(current_state, 1, sum)
# Initialize rate of change vectors
dS <- rep(0, nrow(current_state))
dE <- rep(0, nrow(current_state))
dI <- rep(0, nrow(current_state))
dR <- rep(0, nrow(current_state))
dM <- rep(0, nrow(current_state))
# Loop through each subpopulation
for (i in 1:nrow(current_state)) {
# Calculate total number of individuals in subpopulation i
Ni <- N[i]
# Calculate rates of change for each compartment in subpopulation i
dSi <- -beta * S[i] * I[i] / Ni
dEi <- (beta * S[i] * I[i] / Ni) - sigma * E[i]
dIi <- sigma * E[i] - gamma * I[i] - mu * I[i]
dRi <- gamma * I[i]
dMi <- mu * I[i]
# Update rate of change vectors
dS[i] <- dSi
dE[i] <- dEi
dI[i] <- dIi
dR[i] <- dRi
dM[i] <- dMi
# Calculate movement of individuals between subpopulations
for (j in 1:nrow(current_state)) {
dS[i] <- dS[i] + connectivity_matrix[i, j] * (beta * S[j] * I[j] / Ni)
dE[i] <- dE[i] + connectivity_matrix[i, j] * (sigma * E[j])
dI[i] <- dI[i] + connectivity_matrix[i, j] * (gamma * I[j])
dR[i] <- dR[i] + connectivity_matrix[i, j] * (gamma * R[j])
dM[i] <- dM[i] + connectivity_matrix[i, j] * (mu * I[j])
}
}
# Return a list of rates of change for each compartment
return(list(c(dS, dE, dI, dR, dM)))
})
}
# Example inputs
# Number of subpopulations
num_subpopulations <- 3
# Initial compartment counts for each subpopulation
initial_state <- matrix(c(
S1 = 900, E1 = 10, I1 = 5, R1 = 85,
S2 = 950, E2 = 5, I2 = 2, R2 = 43,
S3 = 850, E3 = 15, I3 = 7, R3 = 100
), ncol = 4, byrow = TRUE)
# Parameters
params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)
# Connectivity matrix (specifies the movement of individuals between subpopulations)
# Example: Fully connected network where individuals can move between all subpopulations
connectivity_matrix <- matrix(0.1, nrow = num_subpopulations, ncol = num_subpopulations)
diag(connectivity_matrix) <- 1 # Individuals stay within their own subpopulation
# Time points
times <- seq(0, 365, by = 1)
# Solve the model
model <- ode(initial_state, times, Metapopulation_SEIR, params, connectivity_matrix)
Which give the following error in the last line:
Error in match.arg(method) : 'arg' must be NULL or a character vector
I tried to name the argument as follows:
model <- ode(initial_state, times, Metapopulation_SEIR, parms=params, connectivity_matrix)
But still getting the same error.
The code contained several mistakes:
ode call must be namedinitial_state must be passed as matrixwithA further issue is, that it is not necessary to use loops. R is a vectorized language, so that states can directly be formulated in matrix form, which makes the code more compact and much faster.
Below an intermediate version of the code, solving only the errors, but still contains loops.
However, I would also strongly recommend to re-formulate the model in matrix form. An example can be found here: https://github.com/tpetzoldt/covid/blob/master/models/sir_metapopulation.R This code contains essentially two versions of the state equations, one outcommented version with separate equations and one that uses matrix multiplication.
library(deSolve)
# Metapopulation SEIR model function
Metapopulation_SEIR <- function(time, current_state, params, connectivity_matrix, n_sub){
# tp: access to names in state **matrix** not possible with "with",
# this works only for vectors like 'parms'
#with(as.list(c(current_state, params)),{
with(as.list(params),{
# tp: reformat state vector as matrix
dim(current_state) <- c(n_sub, 5)
# tp: extract states
S <- current_state[,1]
E <- current_state[,2]
I <- current_state[,3]
R <- current_state[,4]
M <- current_state[,5]
# Calculate total population size in each subpopulation
N <- apply(current_state, 1, sum)
# Initialize rate of change vectors
dS <- rep(0, nrow(current_state))
dE <- rep(0, nrow(current_state))
dI <- rep(0, nrow(current_state))
dR <- rep(0, nrow(current_state))
dM <- rep(0, nrow(current_state))
# tp: loop is not necessary, R is a vectorized language
# Loop through each subpopulation
for (i in 1:nrow(current_state)) {
# Calculate total number of individuals in subpopulation i
Ni <- N[i]
# Calculate rates of change for each compartment in subpopulation i
dSi <- -beta * S[i] * I[i] / Ni
dEi <- (beta * S[i] * I[i] / Ni) - sigma * E[i]
dIi <- sigma * E[i] - gamma * I[i] - mu * I[i]
dRi <- gamma * I[i]
dMi <- mu * I[i]
# Update rate of change vectors
dS[i] <- dSi
dE[i] <- dEi
dI[i] <- dIi
dR[i] <- dRi
dM[i] <- dMi
# Calculate movement of individuals between subpopulations
for (j in 1:nrow(current_state)) {
dS[i] <- dS[i] + connectivity_matrix[i, j] * (beta * S[j] * I[j] / Ni)
dE[i] <- dE[i] + connectivity_matrix[i, j] * (sigma * E[j])
dI[i] <- dI[i] + connectivity_matrix[i, j] * (gamma * I[j])
dR[i] <- dR[i] + connectivity_matrix[i, j] * (gamma * R[j])
dM[i] <- dM[i] + connectivity_matrix[i, j] * (mu * I[j])
}
}
# Return a list of rates of change for each compartment
return(list(c(dS, dE, dI, dR, dM)))
})
}
# Example inputs
# Number of subpopulations
num_subpopulations <- 3
# Initial compartment counts for each subpopulation
# tp: add 5th column for M states
initial_state <- matrix(c(
S1 = 900, E1 = 10, I1 = 5, R1 = 85, M1=0,
S2 = 950, E2 = 5, I2 = 2, R2 = 43, M2=0,
S3 = 850, E3 = 15, I3 = 7, R3 = 100, M3=0
), ncol = 5, byrow = TRUE)
# tp: reformat initial_state into vector
initial_state <- as.vector(initial_state)
# Parameters
params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)
# Connectivity matrix (specifies the movement of individuals between subpopulations)
# Example: Fully connected network where individuals can move between all subpopulations
connectivity_matrix <- matrix(0.1, nrow = num_subpopulations, ncol = num_subpopulations)
diag(connectivity_matrix) <- 1 # Individuals stay within their own subpopulation
# Time points
times <- seq(0, 365, by = 1)
# Solve the model
# tp: add n_sub argument
model <- ode(initial_state, times, Metapopulation_SEIR, params, "lsoda",
connectivity_matrix = connectivity_matrix, n_sub=num_subpopulations)
# tp: alternative: works too, but extra arguments must always be named
model <- ode(initial_state, times, Metapopulation_SEIR, params,
connectivity_matrix = connectivity_matrix, n_sub=num_subpopulations)