# Introduce delayed element in differential equations with Julia's DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794
**Category:** General Usage
**Tags:** question, package
**Created:** [December 3, 2019, 12:36pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794 "2019-12-03T12:36:54Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 3, 2019, 12:36pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/1 "2019-12-03T12:36:54Z")
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:
```julia
function dynamo!(du, u, p, t)
μ, ν, κ, φ, ω, η, β, ρ = p
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
du[4] = resistant
=#
du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2]+u[4])/κ))) - (φ * 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 - ((u[1]+u[2]+u[4])/κ))) - (ω * u[4])
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
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2019, 12:41pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/2 "2019-12-03T12:41:13Z")
You can have a callback change a parameter at time t=1000.
---
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 3, 2019, 1:00pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/3 "2019-12-03T13:00:50Z")
and how would i set up that? would it work also for DDE? Thank you
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2019, 1:08pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/4 "2019-12-03T13:08:34Z")
Message me as a reply so I get an email reminder. I will post some stuff in 40 minutes if I remember.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2019, 2:09pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/5 "2019-12-03T14:09:31Z")
Here’s an example of using callbacks:
[https://docs.juliadiffeq.org/dev/features/callback\_functions/#Example-1:-Interventions-at-Preset-Times-1](https://docs.juliadiffeq.org/dev/features/callback_functions/#Example-1:-Interventions-at-Preset-Times-1)
Everything in DifferentialEquations.jl works across differential equations, so these callbacks work not just for ODEs, but also SDEs, DAEs, DDEs, etc.
---
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 3, 2019, 3:30pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/6 "2019-12-03T15:30:21Z")
Thanks, I switched directly to the DDE and wrote:
```julia
function dynDelay!(du, u, h, p, t)
μ, ν, κ, ϕ, ω, β, ρ, τ = p
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
du[4] = resistant
=#
history = h(p, t - τ)
delay_lysis = ϕ * history[1] * history[3] * exp(-1 * ω * τ)
du[1] = (μ * u[1]) * (1 - (u[1]+u[2])/κ) - (ϕ * u[1] * u[3]) - (ω * u[1])
du[2] = (ϕ * u[1] * u[3]) - delay_lysis - (ω * u[2])
du[3] = (β * delay_lysis) - (ρ * ϕ * u[1] * u[3]) - (ω * u[3])
du[4] = ((ν * u[4]) * (1 - ((u[1]+u[2]+u[4])/κ))) - (ω * u[4])
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[4] += 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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2019, 5:28pm UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/7 "2019-12-03T17:28:10Z")
> [@Luigi\_Marongiu](#):
>
> but the graph does not show the increase in total bacteria at t~1500:
You forgot `soln = solve(prob, algt, callback=cb,tstops=[1000])`, so the callback didn’t trigger. See the discussion in the example.
---
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 4, 2019, 8:12am UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/8 "2019-12-04T08:12:27Z")
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[4]) 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[4] to ~0:
```julia
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]
# add delay
condition(u,t,integrator) = t==1000
affect!(integrator) = integrator.u[4] += 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=[1000])
```
Now the plot is like this:
there is a change in population densities although not quite in the original…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 4, 2019, 8:44am UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/9 "2019-12-04T08:44:17Z")
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[4] ~ 0`, because now you’re saying `u(t)[4] ~ 1` for all `t < 0`. I assume that you want `h(t,p) = [s0, i0, v0, 1e-9]` or whatever, matching your initial condition?
---
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 4, 2019, 9:54am UTC](https://discourse.julialang.org/t/introduce-delayed-element-in-differential-equations-with-julias-differentialequations-jl/31794/10 "2019-12-04T09:54:52Z")
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