I have this model of interactions between bacteria and phages:
using DifferentialEquations function payneJansenODE!(du, u, p, t) μ , φ, η, β, ω = p H = h = 0 #= du = susceptible du = infected du = 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) y = (η * u) println("t=", t, "\tu=", u, " [", χ, "] , \tu=", u, "\tu=", u) du = (μ * u) - (φ * u * u) - (H * u) du = χ*((μ * u) + (φ * u * u) - (η * u) - (H * u)) du = χ*((η * β * u) - (φ * u * u) - (ω * u) - (h * u)) 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 += vϕ # amount of inoculum cb1 = DiscreteCallback(condition1, affect1!) # extintion condition2(u,t,integrator) = u-1 affect2!(integrator) = integrator.u = 0 cb2 = ContinuousCallback(condition2, affect2!) condition3(u,t,integrator) = u-1 affect3!(integrator) = integrator.u = 0 cb3 = ContinuousCallback(condition3, affect3!) condition4(u,t,integrator) = u-1 affect4!(integrator) = integrator.u = 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 is correctly forced to zero when u<1, but not u and u.
cb4 , for some reason, do not take place, I added the term \chi to force everything to zero when u=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, u, and u at the same step (