Adding an if then statement to condition initial value in ODE system; deSolve

Question

I'm trying to add an if then statement to condition the initial value of one of my state variables, and am using deSolve. Essentially, I want to introduce the 3rd ODE (in this case a 3rd species into a population) after the start of the simulation.

Here is what the code looks like without the condition:

Antia_3sp_Model <- function(t,y,p1){
  # Parms
  ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5] 
  # State vars
  Pi <- y[1]; Pj <- y[2]; I <- y[3]
  # ODEs
  dPi = ri*Pi - k*Pi*I
  dPj = rj*Pj - k*Pj*I
  dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
  list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)

Here's what I have so far, after trying to add in an if then statement, saying that before time = 50, the initial value of the 3rd state variable will be 0, and that at or above time = 50, the initial value of the 3rd state variable will be 1.

Antia_3sp_Model <- function(t,y,p1){
  # Parms
  ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5] 
  # State vars
  Pi <- y[1]; Pj <- y[2]; I <- y[3]
  
  if (t[i] < t[50]){
    Pj0 = 0
  }
  else if (t[i] >= t[50]){
    Pj0 = 1
  }
  
  # ODEs
  dPi = ri*Pi - k*Pi*I
  dPj = rj*Pj - k*Pj*I 
  dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
  list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)

Any suggestions?

Please let me know if I should add any additional information, and thank you so much for reading! :)

where you have if(t[i] < t[50]) you don't have a a variable i. I think you want ifelse instead which works on vectors Pj0 <- ifelse(t<t[50], 0, 1). — pdb
@pdb Thanks so much! That makes a lot of sense! I've tried your suggestion, and I'm not getting any error messages, but the output is all NAs after the initial time point––I'm not sure why that would be happening considering that the rest of the model can run when Pj=0. Any ideas? — MadelineJC
@pdb Oops––found a typo––I should have written TT instead of t. After fixing that, I'm getting a strange error message: 'Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (3) must equal the length of the initial conditions vector (2002)'... The function should be returning the same number of derivatives, even when the initial conditions are changed for a portion of the time series... Any ideas? — MadelineJC
This makes no sense. If t surpasses the value 5 (there is no array t available there), the initial values for t=0.1 are past history, they went long ago inactive. Changing them, even if that were possible, does not influence the further integration. Still, one could implement an external input this way, but then you would need to also use that input, for instance to switch on an interaction. — Lutz Lehmann
To introduce a sudden change of states at a certain point of time, you can do this with forcings or events. See the ?events help page or the following short tutorial. — tpetzoldt

tpetzoldt tpetzoldt · Accepted Answer · 2021-03-20T12:13:22

For me, it is not perfectly clear what is meant with the statement that the "initial value of the 3rd state variable" should be 1 for t >= 50. An initial value defines the start of a state variable, that then evolves by the differential equations. In the following, I show the following approaches:

The state variable Pj is initialized to a given value at t = 50. This can be handled by an event.
The state variable Pj receives additional external input at t >= 50. This can be handled with an external signal, also called a forcing function.

The first example shows the event mechanism, implemented as a data frame eventdat. It may also be implemented in a more flexible form with an event function.

Here I increased the "initial" state value at t=50 to 100, to make the effect more pronounced. Rounding of the time vector TT is done to avoid a warning (please ask if you want to know why).

library("deSolve")

Antia_3sp_Model <- function(t, y, p1){
  # Parms
  ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5]
  # State vars
  Pi <- y[1]; Pj <- y[2]; I <- y[3]

  # ODEs
  dPi <- ri*Pi - k*Pi*I
  dPj <- rj*Pj - k*Pj*I
  dI <- p*I*(Pi/(Pi + o) + Pj/(Pj + o))
  list(c(dPi, dPj, dI))
}

parms <- c(ri = 0.3, rj = 0.2, k = 0.001, p = 1, o = 1000)
N0 <- c(Pi = 1, Pj = 1, I = 1)
TT <- round(seq(0.1, 200, 0.1), 1)

## An "initial value" is the value at the beginning. We call the value during
## simulation the "state". If it is meant that the state should be changed at
## a certain point of time, it can be done with an event

# tp: initial value at t=50 set to 100 to improve visibility of effect (was 1)
eventdat <- data.frame(var = "Pj", time = 50, value = 100, method = "rep")

results <- lsoda(N0, TT, Antia_3sp_Model, parms, events=list(data=eventdat), verbose = TRUE)
plot(results, mfcol=c(1, 3))

A forcing function can be used to implement a time dependent parameter or to add a constant value to a state continuously. Note also the compact style of the ODE model. Whether to use the with function or not is a matter of taste. Both have their pros and cons.

But, whether to use an event or a forcing function makes a big difference.


Antia_3sp_Model <- function(t, y, p, import){
  with(as.list(c(y, p)), {
    dPi <- ri*Pi - k*Pi*I
    dPj <- rj*Pj - k*Pj*I + import(t)
    dI <- p*I*(Pi/(Pi + o) + Pj/(Pj + o))
    list(c(dPi, dPj, dI))
  })
}

signal <- approxfun(x=c(0, 50, max(TT)), y=c(0, 1, 1), method="constant", rule=2)

results <- lsoda(N0, TT, Antia_3sp_Model, parms, import=signal, verbose = TRUE)
plot(results, mfcol=c(1, 3))

Adding an if then statement to condition initial value in ODE system; deSolve

1 Answers