# Julia DifferentialEquations: add two delays into ODE

**URL:** <https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362>\
**Category:** Numerics\
**Tags:** question, package, diffeq\
**Created:** [November 18, 2020, 8:06am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362 "2020-11-18T08:06:01Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 18, 2020, 8:06am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/1 "2020-11-18T08:06:01Z")

</div>

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.

```julia
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:

 ![Scenario_3](https://global.discourse-cdn.com/julialang/original/3X/a/2/a249bab7dea6429cc8df6bb6073f62b69ec5d387.png)  
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

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [November 18, 2020, 8:16am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/2 "2020-11-18T08:16:27Z")

</div>

You can create a callback set with multiple callbacks.The doc is @ [Event Handling and Callback Functions · DifferentialEquations.jl](https://diffeq.sciml.ai/latest/features/callback_functions/#CallbackSet)

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 18, 2020, 8:58pm UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/3 "2020-11-18T20:58:54Z")

</div>

I arranged like this:

```julia
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:

 ![Scenario_4](https://global.discourse-cdn.com/julialang/original/3X/0/0/0034591d127ea822502cd91556136652c1305d1d.png)  
If the syntax is correct, then is a problem with the model itself and the case is closed.  
Thanks

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [November 18, 2020, 9:59pm UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/4 "2020-11-18T21:59:52Z")

</div>

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`

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 19, 2020, 7:48am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/5 "2020-11-19T07:48:00Z")

</div>

What would be the syntax? I tried with:

```julia
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

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [November 19, 2020, 8:39am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/6 "2020-11-19T08:39:17Z")

</div>

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!`.

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [November 19, 2020, 8:43am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/7 "2020-11-19T08:43:40Z")

</div>

Should be something like this.

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 19, 2020, 9:41am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/8 "2020-11-19T09:41:00Z")

</div>

Thank you, now it workes as wanted:

 ![Scenario_4](https://global.discourse-cdn.com/julialang/original/3X/c/6/c615754b5bb2679b14bef3019cc06ef0621d976a.png)

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 20, 2020, 9:20am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/9 "2020-11-20T09:20:33Z")

</div>

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:

```julia
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:

 ![Scenario_6](https://global.discourse-cdn.com/julialang/original/3X/a/d/ad50f5235da2af770f0a5fd71b27df57844054e9.png)  
Is there something wrong with the call?

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [November 20, 2020, 9:43pm UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/10 "2020-11-20T21:43:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 21, 2020, 6:47am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/11 "2020-11-21T06:47:30Z")

</div>

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:

```julia
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:

```julia
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:

 ![Scenario_6](https://global.discourse-cdn.com/julialang/original/3X/7/0/70b0857509183b5b4b5f355648e121d352fa3723.png)  
That is, the introduction of the new phage kills everything (but should not affect the blue bacterium and its phage…). Thank you

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [November 21, 2020, 7:33pm UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/12 "2020-11-21T19:33:06Z")

</div>

> [@Luigi\_Marongiu](#):
>
> 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?

You can make `chi` a function

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

```

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

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 22, 2020, 7:01am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/13 "2020-11-22T07:01:08Z")

</div>

I tried with:

```julia
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
julia> soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])
ERROR: TypeError: non-boolean (Int64) used in boolean context

```

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [November 22, 2020, 7:32am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/14 "2020-11-22T07:32:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 22, 2020, 8:13am UTC](https://discourse.julialang.org/t/julia-differentialequations-add-two-delays-into-ode/50362/15 "2020-11-22T08:13:14Z")

</div>

I set:

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

```

and passed it as a parameter:

```julia
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:

 ![Scenario_6](https://global.discourse-cdn.com/julialang/original/3X/f/c/fc7eec0364095beccf7e14c65b4a680baf69cbd5.png)  
Thank you