Julia DifferentialEquations: add two delays into ODE

Hello,
I have a set of ODEs to model the invasion of a bacterium in presence of phages. I have inserted the invading bacterium at t=1000 using a delay.

function multiInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ, ω, β, Β, ρ, τ, η, Η = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage
    =#
    TOT = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * 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])

    du[4] = (ν * u[4]) * (1 - TOT/κ) - (ψ * u[4] * u[6]) - (ω * u[4])
    du[5] = (ψ * u[4] * u[6]) - (Η * u[5]) - (ω * u[5])
    du[6] = (Β * Η * u[5]) - (ρ * ψ * u[4] * u[6]) - (ω * u[6])
end

# parameters
mu    = 0.47        # maximum growth rate susceptible (B. longum)
nu    = 0.72        # maximum growth rate resistant (E. coli)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
psi   = 10.0^-9     # adsorption rate alternative
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
Eta   = 1.0         # lyse rate alternative
beta  = 50.0        # burst size
Beta  = 50.0        # burst size alternative
rho   = 0.6         # reinfection rate
tau   = 3.62        # latency time
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 80.0        # initial phage population
r_s0    = 50000.0     # initial resident susceptible population
r_i0    = 0.0         # initial resident infected population
r_v0    = 1000.0      # initial resident phage population
i_s0    = 0.0         # initial invader susceptible population
i_i0    = 0.0         # initial invader infected population
i_v0    = 100.0       # initial invader phage populationh
# implement
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, r_v0, i_s0, i_i0, i_v0]
parms = [mu, nu, kappa, phi, psi, omega, beta, Beta, rho, tau, eta, Eta]
# delay infection
condition(u, t, integrator) = t==1000           # time of inoculum
affect!(integrator) = integrator.u[4] += 500    # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

which gives:


Would be possible to add a second delay so that the phage of the invader is activated let’s say at t=2000? What would be the syntax?
Thank you

You can create a callback set with multiple callbacks.The doc is @ https://diffeq.sciml.ai/latest/features/callback_functions/#CallbackSet

I arranged like this:

condition_1(u, t, integrator) = t==1000           # time of inoculum
affect_1!(integrator) = integrator.u[4] += 500    # amount of inoculum
cb1 = DiscreteCallback(condition_1, affect_1!)
condition_2(u, t, integrator) = t==1100        # time of activation
affect_2!(integrator) = integrator.u[6] += 1000    # amount of phage
cb2 = DiscreteCallback(condition_2, affect_2!)
delay = CallbackSet(cb1, cb2)
# instantiate model
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=delay, tstops=[1000])

However unless I use the same time span in both cases, the second term does not kick in. For instance, at t=1100 for the phage I get:


If the syntax is correct, then is a problem with the model itself and the case is closed.
Thanks

I think you have two options here:

  1. Add also the second activation time to the tstops array e.g. tstopts=[1000, 1100]
  2. Simply use ContinuousCallback's and forget about tstops

What would be the syntax? I tried with:

julia> condition_1(u, t, integrator) = t==1000           # time of inoculum
condition_1 (generic function with 1 method)

julia> affect_1!(integrator) = integrator.u[4] += 500    # amount of inoculum
affect_1! (generic function with 1 method)

julia> cb1 = DiscreteCallback(condition_1, affect_1!)
DiscreteCallback{typeof(condition_1),typeof(affect_1!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}(condition_1, affect_1!, DiffEqBase.INITIALIZE_DEFAULT, Bool[1, 1])

julia> condition_2(u, t, integrator) = t==1001        # time of activation
condition_2 (generic function with 1 method)

julia> affect_2!(integrator) = integrator.u[6] += 10000    #
affect_2! (generic function with 1 method)

julia> cb = ContinuousCallback(condition_1, affect_1, condition_2, affect_2)
ERROR: UndefVarError: affect_1 not defined
Stacktrace:
 [1] top-level scope at none:1

Tx

Then, you will create two continuous callbacks and combine with callback set. Also their condition is not boolean but should cross zeros like 1000-t. In the documentation you’ll see further details.
Also: Do not forget bang at the end of function names: affect_1!, affect_2!.

1 Like

Should be something like this.

condition_1(u, t, integrator) = t-1000 
condition_2(u, t, integrator) = t-1100
affect_1!(integrator) = integrator.u[4] += 500 
affect_2!(integrator) = integrator.u[6] += 10000
cb1 = ContinuousCallback(condition_1, affect_1!)
cb2 = ContinuousCallback(condition_2, affect_2!)
delay = CallbackSet(cb1, cb2)


prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=delay)

Thank you, now it workes as wanted:

Sorry, just an update. I tried to expand the model by having a resident phage against the invader and adding a second phage at time t=1500:

function multiDoubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ1, ψ2, ω, β, Β1, Β2, η, Η1, Η2 = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage 1
    du[7] = invader phage 2
    =#
    TOT = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * 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])

    du[4] = (ν * u[4]) * (1 - TOT/κ) - (ψ1 * u[4] * u[6]) - (ψ2 * u[4] * u[7]) - (ω * u[4])
    du[5] = (ψ1 * u[4] * u[6]) + (ψ2 * u[4] * u[7]) - (Η1 * u[5]) - (Η2 * u[5]) - (ω * u[5])
    du[6] = (Β1 * Η1 * u[5]) - (ψ1 * u[4] * u[6]) - (ω * u[6])
    du[7] = (Β2 * Η2 * u[5]) - (ψ2 * u[4] * u[7]) - (ω * u[7])
end
psi1    = 10.0^-9       # adsorption rate alternative
Eta1    = 1.0           # lyse rate alternative
Beta1   = 50.0          # burst size alternative
psi2    = 10.0^-9       # adsorption rate alternative
Eta2    = 1.0           # lyse rate alternative
Beta2   = 50.0          # burst size alternative
r_s0    = 50000.0       # initial resident susceptible population
r_i0    = 0.0           # initial resident infected population
r_v0    = 1000.0        # initial resident phage population
i_s0    = 0.0           # initial invader susceptible population
i_i0    = 0.0           # initial invader infected population
i_va0   = 1000.0         # initial invader phage A population
i_vb0   = 0.0           # initial invader phage B population
# implement
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, r_v0, i_s0, i_i0, i_va0, i_vb0]
parms = [mu, nu, kappa, phi, psi1, psi2, omega, beta, Beta1, Beta2,
    eta, Eta1, Eta2]
# delay infection
condition_1(u, t, integrator) = t-1000              # time invader inoculum
condition_2(u, t, integrator) = t-1500              # time phage activation
affect_1!(integrator) = integrator.u[4] += 500      # amount inoculum
affect_2!(integrator) = integrator.u[7] += 500      # amount phage
cb1 = ContinuousCallback(condition_1, affect_1!)
cb2 = ContinuousCallback(condition_2, affect_2!)
delay = CallbackSet(cb1, cb2)
# instantiate model
prob = ODEProblem(multiDoubleInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=delay, tstops=[1000])

However, the second phage is there from the beginning and the resident is not:


Is there something wrong with the call?

Testet it and works perfectly, but look at your model:
u[7] is already in range of 1e9 before t=1500 so adding up 500 would not change anything

I see, that is due to u[4] and u[5]. I must include an extra parameter that neutralizes them until t=1500 is reached! Thank you…
However, this brings me another problem: how can I modify the same equation twice? I re-written the function as follows:

function multiDoubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ1, ψ2, ω, β, Β1, Β2, η, Η1, Η2, χ = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage 1
    du[7] = invader phage 2
    =#
    TOT = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * 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])

    du[4] = (ν * u[4]) * (1 - TOT/κ) - (ψ1 * u[4] * u[6]) - (ψ2 * u[4] * u[7]) - (ω * u[4])
    du[5] = (ψ1 * u[4] * u[6]) + (ψ2 * u[4] * u[7]) - (Η1 * u[5]) - (Η2 * u[5]) - (ω * u[5])
    du[6] = (Β1 * Η1 * u[5]) - (ψ1 * u[4] * u[6]) - (ω * u[6])
    du[7] = χ * (Β2 * Η2 * u[5]) - χ * (ψ2 * u[4] * u[7]) - (ω * u[7])
end

the parameters χ is 0 until the activation is reached, afterwards it takes 1. I then updated the execution as follows:

psi1    = 10.0^-9       # adsorption rate alternative
Eta1    = 1.0           # lyse rate alternative
Beta1   = 50.0          # burst size alternative
psi2    = 10.0^-9       # adsorption rate alternative
Eta2    = 1.0           # lyse rate alternative
Beta2   = 50.0          # burst size alternative
chi     = 0.0           # controller delayed activation
r_s0    = 50000.0       # initial resident susceptible population
r_i0    = 0.0           # initial resident infected population
r_v0    = 1000.0        # initial resident phage population
i_s0    = 0.0           # initial invader susceptible population
i_i0    = 0.0           # initial invader infected population
i_va0   = 1000.0         # initial invader phage A population
i_vb0   = 0.0           # initial invader phage B population
# implement
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, r_v0, i_s0, i_i0, i_va0, i_vb0]
parms = [mu, nu, kappa, phi, psi1, psi2, omega, beta, Beta1, Beta2,
    eta, Eta1, Eta2, chi]
# delay infection
condition_1(u, t, integrator) = t-1000              # time invader inoculum
condition_2(u, t, integrator) = t-1500              # time phage activation
condition_3(u, t, integrator) = t-1500              # control on activation
affect_1!(integrator) = integrator.u[4] += 500      # amount inoculum
affect_2!(integrator) = integrator.u[7] += 1000     # amount phage
affect_3!(integrator) = integrator.u[7] += 1.0      # activation
cb1 = ContinuousCallback(condition_1, affect_1!)
cb2 = ContinuousCallback(condition_2, affect_2!)
cb3 = ContinuousCallback(condition_3, affect_3!)
delay = CallbackSet(cb1, cb2, cb3)
# instantiate model
prob = ODEProblem(multiDoubleInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=delay, tstops=[1000])

So, I modify u[7] at t=1500 with cb2, and I re-modify u[7] with cb3. How does the solver know that at t=1500 the values of i_vb0 changes to 1000? and, more importantly, how does it know to change chi to 1.0?
The model I get is now:


That is, the introduction of the new phage kills everything (but should not affect the blue bacterium and its phage…). Thank you

You can make chi a function

chi = (t) -> 1000 ? 1.0 : 0.0

and call it as χ(t) in your multiDoubleInfect!

1 Like

I tried with:

unction multiDoubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ1, ψ2, ω, β, Β1, Β2, η, Η1, Η2, χ = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage 1
    du[7] = invader phage 2
    =#
    TOT = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * 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])

    du[4] = (ν * u[4]) * (1 - TOT/κ) - (ψ1 * u[4] * u[6]) - (ψ2 * u[4] * u[7]) - (ω * u[4])
    du[5] = (ψ1 * u[4] * u[6]) + (ψ2 * u[4] * u[7]) - (Η1 * u[5]) - (Η2 * u[5]) - (ω * u[5])
    du[6] = (Β1 * Η1 * u[5]) - (ψ1 * u[4] * u[6]) - (ω * u[6])
    du[7] = χ(t) * (Β2 * Η2 * u[5]) - χ(t) * (ψ2 * u[4] * u[7]) - (ω * u[7])
end
chi = (t) -> 1000 ? 1.0 : 0.0

but I got:

julia> soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])
ERROR: TypeError: non-boolean (Int64) used in boolean context

In julia, conditional branches only allow boolean so it should be t>1000 ? 1.0 : 0.0.

1 Like

I set:

Chi(x, y) = (x < y) ? 0.0 : 1.0

and passed it as a parameter:

function multiDoubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ1, ψ2, ω, β, Β1, Β2, η, Η1, Η2, χ, ε = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage 1
    du[7] = invader phage 2
    =#
    sumBact = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - sumBact/κ) - (φ * 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])

    du[4] = (ν * u[4]) * (1 - sumBact/κ) - (ψ1 * u[4] * u[6]) - (ψ2 * u[4] * u[7]) - (ω * u[4])
    du[5] = (ψ1 * u[4] * u[6]) + (ψ2 * u[4] * u[7]) - (Η1 * u[5]) - (Η2 * u[5]) - (ω * u[5])
    du[6] = (Β1 * Η1 * u[5]) - (ψ1 * u[4] * u[6]) - (ω * u[6])
    du[7] = χ(t, ε) * (Β2 * Η2 * u[5]) - χ(t, ε) * (ψ2 * u[4] * u[7]) - (ω * u[7])
end
t_in    = 1500          # time of phage activation 
# implement
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, r_v0, i_s0, i_i0, i_va0, i_vb0]
parms = [mu, nu, kappa, phi, psi1, psi2, omega, beta, Beta1, Beta2,
    eta, Eta1, Eta2, Chi, t_in]

This time there is a kick in:


Thank you