I am modeling the interaction between bacteria and phages with DifferentialEquations and making an interactive plot with GLMakie.
The plot works, but is it possible to catch instability in the ODEs so that I can run a logic test to display a label on the plot? In essence, if the ODEs become unstable, I would like to plot a label “INSTABILITY DETECTED” in the plot itself.
In the example below, instability is detected when the kappa value is increased. How can I catch this event in a variable?
ALSO: in theory, kappa should be displayed with two digits; instead, it goes into a floating point. How can I force two digits?
Thank you and happy Easter
## packages and functions
using GLMakie, DifferentialEquations
function odeSolver!(du, u, p, t)
μ, κ, ω, δ, η, β, λ = p
# μ bacterial growth rate
# κ bacterial carrying capacity
# ω system wash-out rate
# δ phagial infection rate
# η phagial lysis rate (inverse latency)
# β phagial burst size
# λ phagial decay rate
# du[1] = S; du[2] = I; du[3] = V
ρ = 1 - (u[1] + u[2])/κ # rho: logistic factor
ϡ = (δ*u[1]*u[3]) # upsampi : infected bacteria
du[1] = (μ*u[1]*ρ) - ϡ - (ω*u[1])
du[2] = ϡ - (η*u[2]) - (ω*u[2])
du[3] = (β*η*u[2]) - ϡ - (λ*u[3]) - (ω*u[3])
end
## initial parameters
begin
mu = 0.16 # μ: default growth rate
kappa = 2.2e7 # κ: default carrying capacity
omega = 0.05 # ω: default outflow
eta = 0.025 # η: default τ reciprocal
delta = 1e-9 # δ: default adsorption rate
beta = 50 # β: default burst size
lambda = 0 # λ: default decay rate
s0 = 1e5 # default starting amount of naive bacteria
i0 = 0 # default starting amount of infected bacteria
v0 = 1e5 # default starting amount of phages
tm = 2000.0 # default time span
end
## starting values ODE
begin
u0 = [s0, i0, v0]
tspan = [0.0, tm]
parms = [mu, kappa, omega, delta, eta, beta, lambda]
prob = ODEProblem(odeSolver!, u0, tspan, parms)
soln = solve(prob, Rosenbrock23())
end
## set obervables
p_mod = Observable(parms)
u_mod = Observable("")
## plot
# instantiate figure
fig = Figure(size = (1920, 1080))
ax = Axis(fig[1, 1:2])
# ranges
S_values = LinRange(0.0, 1.0e7, 100)
I_values = LinRange(0.0, 1.0e7, 100)
V_values = LinRange(0.0, 1.0e7, 100)
mu_values = LinRange(0.0, 1.0, 100)
kappa_values = LinRange(0.0, 1.0e9, 100)
omega_values = LinRange(0.0, 1.0, 100)
# labels
mu_val = @lift string("mu ", round($p_mod[1], digits=2))
k_val = @lift string("kappa ", round($p_mod[2], digits=2))
o_val = @lift string("omega ", round($p_mod[3], digits=2))
Label(fig[2,1], mu_val, fontsize = 18)
Label(fig[3,1], k_val, fontsize = 18)
Label(fig[4,1], o_val, fontsize = 18)
# sliders
mu_sld = Slider(fig[2, 2], range = mu_values, startvalue = mu)
k_sld = Slider(fig[3, 2], range = kappa_values, startvalue = kappa)
o_sld = Slider(fig[4, 2], range = omega_values, startvalue = omega)
# update parameters
on(mu_sld.value) do val
p_mod.val[1] = val
p_mod[] = p_mod[]
ylims!(ax, 0, max(maximum(X[]), maximum(Y[]))*1.1)
end
on(k_sld.value) do val
p_mod.val[2] = val
p_mod[] = p_mod[]
ylims!(ax, 0, max(maximum(X[]), maximum(Y[]))*1.1)
end
on(o_sld.value) do val
p_mod.val[3] = val
p_mod[] = p_mod[]
ylims!(ax, 0, max(maximum(X[]), maximum(Y[]))*1.1)
end
R = LinRange(tspan[1], tspan[2], 500)
data = lift(p_mod) do p_new
prob = ODEProblem(odeSolver!, u0, tspan, parms)
return solve(prob, Rosenbrock23(), p = p_new, saveat = R)
end
# draw
X = @lift $data[1,:]
Y = @lift $data[2,:]
lines!(ax, R, X, linewidth = 5, color =:blue)
lines!(ax, R, Y, linewidth = 2.5, color =:gold)