Hello,
I am trying to simulate a bacteria/phage infection. I managed, thanks to the forum support, to make a simulation that in principle can change a parameter
in the ODEs – but here I am not doing such a change:
using DifferentialEquations
# function
function multiInfect!(du, u, p, t)
μ, ν, κ, φ, ψ, ω, β, Β, η, Η, ϛ = p
#=
mu = growth rate commensal
nu = growth rate opportunist
kappa = carrying capacity
phi = adsorption rate commensal
psi = adsorption rate opportunist
omega = flow rate
beta = burst size commensal
Beta = burst size opportunist
eta = lyse rate commensal
Eta = lyse rate opportunist
upstigma = reduction infectivity opportunist
du[1] = commensal sensitive
du[2] = commensal infected
du[3] = commensal phage
du[4] = opportunist sensitive
du[5] = opportunist infected
du[6] = opportunist phage
=#
sumBact = u[1]+u[2]+u[4]+u[5]
du[1] = (μ * u[1]) * (1 - sumBact/κ) - (φ * u[1] * u[3]) - (ω * u[1])
du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
du[4] = (ν * u[4]) * (1 - sumBact/κ) - ((ϛ*ψ) * u[4] * u[6]) - (ω * u[4])
du[5] = ((ϛ*ψ) * u[4] * u[6]) - (Η * u[5]) - (ω * u[5])
du[6] = (Β * Η * u[5]) - ((ϛ*ψ) * u[4] * u[6]) - (ω * u[6])
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
tmax = 500.0 # time span 0-tmax
t_in = 250.0
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
o_s0 = 1000.0 # initial opportunist bacteria
o_i0 = 0.0 # initial opportunist infected bacteria
o_v0 = 500.0 # initial opportunist phage
stigma = 1.0 # reduction of infection rate opportunist phage
# implementation
tspan = (0.0, tmax)
u0 = [c_s0, c_i0, c_v0, o_s0, o_i0, o_v0]
parms = [mu, nu, kappa, phi, psi, omega, beta, Beta, eta, Eta, stigma]
# modification
condition(u, t, integrator) = t==t_in # time of alteration
affect!(integrator) = integrator.p[11] = stigma # no change
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[t_in])
which gives me the following figure:
I would like to simulate now the activation of an additional phage. Thus,
function inductInfect!(du, u, p, t)
μ, ν, κ, φ, ψ, ω, β, Β, η, Η, Ψ = p
#=
mu = growth rate commensal
nu = growth rate opportunist
kappa = carrying capacity
phi = adsorption rate commensal
psi = adsorption rate opportunist
omega = flow rate
beta = burst size commensal
Beta = burst size opportunist
eta = lyse rate commensal
Eta = lyse rate opportunist
Psi = adsorption rate induced prophage
du[1] = commensal sensitive
du[2] = commensal infected with lytic phage
du[3] = commensal infected with induced prophage
du[4] = commensal phage
du[5] = commensal prophage
du[6] = opportunist sensitive
du[7] = opportunist infected
du[8] = opportunist phage
dN/dt = rN(1-N+I/K) - φNV - ωN
dI/dt = φNV - ηI - ωI
dV/dt = βηI - φNV - ωV
=#
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[4] * u[8]) - (Η * u[7]) - (ω * u[7])
du[8] = (Β * Η * u[7]) - (ψ * u[6] * u[8]) - (ω * u[8])
end
# additional parameters
tmax = 500
tspan = (0.0, tmax)
t_in = 250.0 # time of induction
p_in = 500.0 # amount of induction
c_s0 = 50000.0 # initial commensal sensitive
c_i0 = 0.0 # initial commensal infected with lytic phage
c_p0 = 0.0 # initial commensal infected with induced prophage
c_v0 = 500.0 # initial commensal phage
p_v0 = 0.0 # initial commensal prophage
o_s0 = 1000.0 # initial opportunist sensitive
o_i0 = 0.0 # initial opportunist infected
o_v0 = 500.0 # initial opportunist phage
Psi = 10.0^-9 # infection rate prophage
u0 = [c_s0, c_i0, c_p0, c_v0, p_v0, o_s0, o_i0, o_v0]
parms = [mu, nu, kappa, phi, psi, omega, beta, Beta, eta, Eta, Psi]
condition(u, t, integrator) = t==t_in
affect1!(integrator) = integrator.u[5] += 0
cb1 = DiscreteCallback(condition,affect1!)
affect2!(integrator) = integrator.u[3] += 0
cb2 = DiscreteCallback(condition,affect2!)
cb = CallbackSet(cb1, cb2)
prob = ODEProblem(inductInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[t_in])
# plot
The idea is that at t=t_in the second phage (u[5]) is added; since this alone did not work, I also added the infected bacteria (u[3]).
Even in this case, I am not adding a new virus but the figure does not look like the one before:
What am I getting wrong? the two figures should be the same since no modification has been introduced…
Thanks
