How to avoid instability in DifferentialEquations?

Hello,
I am trying to model the interaction between bacteria and phages but using certain combinations of parameters (namely the time and amount of phage inoculum) I get instability and the integrator aborts the session.
I am using a function 𝔪 (for Monod) with a delay injection (condition_2) and an extinction cause that sends a bacterium to zero (condition_1). Namely, if the time of injection is too early, I get:

Otherwise I get extinction, as expected:

How do I handle instability and allow the integrator to reach the end of the procedure?
How do I choose the good parameters for the analysis?
Thank you

function 𝔪!(du, u, p, t) # modified MONOD
    μ, ν, k, K, y, Y, σ, ω, ϕ, β = p
    #=
    p = max growth rate A & B, Monod constants A & B, yield A & B,
        [initial substrate], outflow, adsorption rate, burst size
    du[1] = substrate
    du[2] = species A
    du[3] = species B
    du[4] = phages
    =#
    growth_a = (u[1] * u[2]) / (k + u[1])
    growth_b = (u[1] * u[3]) / (K + u[1])
    du[1] = ((σ-u[1])*ω) - ((μ/y)*growth_a) - ((ν/Y)*growth_b)
    du[2] = (μ * growth_a) - (ϕ*u[2]*u[4]) - (ω*u[2])
    du[3] = ν * growth_b - (ω*u[3])
    du[4] = (β*ϕ*u[2]*u[4]) - (ϕ*u[2]*u[4]) - (ω*u[4])
end
mu    = 0.81        # maximum growth rate susceptible (/h)
nu    = 0.91        # maximum growth rate resistant   (/h)
kappa = 3.0*10^-6   # Monod constant (half saturation) (g/L) Su
Kappa = 3.1*10^-4   # Monod constant (half saturation) (g/L) Re
yld   = 2.5*10^10   # yield susceptible (cell/g)
Yld   = 3.8*10^10   # yield resistant (cell/g)
phi   = (3.7^-10)*360  # adsorption rate
beta  = 152.0       # burst size
omega = 6.0*10^-2   # outflow (/h)
s0    = 1*10^-4     # [initial substrate] (g/L)
su0   = 1.0*10^2    # initial susceptible population
re0   = su0*200     # initial resistant population
p0    = 0           # initial phage population
p_in  = 10^4          # phage inoculum
tmax  = 60.0        # time span 0-tmax
t_in  = 10          # time of phage inoculum
# instantiate solver
u0 = [s0, su0, re0, p0]
parms = [mu, nu, kappa, Kappa, yld, Yld, s0, omega, phi, beta]
tspan = (0.0, tmax)
# extintion
condition1(u, t, integrator) = u[2]-1
function affect1!(integrator)
  integrator.u[2] = 0
  integrator.u[4] = 0
end
cb1 = ContinuousCallback(condition1, affect1!)
# inoculum
condition_2(u, t, integrator) = t-t_in              # time 
affect_2!(integrator) = integrator.u[4] += p_in     # amount 
cb2 = ContinuousCallback(condition_2, affect_2!)
# run
delay = CallbackSet(cb1, cb2)
prob = ODEProblem(𝔪!, u0, tspan, parms)
soln = solve( prob, AutoVern7(Rodas5()), callback=delay, tstops=[t_in] )

Start with PSA: How to help yourself debug differential equation solving issues

4 Likes

Thank you! I have tried the new solvers TRBDF2() , KenCarp4() , or QNDF() but I got the same issue reported in the first figure. I then increased the tolerance with soln = solve( prob, AutoVern7(Rodas5()), callback=delay, abstol=1e-20, reltol=1e-20, tstops=[t_in] ) and I got it through:

So the case is closed.
However, I think the problem is not so much the solver or the ODEs (since they work given certain values) but how to pick the correct parameters. This is the same issue as in my previous post: how to get the range of values for the parameters?
Ideally, I would need a program that will output something of this sort:

↓quantity	→time
		    1	2	3	…
10		    No	No	Yes	…
20		    No	Yes	Yes	…
30		    No	Yes	No	…
…

Is it possible to get such a list?
Thanks

You can sample the space to build one.