 # Introduce delayed element in differential equations with Julia's DifferentialEquations.jl

Hello,
I have a set of ordinary differential equations (ODE) that describe the growth of bacteria susceptible and resistant to phage infection:

that generate this plot:

To note that at t=1000 h the resistant strain is introduced.
How can I set the ODEs with Julia’s DifferentialEquations.jl package to accommodate for this delayed introduction? Can I use DDE instead?
I have written the following code:

``````function dynamo!(du, u, p, t)
μ, ν, κ, φ, ω, η, β, ρ = p
#=
du = susceptible
du = infected
du = phages
du = resistant
=#
du = ((μ * u) * (1 - ((u+u+u)/κ))) - (φ * u * u) - (ω * u)
du = (φ * u * u)
- (η * u) - (ω * u)
du = (β * η * u) - (ρ * φ * u * u) - (ω * u)
du = ((ν * u) * (1 - ((u+u+u)/κ))) - (ω * u)
end

# set parameters
mu    = 0.16        # maximum growth rate susceptible strain
nu    = 0.12        # maximum growth rate resistant strain
kappa = 2.2e7       # maximum population density
phi   = 1.0e-9      # adsorption rate
omega = 0.05        # outflow
eta   = 0.25        # lyse rate
beta  = 50.0        # burst size
rho   = 0.6         # reinfection rate
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
r0    = 10000.0     # initial resistant population
i0    = 1.0e-9      # initial infected population
v0    = 80.0        # initial phage population
u0 = [s0, i0, v0, r0]
p = [mu, nu, kappa, phi, omega, eta, beta, rho]
tspan = (0.0, tmax)
prob = ODEProblem(dynamo!, u0, tspan, p)
soln = solve(prob)
``````

but I got this graph: To note that here I simply set the starting value of the resistant strain at 1/5 of the susceptible’s.
Thank you

You can have a callback change a parameter at time t=1000.

and how would i set up that? would it work also for DDE? Thank you

Message me as a reply so I get an email reminder. I will post some stuff in 40 minutes if I remember.

Here’s an example of using callbacks:

Everything in DifferentialEquations.jl works across differential equations, so these callbacks work not just for ODEs, but also SDEs, DAEs, DDEs, etc.

Thanks, I switched directly to the DDE and wrote:

``````function dynDelay!(du, u, h, p, t)
μ, ν, κ, ϕ, ω, β, ρ, τ = p
#=
du = susceptible
du = infected
du = phages
du = resistant
=#
history = h(p, t - τ)
delay_lysis =  ϕ * history * history * exp(-1 * ω * τ)
du = (μ * u) * (1 - (u+u)/κ) - (ϕ * u * u) - (ω * u)
du = (ϕ * u * u) - delay_lysis - (ω * u)
du = (β * delay_lysis) - (ρ * ϕ * u * u) - (ω * u)
du = ((ν * u) * (1 - ((u+u+u)/κ))) - (ω * u)
end
mu    = 0.16        # maximum growth rate susceptible strain
nu    = 0.12        # maximum growth rate resistant strain
kappa = 2.2e7       # maximum population density
phi   = 1.0e-9      # adsorption rate
omega = 0.05        # outflow
beta  = 50.0        # burst size
rho   = 1.0         # 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
r0    = 10000.0     # initial resistant population
h(t, p) = ones(4)
tspan = (0.0, tmax)
u0 = [s0, i0, v0, r0]
parms = [mu, nu, kappa, phi, omega, beta, rho, tau]
condition(u,t,integrator) = t==1000
affect!(integrator) = integrator.u += r0
cb = DiscreteCallback(condition,affect!)
prob = DDEProblem(dynDelay!, u0, h, tspan, p=parms; constant_lags = [tau])
algt = MethodOfSteps(Tsit5())
soln = solve(prob, algt, callback=cb)
``````

but the graph does not show the increase in total bacteria at t~1500: What am I missing?

You forgot `soln = solve(prob, algt, callback=cb,tstops=)`, so the callback didn’t trigger. See the discussion in the example.

I thought the trigger was hardcoded in the call `condition(u,t,integrator) = t==1000`; I updated. Also, I thought that the solver would have put on hold the data under conditioning (u) so I put starting value at 10000.0; actually makes more sense that the starting valu eis modified at t=1000, so I changed the starting value of u to ~0:

``````h(t, p) = ones(4)
tspan = (0.0, tmax)
u0 = [s0, i0, v0, 1.0e-9]
parms = [mu, nu, kappa, phi, omega, beta, rho, tau]
condition(u,t,integrator) = t==1000
affect!(integrator) = integrator.u += r0
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = DDEProblem(dynDelay!, u0, h, tspan, p=parms; constant_lags = [tau])
algt = MethodOfSteps(Tsit5())
soln = solve(prob, algt, callback=cb, tstops=)
``````

Now the plot is like this: there is a change in population densities although not quite in the original…

How come you aren’t starting with `u0 = [s0, i0, v0, 0.0]`? Also, why did you make `h(t, p) = ones(4)`? That contradicts what you’re saying about the initial condition of `u0 ~ 0`, because now you’re saying `u(t) ~ 1` for all `t < 0`. I assume that you want `h(t,p) = [s0, i0, v0, 1e-9]` or whatever, matching your initial condition?

1 Like

You are right, the ones(4) was a residual from a previous example. The use of 1.0e-9 instead of 0.0 was to avoid the risk of divisions by zero. Now the figure shows a better change at t=1000: Thanks