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[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