How to reality-chek a working ODE model made with DifferentialEquations?

Specific Domains Modelling & Simulations
1

Hello,
I have made a simulation of bacterial infection. I am focusing on the induction of phages. The system works but I am not sure about the results. Even if I activate 400 000 hosts at t=250 (that is, 96.6% of the bacteria), there is very little effect on the population, as shown here:

Scenario_4
Scenario_41920×1440 229 KB

It might be that the system is inherently resilient, but I fear there might be simply a problem with the implementation of the ODEs.
Is there a way to discriminate between these two cases? How do I reality-check the system to avoid biased but yet working models?
Thank you.

The code is here:

using DifferentialEquations, PyPlot

function inductInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ, ω, β, Β, η, Η, Ψ = p
sumBact = (u[1]+u[2]+u[3])+(u[6]+u[7])
    ι = (φ * u[1] * u[4]) + (Ψ * u[1] * u[5]) # iota
    du[1] = (μ * u[1]) * (1 - sumBact/κ) - ι - (ω * u[1])
    du[2] = (φ * u[1] * u[4]) - (η * u[2]) - (ω * u[2])
    du[3] = (Ψ * u[1] * u[5]) - (η * u[3]) - (ω * u[3])
    du[4] = (β * η * u[2]) - (φ * u[1] * u[4]) - (ω * u[4])
    du[5] = (β * η * u[3]) - (Ψ * u[1] * u[5]) - (ω * u[5])

    du[6] = (ν * u[6]) * (1 - sumBact/κ) - (ψ * u[6] * u[8]) - (ω * u[6])
    du[7] = (ψ * u[6] * u[8]) - (Η * u[7]) - (ω * u[7])
    du[8] = (Β * Η * u[7]) - (ψ * u[6] * u[8]) - (ω * u[8])
end

# parameters
mu    = 0.47        # maximum growth rate Commensal (B. longum)
nu    = 0.72        # maximum growth rate opportunist (E. coli)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate phage Commensal
psi   = 10.0^-9     # adsorption rate phage opportunist
eta   = 1.0         # lyse rate phage Commensal
Eta   = 1.0         # lyse rate phage opportunist
beta  = 50.0        # burst size phage Commensal
Beta  = 50.0        # burst size phage opportunist
tau   = 3.62        # latency time
omega = 0.05        # outflow
c_s0    = 50000.0   # initial Commensal bacteria
c_i0    = 0.0       # initial Commensal infected bacteria
c_v0    = 500.0     # initial Commensal phage
c_p0    = 0.0       # initial commensal infected with induced prophage
c_vl0   = 500.0     # initial commensal lythic phage
c_vi0   = 0.0       # initial commensal induced prophageo_s0    = 1000.0    # initial opportunist bacteria
o_i0    = 0.0       # initial opportunist infected bacteria
o_v0    = 500.0     # initial opportunist phage
tmax  = 500
tspan = (0.0, tmax)
t_in = 250.0        # time of induction
p_in = 400000       # amount of induction (10% population)
Psi = 10.0^-9       # infection rate prophage
u0 = [c_s0, c_i0, c_p0, c_vl0, c_vi0, o_s0, o_i0, o_v0]
parms = [mu, nu, kappa, phi, psi, omega, beta, Beta, eta, Eta, Psi]
condition1(u, t, integrator) = t==t_in
affect1!(integrator) = integrator.u[5] += p_in
cb1 = DiscreteCallback(condition1, affect1!)
condition2(u, t, integrator) = t==t_in
affect2!(integrator) = integrator.u[3] += p_in
cb2 = DiscreteCallback(condition2, affect2!)
cb = CallbackSet(cb1, cb2)
prob = ODEProblem(inductInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[t_in])
2

This is a modeling question moreso than a package usage question. You have to think about what could be true about the model’s solution under different perturbations, try solving at different parameters, and confirm that these qualitative changes hold. You can get help from things like GlobalSensitivity.jl, DiffEqUncertainty.jl, etc. to quantify different aspects, but ultimately the modeler needs to understand the model in order to evaluate whether it has been correctly described.

1 Like