deSolve: get value at (t-tau)

Question

I am trying to run a published model into deSolve, and I am wondering how I can get the value of B at t-tau for the third equation dH/dt <- phi*B(t-tau)-(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H I can't get it to run so for now I have removed the (t-tau) part and can run the model in R without it, but this is a crucial part for the mechanisms I want to study.

library(deSolve)
# Full model
# rigidode <- function(t, y, parms) {
#   df <- c*N*F-ya*(F+H)-yb*B
#   dB <- L*(f^2/(f^2+b^2))*(H/(H+v))-phi*B
#   dH <- phi*B(t-tau)-(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H
#   dF <- T(a)*(amin+amax*(b^2/(b^2+f^2))-sigma*(F/(F+H)))*H-mr*M(a)*F
#   list(c(dy1, dy2, dy3))
# }

#Now here is the numerical simulation
rigidode <- function(t, y, parms) {
  dy1 <- c*N*y[4]-ya*(y[4]+y[3])-yb*y[2]
  dy2 <- L*(y[1]^2/(y[1]^2+b^2))*(y[3]/(y[3]+v))-phi*y[2]
  dy3 <- phi*y[2]-(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3] #y[2] should be "y[2](t-tau) as B(t-tau)"
  dy4 <- Ta*(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]-m*y[4]
  list(c(dy1, dy2, dy3, dy4))
}

yini <- c(y1=1000, y2=0, y3=16000, y4=8000)
parms <- c(c<-0.033,N<-9,ya<-0.007,yb<-0.0018,L<-2000,b<-500,v<-5000,phi<-1/21,amin<-0.25,amax<-0.25,sigma<-1.3,m<-0.05,Ta<-0.75,tau<-12)
times <- seq(from = 0, to = 600, by = 1)
out <- ode (times = times, y = yini, func = rigidode, parms = parms)
head (out, n = 3)
plot(out)

Answer 1

Honestly, I'm not confident if my answer is correct, but at least give you useful information.
deSolve package have a function lagvalue() to access to past (lagged) values of state variables.
?lagvalue says it is to be used with function dede (but perhaps with ode()).

rigidode <- function(t, y, parms) {

  if (t < tau)
    lagged_y2 <- 0
  else
    lagged_y2 <- lagvalue(t - tau)[2]

  if(t == 100) print(lagged_y2)  # for check

  dy1 <- c*N*y[4]-ya*(y[4]+y[3])-yb*y[2]
  dy2 <- L*(y[1]^2/(y[1]^2+b^2))*(y[3]/(y[3]+v))-phi*y[2]
  dy3 <- phi*lagged_y2-(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]  # modified line
  dy4 <- Ta*(amin+amax*(b^2/(b^2+y[1]^2))-sigma*(y[4]/(y[4]+y[3])))*y[3]-m*y[4]
  list(c(dy1, dy2, dy3, dy4))
}


# parameters are the same as yours.
yini <- c(y1=1000, y2=0, y3=16000, y4=8000)
parms <- c(c<-0.033,N<-9,ya<-0.007,yb<-0.0018,L<-2000,b<-500,v<-5000,phi<-1/21,
           amin<-0.25,amax<-0.25,sigma<-1.3,m<-0.05,Ta<-0.75,tau<-12)
times <- seq(from = 0, to = 600, by = 1)

out <- dede(times = times, y = yini, func = rigidode, parms = parms) # use dede

# printed_value
[1] 37022.74
[1] 37022.74

out[100 - 12 + 1,]
#  time        y1        y2        y3        y4 
# 88.00 157425.09  37022.74  59146.39  12862.78