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

Luigi_Marongiu · December 11, 2020, 9:05am

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

ChrisRackauckas · December 11, 2020, 9:26am

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

Luigi_Marongiu · December 11, 2020, 9:32am

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ϕ])?

ChrisRackauckas · December 11, 2020, 9:43am

Yes.

Luigi_Marongiu · December 11, 2020, 9:54am

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…

ChrisRackauckas · December 11, 2020, 10:06am

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

Luigi_Marongiu · December 11, 2020, 10:43am

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ϕ])

ChrisRackauckas · December 11, 2020, 10:49am

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

Luigi_Marongiu · December 11, 2020, 11:05am

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

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

ChrisRackauckas · December 11, 2020, 11:12am

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

Luigi_Marongiu · December 11, 2020, 11:57am

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…

ChrisRackauckas · December 11, 2020, 12:02pm

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

Luigi_Marongiu · December 11, 2020, 1:52pm

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

ChrisRackauckas · December 11, 2020, 2:49pm

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

Luigi_Marongiu · December 11, 2020, 4:16pm

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?

ChrisRackauckas · December 11, 2020, 5:06pm

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

Luigi_Marongiu · December 11, 2020, 5:26pm

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ϕ])

ChrisRackauckas · December 11, 2020, 5:28pm

Typo. You made a boolean.

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

Luigi_Marongiu · December 11, 2020, 5:48pm

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:

ChrisRackauckas · December 11, 2020, 5:54pm

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.