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

I tried but did not work. I re-arranged as:

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(x, y) = (x < y) ? 0.0 : 1.0
    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(t,ε) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    println("t= ", t, "> ", X(t,ε), ": ", du[3])
end 
# parameters
mu  = 0.5       # growth rate, a
phi = 1e-7      # adsorption rate, b
eta = 5         # lysis rate, k
beta = 100      # burst size, L
omega = 5       #  outflow, m
s0 = 1000       # initial susceptible population, x0
i0 = 0.0        # initial infected population, y0
v0 = 0.0        # initial phage population
t_in = 2.5      # time of inoculum, tφ
v_in = 100      # amount of inoculum, vφ
tmax = 20.0     # duration
# instantiate
tspan   = (0.0, tmax)
parms   = [mu, phi, eta, beta, omega, t_in]
u0      = [s0, i0, v0]
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += v_in   # amount of inoculum
cb = DiscreteCallback(condition,affect!)
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

and I got:

julia>
t= 0.02> 0.0: 0.0
t= 0.0> 0.0: 0.0
t= 0.0011528391772627539> 0.0: 0.0
t= 0.02510627541594441> 0.0: 0.0
t= 0.03765941312391662> 0.0: 0.0
t= 0.10502364904863687> 0.0: 0.0
t= 0.1405332743433084> 0.0: 0.0
t= 0.20382196654005486> 0.0: 0.0
t= 0.21327524779360946> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23500388065939562> 0.0: 0.0
t= 0.32717504217939464> 0.0: 0.0
t= 0.37547864554281657> 0.0: 0.0
t= 0.6346915537961186> 0.0: 0.0
t= 0.7713301294032545> 0.0: 0.0
t= 1.0148606280227241> 0.0: 0.0
t= 1.051236198718852> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.124538601794153> 0.0: 0.0
t= 1.2650322206698672> 0.0: 0.0
t= 1.3386598925940385> 0.0: 0.0
t= 1.7337700218279322> 0.0: 0.0
t= 1.942043929289803> 0.0: 0.0
t= 2.3132498519821483> 0.0: 0.0
t= 2.368695996757698> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.4779124745446595> 0.0: 0.0
t= 2.639783896538048> 1.0: 0.0
t= 2.7246149091335132> 1.0: 0.0
t= 3.1798458226432613> 1.0: 0.0
t= 3.419811123806132> 1.0: 0.0
t= 3.847500516672675> 1.0: 0.0
t= 3.9113834628925255> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.037152472484574> 1.0: 0.0
t= 4.2222816752185315> 1.0: 0.0
t= 4.319301257399964> 1.0: 0.0
t= 4.8399400152491> 1.0: 0.0
t= 5.114383671461072> 1.0: 0.0
t= 5.603523726846436> 1.0: 0.0
t= 5.676585412203474> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.820092888234482> 1.0: 0.0
t= 6.024921293647171> 1.0: 0.0
t= 6.132264522152217> 1.0: 0.0
t= 6.708303337160416> 1.0: 0.0
t= 7.011949956102049> 1.0: 0.0
t= 7.55313835905778> 1.0: 0.0
t= 7.633974382156477> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.792595310607519> 1.0: 0.0
t= 8.015959631677706> 1.0: 0.0
t= 8.133016869457805> 1.0: 0.0
t= 8.761184230157331> 1.0: 0.0
t= 9.092309304086916> 1.0: 0.0
t= 9.682472469951104> 1.0: 0.0
t= 9.77062373676918> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.94339866563282> 1.0: 0.0
t= 10.182812927080112> 1.0: 0.0
t= 10.30828136356051> 1.0: 0.0
t= 10.981585962999713> 1.0: 0.0
t= 11.336504173041954> 1.0: 0.0
t= 11.969073769707178> 1.0: 0.0
t= 12.063559184117928> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.248606138382897> 1.0: 0.0
t= 12.502256710845652> 1.0: 0.0
t= 12.63518588785822> 1.0: 0.0
t= 13.348527257153213> 1.0: 0.0
t= 13.724549999534794> 1.0: 0.0
t= 14.394734137049063> 1.0: 0.0
t= 14.494837945860569> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.690752829766621> 1.0: 0.0
t= 14.95667016398437> 1.0: 0.0
t= 15.096027911328429> 1.0: 0.0
t= 15.84386708814722> 1.0: 0.0
t= 16.23807461659283> 1.0: 0.0
t= 16.940669440458176> 1.0: 0.0
t= 17.045614356315234> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.250894453033002> 1.0: 0.0
t= 17.52739867156927> 1.0: 0.0
t= 17.672304625668488> 1.0: 0.0
t= 18.449917291594815> 1.0: 0.0
t= 18.859819282055923> 1.0: 0.0
t= 19.590386295643068> 1.0: 0.0
t= 19.69950935087289> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.899629072326462> 1.0: 0.0
t= 19.91010889559556> 1.0: 0.0
t= 19.915600995490593> 1.0: 0.0
t= 19.94507339686609> 1.0: 0.0
t= 19.960609146882575> 1.0: 0.0
t= 19.988298464713438> 1.0: 0.0
t= 19.992434352185413> 1.0: 0.0
t= 20.0> 1.0: 0.0
t= 20.0> 1.0: 0.0
julia> 

Looks like the test worked but the condition did not…

However you rearrange the code, the way you’re using condition and affect! simply won’t work. In my last post, I showed you a different way that does work.

Sorry, but I do not understand the code. For starter, why affect and condition don’t work, since they work in other instances. Then, if I understand the code, you are splitting the analysis in two: one before t=2.5 and one after. Why not a single one with a split in the middle (the role of having condition)? And what is the role of @assert? Thank you

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

How would I do that?

I wrote this:

using DifferentialEquations
using PyPlot

# function
function baseInfect!(du, u, p, t)
    μ, κ, φ, ω, η, β = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end

# parameters
mu    = 0.47        # maximum growth rate susceptible (B. longum)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
beta  = 50.0        # burst size
tmax  = 4000.0      # time span 0-tmax
t_in  = 1000.0      # injection time
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 0.0         # initial phage population
v_in  = 80.0        # amount of injection

#=
SCENARIO 1
    one host one predator
=#
# instantiate solver
u0 = [s0, i0, v0]
parms = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)
# delay infection
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += v_in    # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = ODEProblem(baseInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])


In this case, condition and affect do work. Is the problem of time scale? In this example, the injection is at t=1000, in the original problem it was 2.5. In that case, would [tomaklutfu]'s suggestion be the easiest?

Interestingly, I tried with t=900 or 1100, and I get this:

So I got lucky the first time choosing 1000 as t_in…

That way could work too, when you set tstops=[t_in].

1 Like

Looks like it is fine this way. Thanks