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)