# Fix incongruencies in ODE models

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
omega = flow rate
beta  = burst size commensal
Beta  = burst size opportunist
eta   = lyse rate commensal
Eta   = lyse rate opportunist
upstigma = reduction infectivity opportunist

du = commensal sensitive
du = commensal infected
du = commensal phage

du = opportunist sensitive
du = opportunist infected
du = opportunist phage
=#
sumBact = u+u+u+u
du = (μ * u) * (1 - sumBact/κ) - (φ * u * u) - (ω * u)
du = (φ * u * u) - (η * u) - (ω * u)
du = (β * η * u) - (φ * u * u) - (ω * u)

du = (ν * u) * (1 - sumBact/κ) - ((ϛ*ψ) * u * u) - (ω * u)
du = ((ϛ*ψ) * u * u) - (Η * u) - (ω * u)
du = (Β * Η * u) - ((ϛ*ψ) * u * u) - (ω * u)
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 = 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
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 = commensal sensitive
du = commensal infected with lytic phage
du = commensal infected with induced prophage
du = commensal phage
du = commensal prophage

du = opportunist sensitive
du = opportunist infected
du = opportunist phage

dN/dt = rN(1-N+I/K) - φNV - ωN
dI/dt = φNV - ηI - ωI
dV/dt = βηI - φNV - ωV
=#
sumBact = (u+u+u)+(u+u)
ι = (φ * u * u) + (Ψ * u * u) # iota
du = (μ * u) * (1 - sumBact/κ) - ι - (ω * u)
du = (φ * u * u) - (η * u) - (ω * u)
du = (Ψ * u * u) - (η * u) - (ω * u)
du = (β * η * u) - (φ * u * u) - (ω * u)
du = (β * η * u) - (Ψ * u * u) - (ω * u)

du = (ν * u) * (1 - sumBact/κ) - (ψ * u * u) - (ω * u)
du = (ψ * u * u) - (Η * u) - (ω * u)
du = (Β * Η * u) - (ψ * u * u) - (ω * u)
end
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 += 0
cb1 = DiscreteCallback(condition,affect1!)
affect2!(integrator) = integrator.u += 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) is added; since this alone did not work, I also added the infected bacteria (u).
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

In order to make the process clear to us, can you add the new states to the end of the state vector(like make them `u`, `u`)? It will be easier to separate out new ones. It seems that, in the second one, opportunist phage` is not correct from the start if we consider the first plot ground truth.

Thanks,
if the call of the function is OK, then the problem is with the function itself. Actually, u ends up negative, iota is always 0. This is part of the output of t, u, u, u, and u (at about t=59 the system crashes):

Also, there was a typo on du; it should be

``````    du = (ψ * u * u) - (Η * u) - (ω * u)

``````

however, u starts lower than u even if the values are the same:

1 Like