Hello,
I have this model of interactions between bacteria and phages:
using DifferentialEquations
function payneJansenODE!(du, u, p, t)
μ , φ, η, β, ω = p
H = h = 0
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
μ = growth rate (replication coef), a
φ = adsorption rate (transmission coef), b
η = lysis rate, k
β = burst size, L
ω = outflow (decay rate phage), m
H = host response against bacterium
h = host response against phage
=#
fx(x) = (x == 0) ? 0.0 : 1.0
χ = fx(u[1])
y = (η * u[2])
println("t=", t, "\tu[1]=", u[1], " [", χ, "] , \tu[2]=", u[2], "\tu[3]=", u[3])
du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
du[2] = χ*((μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2]))
du[3] = χ*((η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3]))
end
println(); println()
# parameters
mu = 0.5 # growth rate, a
phi = 1e-7 # adsorption rate, b
eta = 5 # lysis rate, k
beta = 100 # burst size, L
omega = 5 # outflow, m
s0 = 1000 # initial susceptible population, x0
i0 = 0.0 # initial infected population, y0
v0 = 0.0 # initial phage population
tϕ = 2.5 # time of inoculum, tφ
vϕ = 1e9 # amount of inoculum, vφ
tmax = 20.0 # duration
# model
u0 = [s0, i0, v0]
parms = [mu, phi, eta, beta, omega]
tspan = (0.0, tmax)
# inoculum
condition1(u, t, integrator) = t==tϕ # time of inoculum
affect1!(integrator) = integrator.u[3] += vϕ # amount of inoculum
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u,t,integrator) = u[1]-1
affect2!(integrator) = integrator.u[1] = 0
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u,t,integrator) = u[1]-1
affect3!(integrator) = integrator.u[2] = 0
cb3 = ContinuousCallback(condition3, affect3!)
condition4(u,t,integrator) = u[1]-1
affect4!(integrator) = integrator.u[3] = 0
cb4 = ContinuousCallback(condition4, affect4!)
modification = CallbackSet(cb1, cb2, cb3, cb4)
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])
I would like to account for the fact that bacteria can go extinct (when the population grows to numbers less than 1), so I am forcing the variables to zero. The tabulation shows that u[1] is correctly forced to zero when u[1]<1, but not u[2] and u[3].
Since cb3 and cb4 , for some reason, do not take place, I added the term \chi to force everything to zero when u[1]=0. But even this does not work.
Hence, I would like to ask whether is possible to modify multiple states of the integrator at once. That is, modify u[1], u[2], and u[3] at the same step (condition2 and affect2).
Thank you
