Julia DifferentialEquations: add delay into ODE

[Luigi_Marongiu]

Hello,
I have a set of ODE to model the growth of bacteria nad phages:

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[2]) - (ρ * ψ * 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    = 1.0e-9      # initial infected population
v0    = 80.0        # initial phage population
r_s0    = 50000.0     # initial resident susceptible population
r_i0    = 1.0e-9      # initial resident infected population
r_v0    = 80.0        # initial resident phage population
i_s0    = 500.0       # initial invader susceptible population
i_i0    = 1.0e-9      # initial invader infected population
i_v0    = 80.0        # initial invader phage population
# 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]
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()) )

which gives me this figure

I would like to introduce the invading bacteria after a certain time lag, let’s say 1000 hours, so I added this modificaion:

condition(u, t, integrator) = t==1000
affect!(integrator) = integrator.u[4] += i_s0 
cb = DiscreteCallback(condition,affect!)
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

but the solution is the same and it is evident that the invader was present from the beginning.
How can I properly set the delay?
Thank you

[klaff]

You are setting that population to 500 initially in u0 (so it is not starting at zero), and adding 500 again at t==1000, but by then the population appears to be so large that 500 would not make a visible difference (if I’m lining up the state variables to the plot lines correctly).

[Luigi_Marongiu]

I see… the initial inoculum must be zero, of course. Thanks.
I changed i_s0 to 1e-9 (I was told it is better to have very small numbers than zero altogether) and

condition(u, t, integrator) = t==1000           # time of inoculum
affect!(integrator) = integrator.u[4] += 500    # amount of inoculum

but again the inoculum does not happen at t=1000 but at the beginning:

What am I missing?
Thank you

[Luigi_Marongiu]

Actually, if I use zero instead of the suggested 1e-9, the system works:

I think this solves the problem!
Thank you