Incorporating time step dependent rate/parameter in a metapopulation model

Question

I would like to incoparate a rate that changes with times in a metapopulation model using code by B. Raynor. The code was originally published at https://rpubs.com/bhraynor/MetapopulationModel

In the code below, I want mu to take values in the vector mu = seq(from = 0, to = 1, length = 100)

Below is the code that runs well, but wont run if I uncomment the mu = seq(from = 0, to = 1, length = 100) vector and comment the the vector mu = c(1/1000, 1/1000, 1/1000).

I get error saying: Error in eval(substitute(expr), data, enclos = parent.frame()) : dims [product 3] do not match the length of object [100]

#load libraries
library(dplyr)
library(deSolve)

MODEL <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
  
    #Define initial conditions
    S = matrix(state[1:3], ncol=1)
    I = matrix(state[(4:6)], ncol=1)
    R = matrix(state[(7:9)], ncol=1)
    
    #Define params
    gamma <- matrix(parameters[paste0("gamma", 1:3)], ncol=1)
    alpha <- matrix(parameters[paste0("alpha", 1:3)], ncol=1)
    mu <- matrix(parameters[paste0("mu", 1:3)], ncol=1)

    #mu = seq(from = 0.001, to = 0.0015, length = 100)
    
    dS <- -mu*S + alpha*R
    dI <-  mu*S - gamma*I 
    dR <-  gamma*I - alpha*R
    
    output <- c(dS, dI, dR)
    list(output)
  })
}

init <- c(S1 = 100,  S2 = 100,  S3 = 100,  I1 = 10,  I2 = 10, I3 = 0, R1 = 0,  R2 = 0,  R3 = 0 )

#parameters
parms <- c(gamma1 = 0.2000000, gamma2 = 0.2500000, gamma3 = 0.3333333, mu1 = 0.0010000,
           mu2 = 0.0010000, mu3 = 0.0010000,alpha1 = 0.200000, alpha2 = 0.3000000, alpha3 = 0.4000000)

Time = 100
dt = 1 #Step size dt
times <- seq(0, Time, by = dt)

#Run simulation
out <- ode(y=init, times=times, func=MODEL, parms=parms)

tpetzoldt · Accepted Answer · 2024-05-11 08:06:21Z

2

A parameter is usually a fixed value. A time-dependent value can be implemented by several methods:

Integrate the dependency into the model function, e.g. by adding an additional differential equation for the parameter,
Use a forcing function and interpolate from an external "signal",
Use events, to modify the state variables. This method can be combined with method #1.

The methods #2 and #3 are described at the ?forcings and ?events help pages of the deSolve package and an online tutorial page.

If parameters are only a function of time, i.e. changes independently of the state variables method #3 (forcing function) can be used. Then it depends on how many parameters are time dependent and how. If only a small number of parameters is time variable, they can be interpolated with approxfun. Here approxfun is a special function that returns a function object, containing the data.

For performance reasons, only the nodes (= data points) should be given.

library(deSolve)

MODEL <- function(time, state, parameters, mu_t) {
  with(as.list(c(state, parameters)), {
  
    #Define initial conditions
    S = matrix(state[1:3], ncol=1)
    I = matrix(state[(4:6)], ncol=1)
    R = matrix(state[(7:9)], ncol=1)
    
    #Define params
    gamma <- matrix(parameters[paste0("gamma", 1:3)], ncol=1)
    alpha <- matrix(parameters[paste0("alpha", 1:3)], ncol=1)
    mu <- matrix(parameters[paste0("mu", 1:3)], ncol=1)

    #mu = seq(from = 0.001, to = 0.0015, length = 100)
    mu <- mu_t(time)
    
    dS <- -mu*S + alpha*R
    dI <-  mu*S - gamma*I 
    dR <-  gamma*I - alpha*R

    ## store mu as auxiliary variable for debugging and plotting
    list(c(dS, dI, dR), mu=mu)
  })
}

init <- c(S1 = 100,  S2 = 100,  S3 = 100,  I1 = 10,  I2 = 10, I3 = 0, R1 = 0,  R2 = 0,  R3 = 0 )

#parameters
parms <- c(gamma1 = 0.2000000, gamma2 = 0.2500000, gamma3 = 0.3333333, mu1 = 0.0010000,
           mu2 = 0.0010000, mu3 = 0.0010000,alpha1 = 0.200000, alpha2 = 0.3000000, alpha3 = 0.4000000)

## Method A: many time steps
## rule = 2 avoids NA if solver overshoots end of time
mu_t <- approxfun(x = 0:100, y = seq(0.01, 0015, length.out=101), rule=2)

## Method B: if the change is linear, only the nodes are necessary, 
## e.g. two data points in the given example
mu_t <- approxfun(x = c(0, 100), y = c(0.01, 0015), rule=2)

times <- seq(0, 100, by = 1)

#Run simulation
out <- ode(y=init, times=times, func=MODEL, parms=parms, mu_t=mu_t)

## plot one set of variables as example, plus the values of mu
plot(out, which=c("S1", "I2", "R1", "mu"))

This works also for more than one time dependent parameter, but if a large number of parameters is to be interpolated, more elegant approaches exist. Using matrix formulations, it is also possible to streamline the whole system of equations, but that's another question.

edited May 11, 2024 at 8:06

answered May 10, 2024 at 5:35

tpetzoldt

5,8382 gold badges14 silver badges33 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

ccc Over a year ago

Thanks, @tpetzoldt, I have been able to implement using the approx function in a single patch model. However, the same approach fails to work on a multi-patch metapopulation model, stating that Error in eval(substitute(expr), data, enclos = parent.frame()) : dims [product 3] do not match the length of object [100] I suspect that ode is failing to extract vectors with 3 elements each time. The code works if I manually specify the time, e.g. by saying mu[6]. Any idea how I can create a compartment that counts with time? so I can use it as mu[counter].

tpetzoldt Over a year ago

Please create a minimum reproducible example, e.g. a toy model with 1-3 states and 2-3 parameters and without special rfun function.

ccc Over a year ago

thanks once more for your feedback. I have edited the code in the question by removing special functions.

Collectives™ on Stack Overflow

Incorporating time step dependent rate/parameter in a metapopulation model

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related