a long time has elapsed since our last correspondence. I was wondering, would it be possible to formulate a model for "fitode" that accepts a time dependent variable such as exploitation rates? My aim would be to elucidate (as you did with the 4 species example above), the parameter sets (if any) that "fitode" returns such that co-existence into the future is granted for a simple two-species modified Lotka-Volterra model. Essentially I would like to reproduce your answer for a simple model (2 vs 4 species) but adding harvesting rates as per the following example:

# Assessment data from the Baltic Sea 
ObsSSB <- data.frame(Herring = c(1932041, 1864349, 1672273, 1944689, 1905001, 1814679, 1626604, 1438139, 1520902, 1421057, 1266502, 1177136, 1090921, 1011718, 1013695, 856321, 714642, 647277, 675487, 649195, 651524, 539582, 483859, 452706, 417827, 363049, 354724, 339899, 332746, 367432, 377935, 424888, 461102, 478633, 477803, 536040, 565200, 557749, 598115, 629339, 695343, 643874, 579087, 597533, 581015, 460378, 364981),
                     Sprat = c(940000, 726000, 625000, 1044000, 695000, 377000, 227000, 199000, 254000, 394000, 616000, 605000, 570000, 461000, 403000, 423000, 556000, 775000, 1045000, 1360000, 1375000, 1429000, 1811000, 1777000, 1354000, 1353000, 1319000, 1196000, 942000, 829000, 1040000, 1324000, 1055000,  898000, 921000, 827000, 948000, 752000, 694000, 706000, 620000, 680000, 1077000, 1095000, 982000, 855000, 817000))
rownames(ObsSSB) <- 1974:2020
Fdata <- data.frame(Herring = c(0.1603, 0.17, 0.1579, 0.1469,0.1263, 0.1533, 0.1593, 0.1839, 0.1648, 0.2255, 0.2369, 0.2524, 0.2265, 0.2626, 0.2547, 0.3426, 0.3324, 0.3416, 0.299, 0.3367, 0.4076, 0.3872, 0.4094, 0.4564, 0.4848, 0.4127, 0.49, 0.4336, 0.3995, 0.3179, 0.2713, 0.2464, 0.2712, 0.2729, 0.2744, 0.2485, 0.2943, 0.2299, 0.1658, 0.1492, 0.2099, 0.2886, 0.3496, 0.352, 0.4474, 0.4951, 0.46),
                    Sprat = c(0.3696, 0.396, 0.3777, 0.336, 0.3408, 0.2627, 0.3005, 0.1828, 0.3033, 0.1358, 0.1796, 0.1659, 0.2091, 0.2689, 0.2391, 0.2155, 0.1369, 0.1734, 0.2004, 0.1625, 0.2638, 0.3485, 0.2984, 0.416, 0.3981, 0.37, 0.3193, 0.2866, 0.4008, 0.3958,  0.4933, 0.4536, 0.3848, 0.3793, 0.4115, 0.5091, 0.4403, 0.3968, 0.3462, 0.4192, 0.4341, 0.4538, 0.3455, 0.3961, 0.395, 0.3918, 0.368))
rownames(Fdata) <- 1974:2020
### Fitting a modified 2 species Lotka-Volterra model to observed SSB data ###
library(deSolve)
spp <- ObsSSB/10^6
# First guess of parameters
nospp <- ncol(spp)
data_lagged <- rbind(NA, spp) # Lagged data
db.bdt <- list()
for (i in 1:ncol(spp)) {
  db.bdt[[i]] <- log(append(spp[,i], NA)/data_lagged[,i])
}
modlst <- NULL # Linear models
for(i in 1:nospp) {
  moddat <- data.frame(db.bdt = db.bdt[[i]], 
                       data_lagged)[-nrow(data_lagged),]
  modlst[[i]] <- eval(parse(text = paste("lm(db.bdt~", 
                                         paste(colnames(data_lagged), 
                                               collapse = "+"), ", data=moddat)", sep = "")))
}
r_start <- unname(sapply(modlst, function(dat){coef(dat)["(Intercept)"]}))
A_start <- unname(sapply(modlst, function(dat){coef(dat)[-1]}))
# Model input
params <- c(r1 = r_start[1], r2 = r_start[2],
            a11 = A_start[1,1], a12 = A_start[2,1],
            a21 = A_start[1,2], a22 = A_start[2,2])
ini <- c(b1 = spp[1,1], b2 = spp[1,2])
tmax <- nrow(spp)
t <- seq(1,tmax,1)
# Model formulation
inv <- 1e-5 # Numerical fudge to avoid biomasses becoming negative.
HS2LLV <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    H1 <- H1(time)
    H2 <- H2(time)
    db1.dt = b1*(r1+a11*b1+a12*(b2^2))-H1*b1+inv
    db2.dt = b2*(r2+a22*b2+a21*(b1^2))-H2*b2+inv
    list(c(db1.dt, db2.dt))
  })
}
# Fdata
matplot(as.numeric(rownames(Fdata)),Fdata, type = "l", xlab = "Time (yr)", 
        ylab = "Fishing Mortality (1/yr)", lty = 1)
# Fit two splines to Fdata to smooth out bumps for fit
library(splines)
Sp.fit <- array(dim = c(nrow(Fdata),ncol(Fdata)))
for (i in 1:ncol(Fdata)) {
  x <- round(seq(1974,2020,length.out = 4), digits = 0)
  y <- Fdata[as.character(x),i]
  ## Fit a couple of quadratic splines with different degrees of freedom
  fit <- lm(y ~ bs(x, degree = 2, df=length(x)-1))
  x0 <- seq(min(x), max(x), by = 1)
  Sp.fit[,i] <- predict(fit, data.frame(x = x0))
}
matlines(1974:2020,Sp.fit, lty = 2, lwd = 2)
H1 <- approxfun(Sp.fit[,1])
H2 <- approxfun(Sp.fit[,2])
# Optimiser
LVmse = function(parms) {
  out = as.matrix(deSolve::ode(ini, 1:tmax, HS2LLV, parms, method="rk4")[,-1])
  RSS = sum((spp-out)^2, na.rm = TRUE)
  return(RSS)
}
# Running model and optimising parameters
optimres1 <- optim(par = params, fn = LVmse)
fit1 <- deSolve::ode(y = ini, time = 1:47, func = HS2LLV, parms = optimres1$par)
# Plot
extrafont::loadfonts(device = "all")
matplot(1974:2020, ObsSSB/10^6, type = "n", xlab = "Time (years)", 
        ylab = "SSB (1000 x tons)", ylim = c(0,2.025))
grid()
matplot(1974:2020, ObsSSB/10^6, pch = 16, col = c("gray", "steelblue"), add = TRUE)
matplot(1974:2020, fit1[,2:3], type = "l", add = TRUE, col = c("gray", "steelblue"), lwd = 2)
legend("topright", legend = c("Herring","Sprat"), 
       lty = c(1,2), col = c("gray", "steelblue"), lwd = 2, bty = "n")