Using dede from the R deSolve/FME packages to fit data to a compartment model

Question

I'm trying to fit the tetracycline data set from Bates & Watts to a compartment model which forms a system of first order differential equations. The system has an analytic solution but I want to use the dede function to estimate the parameters numerically.

I can get parameter estimates which are close to the ones published in Bates and Watts but I'm wondering if I have coded the problem correctly. Specifically, since Bates & Watts account for dead time in their solution, I'm concerned about whether I have coded the use of lagvalue() in the function called DiffEqns correctly.

My programming question relates to coding of the derivatives with lag time. They are currently coded as:

dy1 <- -theta1*y1lag dy2 <- theta1*y1lag - theta2*y2lag

However, I wonder if the derivatives should be coded instead as:

dy1 <- -theta1*y1lag*y[1] dy2 <- theta1*y1lag*y[1] - theta2*y2lag*y[2]

Thanks and regards,

# Analyze the tetracycline data set as a two-compartment model 
# (see Bates & Watts, "Nonlinear Regression Analysis and Its Applications")
## Note: the differential equations for the compartment model are:
##   dy1/dt = -theta1*y1 
##   dy2/dt =  theta1*y1 - theta2*y2
##   (see p. 169 in Bates & Watts)
# Load packages
library(FME)
# Create the tetracycline dataset (see p. 281 in Bates & Watts)
tetra <- structure(list(time = c(1, 2, 3, 4, 6, 8, 10, 12, 16), 
                        conc = c(0.7,1.2, 1.4, 1.4, 1.1, 0.8, 0.6, 0.5, 0.3)), 
                   row.names = c(NA, 9L), class = "data.frame")
# Observe that: A) "conc" = data for y2; B) there is no data for y1; C) data start at time =  1 instead of time = 0
# Create a differential equation model with dead time 
DiffEqns <- function(t, y, parms) {
  theta1 <- parms[1] # rate constant for y1
  theta2 <- parms[2] # rate constant for y2
  theta3 <- parms[3] # amount of y1 at time = 0
  theta4 <- parms[4] # parameter that accounts for dead time
  y1lag <- ifelse(t - theta4 < 0, 0, lagvalue(t - theta4, 1))
  y2lag <- ifelse(t - theta4 < 0, 0, lagvalue(t - theta4, 2))
  
  dy1 <- -theta1*y1lag
  dy2 <-  theta1*y1lag - theta2*y2lag
  
  return(list(c(dy1, dy2), y1lag = y1lag, y2lag = y2lag))
}
# Find a numerical solution for the system of delay differential equations using dede() from deSolve
time <- seq(from = 0, to = 16, by = 0.1)
Cost <- function(P) {
  theta1 <- P[1]
  theta2 <- P[2]
  theta3 <- P[3]
  theta4 <- P[4]
  theta <- c(theta1, theta2, theta3, theta4)
  yinit <- c(y1 = theta3, conc = 0)
  out <- dede(y = yinit, times = time, func = DiffEqns, parms = theta)
  modCost(model = out, obs = tetra)
}
theta <- c(theta1 = 0.1, theta2 = 0.2, theta3 = 5, theta4 = 0.2) # starting values for the parameters
yinit <- c(y1 = theta[3], conc = 0)
CompModFit2 <- modFit(f = Cost, p = theta, lower = c(0,0,0,0))
FMEtheta <-  coef(CompModFit2)
# Compare data to numerical model solution using parameters from modFit
dedeFitted <- dede(times = time,y = c(y1 = FMEtheta[3], conc = 0), func = DiffEqns, parms = FMEtheta)
plot(dedeFitted, obs=tetra)
# Parameters from FME are:
# theta1     theta2     theta3     theta4 
#0.1193617  0.6974401 10.7188251  0.2206997 
# Compare FME parameters to the parameter estimates published in Bates & Watts:
# theta1    theta2     theta3     theta4 
# 0.1488    0.7158    10.10    0.4123