# Fix incongruencies in ODE models
**URL:** https://discourse.julialang.org/t/fix-incongruencies-in-ode-models/72880
**Category:** Numerics
**Tags:** question
**Created:** [December 10, 2021, 7:45am UTC](https://discourse.julialang.org/t/fix-incongruencies-in-ode-models/72880 "2021-12-10T07:45:07Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 10, 2021, 7:45am UTC](https://discourse.julialang.org/t/fix-incongruencies-in-ode-models/72880/1 "2021-12-10T07:45:07Z")
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:
```julia
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,
```julia
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
---
### Author:  [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)
#### Post date: [December 10, 2021, 8:13am UTC](https://discourse.julialang.org/t/fix-incongruencies-in-ode-models/72880/2 "2021-12-10T08:13:53Z")
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[7]`, `u[8]`)? 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.
---
### Author:  [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)
#### Post date: [December 10, 2021, 9:00am UTC](https://discourse.julialang.org/t/fix-incongruencies-in-ode-models/72880/3 "2021-12-10T09:00:36Z")
Thanks,
if the call of the function is OK, then the problem is with the function itself. Actually, u[4] ends up negative, iota is always 0. This is part of the output of t, u[4], u[5], u[7], and u[8] (at about t=59 the system crashes):
Also, there was a typo on du[7]; it should be
```julia
du[7] = (ψ * u[6] * u[8]) - (Η * u[7]) - (ω * u[7])
```
however, u[8] starts lower than u[4] even if the values are the same:
