rmatrixmodelingodedesolve

Incorporating time step dependent rate/parameter in a metapopulation model

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)
Solution

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

    1. Integrate the dependency into the model function, e.g. by adding an additional differential equation for the parameter,
    2. Use a forcing function and interpolate from an external "signal",
    3. 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.