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