Unequal derivative number and length of initial condition vectors when working with ODEs

Question

I am have issues with the ODE function in R. I keep getting the error "Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (3006) must equal the length of the initial conditions vector (6)". The function seems to be creating derivatives of the size of the tspan plus the number of initial condition vectors. Any advice on how to procede?

library(deSolve)
tend <- 300
dt <- 1
tspan <- seq(0, tend, by = dt)

# Parameters
rHI <- 0.4 # PI producers
rHN <- 0.4 # X2 free riders
rLI <- rHI * 0.1 # Xr resistant
rLN <- rHI * 0.1
K <- 10e13 # carrying capacity of tumor

delta0 <- 2
rhoS <- 1
rhoI <- 1
lambdaA <- 0.5
lambdaS <- 0.5
beta <- 1
eta <- 1e3
lambdaAS <- 0.02
BCG_dose <- 4e6
t0 <- 0      
tnow <- times[1]


## Initial Conditions
TV <- 1.4 * 10e10
HI0 <- 0.5 * TV
HN0 <- 0.0 * TV
LI0 <- 0.25 * TV
LN0 <- 0.25 * TV
A0 <- 0
S0 <- 0


##Setting Up 
initcond <- c(HI0, HN0, LI0, LN0, A0, S0)
names(initcond) <- c("HI", "HN", "LI", "LN", "A", "S")
##Setting Up BCG dose
interval <- 7
BCG <- rep(0, length(tspan))
BCG[seq(interval, length(tspan), by = interval)] <- BCG_dose
for (i in 2:length(tspan)) {
  if (tspan[i] - tnow >= interval) {
    # update tnow and BCG0 for exponential decay
    BCG0 <- BCG[i-1]
    t0 <- tnow
    tnow <- tspan[i]
    BCG[i] <- BCG0 * exp(-0.3 * (tnow - t0))
  } else {
    BCG[i] <- BCG[i-1]
  }
}
bladder_ODE <- function(t, x, params, BCG) {
  HI <- x[1]
  HN <- x[2]
  LI <- x[3]
  LN <- x[4]
  A <- x[5]
  S <- x[6]
  

  
  BCG0 <- BCG[which.max(tspan <= t)]
  dBCG <- -BCG0 * exp(-0.3 * (t - t0))
  BCG[which.max(tspan <= t)] <- BCG0 + dBCG
  
  delta <- delta0 * max((A/(S+A) - 0.5), 0)
  
  dHI <- rHI * HI * (1 - (HI + HN + LI + LN) / K) - delta * HI
  
  dHN <- rHN * HN * (1 - (HI + HN + LI + LN) / K)
  
  dLI <- rLI * LI * (1 - (HI + HN + LI + LN) / K)
  
  dLN <- rLN * LN * (1 - (HI + HN + LI + LN) / K)
  
  dA <- rhoS * (HI + LI) + beta * BCG / (eta + BCG) - lambdaA * A
  
  dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
  
  return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
}

out <- ode(y = initcond, times = tspan, func = bladder_ODE)

I have tried changing the position of where tspan is defined, however, I am still running into the same issue. I have also changed the size of the iteration from 0.1 to 1, however, it still does not work.

Answer 1

The problem is, that in the model bladder_ODE a variable BCG is used in the dA-equation. As BCG is a vector with 301 values, the equation is vectorized and dA is extended to 301 elements.

The code contained also two other small issues:

• in line 5, tspan is defined, but some lines below, tnow <- times[1] appeared, that may be tnow <- tspan[1]
• the function bladder_ODE <- function(t, x, params, BCG)contains an argument BCG that needs to be also used in the ode function call, see below.

With these changes, it is now possible to reproduce the error. I have put two cat-lines in the code to print the length of the affected variables.

bladder_ODE <- function(t, x, params, BCG) {
  HI <- x[1]
  HN <- x[2]
  LI <- x[3]
  LN <- x[4]
  A <- x[5]
  S <- x[6]

  BCG0 <- BCG[which.max(tspan <= t)]
  dBCG <- -BCG0 * exp(-0.3 * (t - t0))
  BCG[which.max(tspan <= t)] <- BCG0 + dBCG
  
  cat(length(BCG), "\n") # <----------
  
  delta <- delta0 * max((A/(S+A) - 0.5), 0)
  dHI <- rHI * HI * (1 - (HI + HN + LI + LN) / K) - delta * HI
  dHN <- rHN * HN * (1 - (HI + HN + LI + LN) / K)
  dLI <- rLI * LI * (1 - (HI + HN + LI + LN) / K)
  dLN <- rLN * LN * (1 - (HI + HN + LI + LN) / K)
  dA <- rhoS * (HI + LI) + beta * BCG / (eta + BCG) - lambdaA * A
  dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
  
  cat(length(dA), "\n") # <-----------
  
  return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
}

## aded argument BCG
out <- ode(y = initcond, times = tspan, func = bladder_ODE, BCG=BCG)

That then returns:

301 
301 
Error in checkFunc(Func2, times, y, rho) : 
  The number of derivatives returned by func() (306) must equal the length of the initial conditions vector (6)

Besides this, the complete BCG[which.max(....)] looks quite unusual. You may search StackOverflow or the whole internet for: desolve forcings

Another comment: The order of magnitude of the initial states is quite large (1e10), which an lead to numerical problems. In principle, it would be possible to set individual tolerances for each of the state variables. However, it is often easier to re-scale the problem. I usually recommend to rescale it to a range between 0.0001 to 10000 or even better between 0.001 and 1000.