deSolve ODE Integration Error, am I using the wrong function?

Question

I'm attempting to solve a set of equations related to biological processes. One equation (of about 5) is for a pharmacokinetic (PK) curve of the form C = Co(exp(k1*t)-exp(k2*t). The need is to simultaneously solve the derivative of this equation along with some enzyme binding equations and initial results where not as expected. After troubleshooting, realized that the PK derivative doesn't numerically integrate by itself, if k is negative using the desolve ode function. I've attempted every method (lsode, lsoda, etc) in the ode function, with no success. I've tried adjusting rtol, it doesn't resolve.

Is there an alternative to the deSolve ode function I should investigate? Or another way to get at this problem?

Below is the code with a simplified equation to demonstrate the problem. When k is negative, the integrated solution does not match the analytical result. When k is positive, results are as expected.

First Image, result with k=0.2: Analytical and Integrated results match when k is positive

Second Image, result with k=-0.2: Integrated result does not match analytical when k is negative

library(deSolve)

abi <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        
        dI <- k*exp(k*t)
                list(c(dI))
    })
}

k <- c(-0.2)

times <- seq(0, 24, by = 1)

I_analytical <- exp(k*times)

parameters <- c(k)
state <- c(I = 0)

out <- ode(y = state, times = times, func = abi, parms = parameters)

plot(out)
points(I_analytical ~ times)

It was pointed out that the initial condition easily resolves the above example, which is very helpful. Here is the equation I can't accurately integrate, I've tried a few different initial conditions without real success.

library(deSolve)

## Chaos in the atmosphere
CYP <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        #dE <- ksyn - (kdeg * E) + (k2 * EI) - (k1 * E * I)
        #dEI <- (k1 * E * I) - (k2 * EI) + (k4 * EIstar) - (k3 * EI)
        #dEIstar <- (k3 * EI) - (k4 * EIstar)
        #dOcc <- dEI + dEIstar
        dI <- a*tau1*exp(tau1*t) + b*tau2*exp(tau2*t) + c*tau3*exp(tau3*t)
            #list(c(dE, dEI, dEIstar, dOcc, dI))
          list(c(dI))
    })
}

ifit <- c(-0.956144311,0.82619445,0.024520276,-0.913499862,-0.407478829,-0.037174745)
a = ifit[1]
b = ifit[2]
c = ifit[3]
tau1 = ifit[4]
tau2 = ifit[5]
tau3 = ifit[6]


parameters <- c(ksyn = 0.82, kdeg = 0.02, k1 = 2808, k2 = 370.66, k3 = 2.12, k4 = 0.017, a, b, c, tau1, tau2, tau3)

#state <- c(E = 41, EI = 0, EIstar = 0, Occupancy = 0, I = 0.0)
state <- c(I=-0.01)
times <- seq(0, 24, by = .1)
out <- ode(y = state, times = times, func = CYP, parms = parameters)

I_analytical <- a*exp(tau1*times) + b*exp(tau2*times) + c*exp(tau3*times)

plot(out)
points(I_analytical ~ times)

Target curve and the ode solution line.

Answer 1

The initial value should be

state <- c(I= a + b + c)
#state <- c(I = 1)

Answer 2

The first script contains several issues. The most important two are that (1) the model function (abi) must contain the derivative, not an integrated function, while (2) the analytically integrated model missed I_0 that results from the integration constant.

Let's assume a first order decay model

dI/dt = k I

then analytical integration yields

I_t = I_0 exp(kt)

The code is then:

library(deSolve)

abi <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    # dI <- k*exp(k*t) # original
    dI <-  k * I       # corrected, should be the dervivative
    list(c(dI))
  })
}

k <- -0.2 # simplified, c() was not necessary
times <- seq(0, 24, by = 1)

# correction: set I0 to a value > zero
I0 <- 10

# I_analytical <- exp(k*times)    # original
I_analytical <- I0 * exp(k*times) # corrected, multiplied with I0

#state <- c(I = 0) # original
state <- c(I = I0) # corrected
parameters <- c(k = k)

out <- ode(y = state, times = times, func = abi, parms = parameters)

plot(out)
points(I_analytical ~ times)

This code can be further simplified if you want.