Setting initial (state) values for ODE system in compiled model (deSolve, Rcpp)

Question

I am struggling with a probably minor problem while calling compiled ODEs to be solved via the R package 'deSolve' and I seeking advice from more expert users.

Background

I have a couple of ODE systems to be solved with 'deSolve'. I have defined the ODEs in separate C++ functions (one for each model) I am calling through R in conjunction with 'Rcpp'. The initial values of the system change if the function takes input from another model (so basically to have a cascade).

This works quite nicely, however, for one model I have to set the initial parameters for t < 2. I've tried to do this in the C++ function, but it does not seem to work.

Running code example

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {
  
  List dn(3);
  
  double tau2 = parameters["tau2"]; 
  double Ae2_4 = parameters["Ae2_4"]; 
  double d2 = parameters["d2"]; 
  double N2 = parameters["N2"];
  
  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];
  
  // change starting conditions for t < 2
  if(t < 2) {
    n2 = (n2 * m4) / N2;
    m4 = n2;
    ne = 0;
    
  }
  
  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;    
  
  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

results <- deSolve::lsoda(
  y = state,
  times = 1:10,
  func = set_ODE,
  parms = parameters
)

print(results)
*/

The output reads (here only the first two rows):

  time            ne           n2            m4
1     1  1.000000e+01 0.000000e+00  0.000000e+00
2     2  1.000000e+01 2.169236e-07 -1.084618e-11

Just in case: How to run this code example?

My example was tested using RStudio:

• Copy the code into a file with the ending *.cpp
• Click on the 'Source' button (or <shift> + <cmd> + <s>)

It should work also without RStudio present, but the packages 'Rcpp' and 'deSolve' must be installed and to compile the code it needs Rtools on Windows, GNU compilers on Linux and Xcode on macOS.

Problem

From my understanding, ne should be 0 for time = 1 (or t < 2). Unfortunately, the solver does not seem to consider what I have provided in the C++ function, except for the ODEs. If I change state in R to another value, however, it works. Somehow the if-condition I have defined in C++ is ignored, but I don't understand why and how I can calculate the initial values in C++ instead of R.

Answer 1

I was able to reproduce your code. It seems to me that this is indeed elegant, even if it does not leverage the full power of the solver. The reason is, that Rcpp creates an interface to the compiled model via an ordinary R function. So back-calls from the slovers (e.g. lsoda) to R are necessary in each time step. Such back-calls are not for the "plain" C/Fortran interface. Here communication between solver and model takes place at the machine code level.

With this informational, I can see that we don't need to expect initialization issues at the C/C++ level, but it looks like a typical case. As the model function is simply the derivative of the model (and only this). The integration is done by the solver "from outside". It calls the model always with the actual integration state, derived from the time step before (roughly speaking). Therefore, it is not possible to force the state variables to fixed values within the model function.

However, there are several options how to resolve this:

• chaining of lsoda calls
• use of events

The following shows a chained approach, but I am not yet sure about the initialization of the parameters in the first time segment, so may only be part of the solution.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {

  List dn(3);

  double tau2 = parameters["tau2"];
  double Ae2_4 = parameters["Ae2_4"];
  double d2 = parameters["d2"];
  double N2 = parameters["N2"];

  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];

  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;

  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

## the following is not yet clear to me !!!
## especially as it is essentially zero
y1 <- c(ne = 0,
       n2 = unname(state["n2"] * state["m4"]/parameters["N2"]),
       m4 = unname(state["n2"]))


results1 <- deSolve::lsoda(
  y = y,
  times = 1:2,
  func = set_ODE,
  parms = parameters
)

## last time step, except "time" column
y2 <- results1[nrow(results1), -1]

results2 <- deSolve::lsoda(
  y = y2,
  times = 2:10,
  func = set_ODE,
  parms = parameters
)

## omit 1st time step in results2
results <- rbind(results1, results2[-1, ])

print(results)
*/

The code has also another potential problem as the parameters span several magnitudes from 1e-8 to 1e17. This can lead to numerical issues, as the relative precision of most software, including R covers only 16 orders of magnitude. Can this be the reason, why the results are all zero? Here it may help to re-scale the model.