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:

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