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

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.