How to force variable to zero upon meeting condition in Julia DifferentialEquations

Hello,
I have this model of bacteria and phage interactions:

using DifferentialEquations
function payneJansenODE!(du, u, p, t)
    μ , φ, η, β, ω = p
    H = h = 0
    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
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
u0 = [s0, i0, v0]
parms   = [mu, phi, eta, beta, omega]
tspan = (0.0, tmax)
condition(u, t, integrator) = t==tϕ            # time of inoculum
affect!(integrator) = integrator.u[3] += vϕ    # amount of inoculum
cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[tϕ])

The model works but from the plot (the horizontal line is y=1):

it is apparent that the bacterium went extinct since there is no such thing as 1/10 of a bacterium.
It is possible to force u[1] to zero (which will bring also u[2] = 0) when a condition is reached? (in this case, when u[1] <1).
Thank you

1 Like

Yes, that would be a ContinuousCallback where condition(u,t,integrator) = u[1]-1 and affect!(integrator) = integrator.u[1] = 0.

2 Likes

Thank you, that is I need to create cb1 = ContinuousCallback(condition(u,t,integrator) = u[1]-1), modification = CallbackSet(cb, cb1) and then soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])?

Yes.

gave an error but I ran:

julia> # inoculum
       condition(u, t, integrator) = t==tϕ            # time of inoculum
condition (generic function with 1 method)
julia> affect!(integrator) = integrator.u[3] += vϕ    # amount of inoculum
affect! (generic function with 1 method)
julia> cb1 = DiscreteCallback(condition, affect!)
DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}(condition, affect!, DiffEqBase.INITIALIZE_DEFAULT, Bool[1, 1])
julia> # extintion
       affect!(integrator) = integrator.u[1] = 0
affect! (generic function with 1 method)
julia> cb2 = ContinuousCallback(condition(u,t,integrator) = u[1]-1)
ERROR: syntax: invalid keyword argument name "condition(u, t, integrator)"
Stacktrace:
 [1] top-level scope at none:0
julia> cb2 = ContinuousCallback(condition, affect!)
ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}(condition, affect!, affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0)
julia> # combine
       modification = CallbackSet(cb1, cb2)
CallbackSet{Tuple{ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}},Tuple{DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}}}((ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}(condition, affect!, affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0),), (DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}(condition, affect!, DiffEqBase.INITIALIZE_DEFAULT, Bool[1, 1]),))
julia> # run
       prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: [1000.0, 0.0, 0.0]
julia> soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])
retcode: Success


However, there are no u[2] and u[3]:

julia> getindex.(soln.u, 2)
12-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> getindex.(soln.u, 3)
12-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

The condition u[1]=0 must have kicked in too soon…

This has nothing to do with DiffEq… don’t define all of your functions with the same name :man_facepalming:. condition1(...) = ..., condition2(...) = ...

ops… I tried with this, but the condition u[1]=0 did not work:

# 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!)
# combine
modification = CallbackSet(cb1, cb2) 
# run
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])

Does u[1] ever equal 1? You didn’t clearly plot u[1].

I tried with condition2(u,t,integrator) = u[1]<1
and got:


where susceptible = u[1], infected=u[2] and `phage=u[3].

Yes, u[1]<1 isn’t u[1]-1.

But even if u[1]<1 is correct, the condition does not happen. The blue and red lines should go down and then disappear (since the axis are log). The model instead continues till the end and in fact the bacterium resurrects…

u[1]<1 is not correct and would not cause the condition to happen.

So what would be the condition that is triggered when u[1] gets decimal?

You have to pick a value to trigger it at. 0.9999999?

Sorry, I don’t understand. So it has to be u[1]<0.9999999? What would be the difference with u[1]<1 or [u]= 0.0001?

That’s why I said u[1] - 1

OK, so why it does not work?

condition1(u, t, integrator) = t==tϕ            # time of inoculum
affect1!(integrator) = integrator.u[3] += vϕ    # amount of inoculum
cb1 = DiscreteCallback(condition1, affect1!)
condition2(u,t,integrator) = u[1]-1    # extintion
affect2!(integrator) = integrator.u[1] == 0
cb2 = ContinuousCallback(condition2, affect2!)
modification = CallbackSet(cb1, cb2)
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])

Typo. You made a boolean.

affect2!(integrator) = integrator.u[1] = 0
1 Like

:woozy_face: gosh! Sorry about that. Now looks much better, in fact. However, u[2] and u[3] don’t go to zero so I also forced u[2] to zero but the phage does not disappear. But, if the syntax is OK, then is a problem with the model:

Yeah no worries. I was confused why you were looking for other things and did see that you had tried it with an error haha. My inner voice was “what the heck is this guy asking for?”

But yeah, getting “exact” extinctions in an ODE-based model is always hard. ODEs generally don’t go extinct: it’s a downside to that modeling choice. If you want that kind of behavior, continuous-time Markov chains are where people go. Catalyst.jl is made for that.

https://catalyst.sciml.ai/dev/

Though then the results are stochastic and require more time to interpret.

2 Likes