Search code examples
rif-statementodedesolve

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

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! :)

Solution

  • 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:

    1. The state variable Pj is initialized to a given value at t = 50. This can be handled by an event.
    2. 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))