Variable always equal to 0 in my ode model

Question

I have written an ODE model in R using the ode() function from the deSolve package. The equation predicts the evolution of the mass of pesticide on foliar surface. app1 is equal to the pesticide's mass applied (on day 156) and the function should track the evolution of the toxicant. However, I noticed that the output for the variable mf (foliar mass) is always equal to 0, even though it should be changing over time (fram the day 156). I am not sure what is causing this issue. Can someone help me understand what I am doing wrong?

library(deSolve)

mod <- function(t, yini, param, Precip) {   
  app1 <- ifelse(t == 156, 1.12, 0)
  jgrow <- 32
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t > 168, 0.9, ifelse(t < 136, 0, NA)))
  cover <- 0.9 * igrow / jgrow
  P = Precip[t]
  
  with(as.list(c(yini, param)), { 
    
    mharv = ifelse(t== 201, mf, 0)
    
    
    d_MF <- (app1 * sa * (1-drift)*cover*sa)- (mf* exp(-kf*dt)) + (yf * (mf* exp(-kf*dt))) - mf * (1-exp(-fet * P * dt)) - mharv
    return(list(c(d_MF))) 
    
  })} 

params <- list(drift = 0.1, dt = 1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)
t <- seq(1, 365, by = 1)
Precip=rep(0, 365)
Precip[155:205] =c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)
y0 <- c(mf = 0)
sol <- ode(y = y0, times = t, func = mod, parms = params, Precip = Precip)

Answer 1

Here a part of a solution. The time was rounded with floor and precip converted to a forcing function with constant interpolation.

The igrow <- ifelse(...)-function looks strange anyway, especially as a value in the model should never be NA.

library(deSolve)

mod <- function(t, yini, param, f_precip) {   
  app1 <- ifelse(floor(t) == 156, 1.12, 0)   # use round, floor or trunc
  jgrow <- 32
  
  ## the following can also be converted to a forcing function
  ## but I suspect it may have some errors
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t > 168, 0.9, ifelse(t < 136, 0, NA)))
  
  cover <- 0.9 * igrow / jgrow
  P <- f_precip(t)
  
  with(as.list(c(yini, param)), { 
    mharv <- ifelse(floor(t) == 201, mf, 0) # use floor, trunc or round
    d_MF <- (app1 * sa * (1 - drift) * cover * sa) - (mf * exp(-kf * dt)) + 
            (yf * (mf * exp(-kf * dt))) - mf * (1 - exp(-fet * P * dt)) - mharv
    return(list(c(d_MF))) 
  })
} 

params <- list(drift = 0.1, dt = 1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)

t <- seq(1, 365, by = 1)
Precip <- rep(0, 365)
Precip[155:205] <- c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 
                     0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 
                     0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 
                     3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)

## define a forcing function by tabular interpolation
f_precip <- approxfun(t, Precip, method = "constant", rule = 2)

y0 <- c(mf = 0)
sol <- ode(y = y0, times = t, func = mod, parms = params, 
           f_precip = f_precip)

plot(sol)

Edit

To answer the comment of the OP, why "the function reaches its maximum on day 2 after application", we can apply the following modifications:

I the model function, replace app1 as follows. It is also a good idea to add some internal variables as extra outputs in return after the c() of the state vector.

mod <- function(t, yini, param, f_precip, f_app) {   
  #app1 <- ifelse(floor(t) == 156, 1.12, 0)
  app1 <- f_app(t)
    .....
    return(list(c(d_MF), P = P, app = app1, igrow=igrow))
  })
} 

Then define a new forcing:

f_app <- approxfun(c(0, 156,        156+1/24, 365),
                   c(0, 1.12 * 24,    0,        0), method="constant", rule=2)

We see that only the relevant time steps are needed.

The call to ode must then be modified to contain the new forcing. And it is very important to reduce the time step, so that the external forcing is not overlooked by the automatic step size algorithm of the solver. It can be done by setting smaller time steps in the t-vector or by an additional argument hmax.

sol <- ode(y = y0, times = t, func = mod, parms = params, 
           f_precip = f_precip, f_app = f_app, hmax=1/24)

plot(sol, mfrow=c(2, 2))

A small time step, however, slows down the simulation. To avoid this, the event mechanism can be used. Then you have to re-formulate the model that it does not contain app1 in the derivative, but instead calculate the effect of it on mf. It is then directly and instantaneously added to mf at the event time.