Multiple forcings in a multi-patch ode model - R package desolve and compiled C code

Question

I am trying to create an SEIR model with multiple patches using the package deSolve in R. At each time step, there is some movement of individuals between patches that can infect individuals in other patches. I also have an external forcing parameter that is specific to each patch (representing different environmental conditions). I've been able to get this working in base R, but given the number of patches and compartments and the duration of the model, I'm trying to convert it to compiled code to speed it up.

I've gotten the different patches working, but am struggling with how to incorporate a different forcing parameter for each patch. When forcings are provided, there is an automatic check checkforcings (https://rdrr.io/cran/deSolve/src/R/forcings.R) that doesn't allow for a matrix with more than two columns, and I'm not quite sure what the best workaround is for this. Write my own ode and checkforcings functions to override this? Restructure the forcings data once it gets into C? My final model has 195 patches so I'd prefer to be to automate it somehow so I am not writing out thousands of equations or hundreds of functions.

Also fine if the answer is just, do this in a different language, but would appreciate insight into what language I should switch to. Julia maybe?

Below is code for a very simple example that just highlights this "different forcings in different patches problem".

R Code

# Packages #########################################################

library(deSolve)
library(ggplot2); theme_set(theme_bw())
library(tidyr)
library(dplyr)


# Initial Parameters and things ####################################

times <- 1:500
n_patch <- 2
patch_ind <- 100
state_names <- (c("S", "I"))
n_state <- length(state_names)

x <-rep(0, n_patch*n_state)
names(x) <- unlist(lapply(state_names, function(x) paste(x, 
                                                         stringr::str_pad(seq(n_patch), width = 3, side = "left", pad =0), 
                                                         sep = "_")))

#start with infected individuals in patch 1
x[startsWith(names(x), "S")] <- patch_ind
x['S_001'] <- x['S_001'] - 5
x['I_001'] <- x['I_001'] + 5
x['I_002'] <- x['I_002'] + 20

params <- c(gamma = 0.1, betam = 0.2)

#seasonality
forcing <- data.frame(times = times,
                      rain = rep(rep(c(0.95,1.05), each = 50), 5))
new_approx_fun <- function(rain.column, t){
  approx_col <- approxfun(rain.column, rule = 2)
      return(approx_col(t))
    }
    
    rainfall2 <- data.frame(P1 = forcing$rain, 
                            P2 = forcing$rain+0.01)

# model in R

    r.mod2 <- function(t,x,params){
      # turn state.vec into matrix
  # columns are different states, rows are different patches
  states <- matrix(x,
                   nrow = n_patch,
                   ncol = n_state, byrow = F)
  
  S <- states[,1]
  I <- states[,2]
  N <- rowSums(states[,1:2])
  
  with(as.list(params),{
    #seasonal forcing
    rain <- as.numeric(apply(as.matrix(rainfall2), MARGIN = 2, FUN = new_approx_fun, t = t))
    
    dS <- gamma*I - rain*betam*S*I/N 
    dI <- rain*betam*S*I/N  - gamma*I
    
    return(list(c(dS, dI), rain))
  })
}

out.R2 <- data.frame(ode(y = x, times =times, func = r.mod2,
                        parms = params))

#create seasonality for C
ftime <- seq(0, max(times), by = 0.1)
rain.ft <- approx(times, rainfall2$P1, xout = ftime, rule = 2)$y
forcings2 <- cbind(ftime, rain.ft, rain.ft +0.01)

# C model
system("R CMD SHLIB ex-patch-season-multi.c")
dyn.load(paste("ex-patch-season-multi", .Platform$dynlib.ext, sep = ""))
out.dll <- data.frame(ode(y = x, times = times, func = "derivsc",
                          dllname = "ex-patch-season-multi", initfunc = "parmsc",
                          parms = params, forcings = forcings2, 
                          initforc = "forcc", nout = 1, outnames = "rain"))

C code

#include <R.h>
#include <math.h>
#include <Rmath.h>
// this is for testing to try and get different forcing for each patch //

/*define parameters, pay attention to order */

static double parms[2];
static double forc[1];

#define gamma parms[0]
#define betam parms[1]
//define forcing
#define rain forc[0]



/* initialize parameters */
void parmsc(void (* odeparms)(int *, double *)){
  int N=2;
  odeparms(&N, parms);
}

/* forcing */
void forcc(void (* odeforcs)(int *, double *))
{
  int N=1;
  odeforcs(&N, forc);
}

/* model function */
void derivsc(int *neq, double *t, double *y, double *ydot, double *yout, int *ip){
  
  //use for-loops for patches
  //define all variables at start of block
  int npatch=2;
  double S[npatch]; double I[npatch]; double N[npatch];
  int i; 
  for(i=0; i<npatch; i++){
    S[i] = y[i];
  };
  for(i=0; i <npatch; i++){
    int ind = npatch+i;
    I[i] = y[ind];
  };
  for(i=0; i<npatch; i++){
    N[i] = S[i] + I[i];
  };
  
  //use for loops for equations 
  {
    // Susceptible
    for(i=0; i<npatch; i++){
      ydot[i] = gamma*I[i] - rain*betam*I[i]*S[i]/N[i] ;
    };
    //infected
    for(i=0; i<npatch; i++){
      int ind=npatch+i;
      ydot[ind] = rain*betam*I[i]*S[i]/N[i] - gamma*I[i];
    };
  };
  
  yout[0] = rain;
  
}

Answer 1

The standard way for multiple forcings in compiled code of the deSolve package is described in the lsoda help page:

forcings only used if ‘dllname’ is specified: a list with the forcing function data sets, each present as a two-columned matrix

Such a list can be created automatically in a script.

There are also other ways possible with some creative C or Fortran programming.

For more complex models, I would recommend to use the rodeo package. It allows to specify dynamic models in a tabular form (CSV, LibreOffice, Excel), including parameters and forcing functions. The code generator of the package creates then a fast Fortran code, that can be solved with deSolve. An overview can be found in a paper of Kneis et al (2017), https://doi.org/10.1016/j.envsoft.2017.06.036 and a more extended tutorial at https://dkneis.github.io/ .