Hello,
I would like to find the minimal amount of phages required to distinguish a bacterial population. I have run this model:
using DifferentialEquations
function growth!(du, u, p, t)
μ, κ, φ, ω, η, β = p
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
=#
du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * 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])
end
# parameters
mu = 0.47 # maximum growth rate susceptible (B. longum)
kappa = 2.2*10^7 # maximum population density
phi = 10.0^-9 # adsorption rate resident commensal phage
eta = 1.0 # lyse rate resident commensal phage
beta = 50.0 # burst size resident commensal phage
omega = 0.05 # outflow
tmax = 4000.0 # time span 0-tmax
r_s0 = 1.0*10^7 # initial resident susceptible population u[1]
r_i0 = 0.0 # initial resident infected population u[2]
r_v0 = 1000.0 # initial resident phage population u[3]
T_inf = 1000 # time infection
amount_inf = 500 # amount infection
Vϕ = (mu/phi) * (1-(r_s0/kappa)) - omega/phi
Tϕ = T_inf
# apply
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, 0]
parms = [mu, kappa, phi, omega, eta, beta]
# inoculum
condition1(u, t, integrator) = t==T_inf
affect1!(integrator) = integrator.u[3] += amount_inf
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u[1]-1
function affect2!(integrator)
integrator.u[1] = 0
integrator.u[2] = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u[3]-1
function affect3!(integrator)
integrator.u[3] = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(growth!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[T_inf])
where I have a condition that simulates extinction when the number of elements falls below 1. I tried to calculate the best amount of phage with the formula Vϕ = (mu/phi) * (1-(r_s0/kappa)) - omega/phi which gives 2·10⁸ phages, but the value is not good because the bacterium does not go extinct:
Essentially, I would like to try several values for both
Vϕ (amount of phages) and Tϕ (time of inoculum), so Parallel Ensemble seems a good tool for the job. But how do I implement it?In the example reported on the SDE page, the Lotka-Volterra equation is modified by a function simulating noise:
function g(du, u, p, t)
du[1] = p[3] * u[1]
du[2] = p[4] * u[2]
end
but in my case, I need to modify Vϕ and Tϕ in growth’s u[3]. I also need to define the range of value to test with the additional function:
function prob_func(prob, I, repeat)
x = 0.3*rand(2)
remake(prob, p=[p[1:2]; x]
end
but what is the role of this function (that in the example simulates random numbers for function f implementing the Lotka-Volterra model.
What is the syntax in my case?
Let’s say I would like to test a range of Vϕ from 1 to 10⁹ with steps of 10 millions and Tϕ=tspan with steps of 1 hour. How would I do that?
Thank you
