How to debug ODE

Hello,
I have a set of ODEs to model a bacterial infection:

function kinInfect!(du, u, p, t)
    μ , φ, η, β, ω, χ, ε = p
    H = h = 0
    #=
    du[1]   = susceptible
    du[2]   = infected
    du[3]   = phages
    μ       = growth rate (replication coef), a
    φ       = adsorption rate (transmission coef), b
    η       = lysis rate, k
    β       = burst size, L
    ω       = outflow (decay rate phage), m
    H       = host response against bacterium
    h       = host response against phage
    =#
    x = χ(t, ε)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = (x * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end
Chi(x, y) = (x < y) ? 0.0 : 1.0
# parms
tmax    = 20.0    # time span 0-tmax
mu      = 0.5     # growth rate (replication coef), a
phi     = 1e-7    # adsorption rate (transmission coef), b
eta     = 5       # lysis rate, k
beta    = 100     # burst size, L
omega   = 5       # outflow (decay rate phage), m
tmax    = 20.0    # time span 0-tmax
t_in    = 2.5     # time of inoculum
v_in    = 1e8     # amount of inoculum
s0      = 1000.0  # initial susceptible population
i0      = 0.0     # initial infected population
v0      = 0.0     # initial phage population
tspan   = (0.0, tmax)
parms   = [mu, phi, eta, beta, omega, Chi, t_in]
u0      = [s0, i0, v0]
# delay infection
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += 1e8   # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# implement model
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

the function Chi is there to neutralize the first term of du[3] until the inoculum time is reached. The first time I ran it I got:


But the second time I got:

I don’t remember having changed anything between the first and second time (apart from removing the vertical line).
How can I debug an ODE function? Are there some techniques?
Thank you

I guess that it’s the plot(...) command messing it up (probably an erroneous assignment t_in=... )

I don’t know the reason but one thing to care about is that you have to stop at 2.5sec(or any time with this kind of discontinuity) to correctly get the changes in Chi(x, y) = (x < y) ? 0.0 : 1.0. Otherwise, the solver allowed to pass that time by large margin if accepted by its error metrics.

1 Like

that may be truly the thing. How would I do that? Tx

The t==t_in looks dodgy. Unless you get lucky with rounding errors, that will never be true at a discrete time step, and the inoculation won’t happen. The plotted phage number never gets to 1e8, so it looks like the inoculation isn’t happening.

If I were doing this, I would solve for 2.5s, add 1e8 to the returned solution, then solve again for the remaining time, with that as the initial condition, and with the parameter chi set to 1 instead of 0.

Simply adding more to tstops=[2.5,1000 #=And more of the discontinuity time=#] or by continuous callbacks.

Here’s the simplest way to do what I think you want to do.

tmax    = 20.0    # time span 0-tmax
t_in    = 2.5     # time of inoculum
v_in    = 1e8     # amount of inoculum

μ      = 0.5     # growth rate (replication coef), a
φ     = 1e-7    # adsorption rate (transmission coef), b
η     = 5       # lysis rate, k
β    = 100     # burst size, L
ω   = 5       # outflow (decay rate phage), m
H = 0      # host response against bacterium
h = 0    # host response against phage

s0      = 1000.0  # initial susceptible population
i0      = 0.0     # initial infected population
v0      = 0.0     # initial phage population

function kinInfect!(du, u, _, t)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = ((t > t_in) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    @assert du[3] ≈ (η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end

prob = ODEProblem(kinInfect!, [s0, i0, v0], (0, t_in))
soln1 = solve(prob, AutoVern7(Rodas5()))
s1, i1, v1 = soln1[end]
v1 += v_in
prob = ODEProblem(kinInfect!, [s1, i1, v1], (t_in, tmax))
soln2 = solve(prob, AutoVern7(Rodas5()))

tt = [soln1.t; soln2.t]
uu = zeros(length(tt), 3)
for (j,u) in pairs([soln1.u; soln2.u])
    uu[j,:] = u
end
plot(tt, uu)

Note:

  • Parameters defined at the top level are still available in the body of kinInfect!

  • The boolean expression t>t_in can be promoted to a number, 0 if false, 1 if true. This makes a lot of simulation code simpler. Edit: But it can be left out of this code. If the initial state is in the subspace with no phages and no infections, u[2] == u[3] == 0, then it stays there.

  • You can use the solution of one problem as the initial value for another.

  • The plotting code should be plot(soln1); plot!(soln2), but plot! has a bug and replaces the existing plot instead of adding to it.

1 Like