Define an ODE which depends on an external function with the @ode_def macro

Specific Domains Numerics
,
1

Hello! I want to solve a system of ODEs that includes a parameter that varies over time. To obtain the value of the parameter I use an external function that uses the solution of another system of ODEs that has already been solved. I obtained the desired result by defining the second system with a generic function, but I want to see if it is possible to achieve it by defining it with the macro @ode_def of ParameterizedFunctions.jl.

I have thought of two ways to do this:

  1. Just define the second system with @ode_def and include the function in the respective equation, but I think this is not possible.
  2. Use a callback. This does not throw me any error but the result is not correct and I don’t know why.

Thank you for your time. I hope I have been clear enough. I attach the code below.

Code that generates the desired result using a generic function to define the second system of ODEs 
using DifferentialEquations
using ParameterizedFunctions

# Solve the first system of ODEs 
odeFB = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

x0 = ones(length(odeFB.syms))
p  = (1400, 24.15, 1100, 5.3, 4, 2, 2.3, 5.3, 2, 2.3)
tstops = 80000
tspan = (0,81074)
solFB = solve(ODEProblem(odeFB,x0,tspan,p); reltol = 1e-12, tstops = tstops);

# External function to calculate a constant of the equation du[2] of the second system 
function f(t)
    if t > 0
        i = findfirst(times .>= t)
        return ((((sol[i] - sol[i-1]) / (times[i] - times[i-1])) * (t .- times[i-1])) .+ sol[i-1])
    else 
        return x0[1]
    end
end 

# Solve the second system of ODEs 
function odeNF(du, u, p, t)
    Vmax, Ge, Km, k2, k3, k4, k5, km2, km4, km5 = p 
    du[1] = ((Vmax * Ge * u[1])/(Km + u[1])) - (k4 * u[1]) - (k2 * u[3] * u[1]) + (km2 * u[3]) + (km4 * u[2])
    du[2] = (k4 * f(t))- (k5 * u[2]) - (km4 * u[2]) + (km5 * u[3])
    du[3] =  - (k3 * u[3]) + (k5 * u[2]) - (km5 * u[3])
end 

k = findfirst(solFB.t .== tstops)
x0 = solFB[k]
sol = solFB[1, k:end]
times = solFB.t[k:end] .- solFB.t[k]
tspan = times[end]
solNF = solve(ODEProblem(odeNF, x0, tspan, p); reltol = 1e-12, saveat = times);
Just define the second system with @ode_def

The code would be exactly the same as the previous one, only that the second system is defined in this way:

odeNF = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * f(t)) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 km2 km4 km5

However, doing this throws an error → LoadError: TypeError: non-boolean (Term{Bool, Nothing}) used in boolean context in expression starting at Untitled-1.ipynb:28

Callback 
using DifferentialEquations
using ParameterizedFunctions

# Solve the first system of ODEs 
odeFB = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * Gp) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 k6 km2 km4 km5

x0 = ones(length(odeFB.syms))
p  = [1400, 24.15, 1100, 5.3, 4, 2, 2.3, NaN, 5.3, 2, 2.3]
tstops = 80000
tspan = (0,81074)
solFB = solve(ODEProblem(odeFB,x0,tspan,p); reltol = 1e-12, tstops = tstops);

# Set callback 
condition(u, t, integrator) = t > 0

function affect!(integrator)
    i = findfirst(times .>= integrator.t)
    integrator.p[8] = integrator.p[6] * ((((sol[i] - sol[i-1]) / (times[i] - times[i-1])) * (integrator.t .- times[i-1])) .+ sol[i-1])
end 

cb = DiscreteCallback(condition, affect!, save_positions=(false,false))

# Solve the second system of ODEs 
odeNF = @ode_def begin 
    dGp = ((Vmax * Ge * Gp)/(Km + Gp)) - (k4 * Gp) - (k2 * A * Gp) + (km2 * A) + (km4 * Gi)
    dGi = (k4 * k6) - (k5 * Gi) - (km4 * Gi) + (km5 * A)
    dA = - (k3 * A) + (k5 * Gi) - (km5 * A)
end Vmax Ge Km k2 k3 k4 k5 k6 km2 km4 km5 

k = findfirst(solFB.t .== tstops)
x0 = solFB[k]
p[8] = x0[1]
sol = solFB[1, k:end]
times = solFB.t[k:end] .- solFB.t[k]
tspan = times[end]
solNF = solve(ODEProblem(odeNF, x0, tspan, p); reltol = 1e-12, saveat = times, callback = cb);

This does not throw an error but the result is not correct.

2

This is a duplicate of another thread which has the answer: