I have a population system (bacteria/phages) that I have described with these ODEs that include the logistic factor:
function sir!(du, u, p, t) # SIRl with logistic term μ, κ, φ, ω, η, β = p #= du = susceptible du = infected du = phages =# du = ((μ * u) * (1 - ((u+u)/κ))) - (φ * u * u) - (ω * u) du = (φ * u * u) - (η * u) - (ω * u) du = (β * η * u) - (φ * u * u) - (ω * u) end # parms mu = 0.47 # maximum growth rate susceptible 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 tmax = 4000.0 # time span 0-tmax r_s0 = 50000.0 # initial susceptible population u r_i0 = 0.0 # initial infected population u r_v0 = 1000.0 # initial infectious agent population u Tp = 1000 # time of infection # execute Vp = (mu/phi) * (1-(r_s0/kappa)) - omega/phi tspan = (0.0, tmax) u0 = [r_s0, r_i0, 0] parms = [mu, kappa, phi, omega, eta, beta] # inoculum condition1(u, t, integrator) = t==Tp affect1!(integrator) = integrator.u += Vp cb1 = DiscreteCallback(condition1, affect1!) # extintion condition2(u, t, integrator) = u-1 function affect2!(integrator) integrator.u = 0 integrator.u = 0 end cb2 = ContinuousCallback(condition2, affect2!) condition3(u, t, integrator) = u-1 function affect3!(integrator) integrator.u = 0 end cb3 = ContinuousCallback(condition3, affect3!) # run modification = CallbackSet(cb1, cb2, cb3) prob = ODEProblem(SIR!, u0, tspan, parms) soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])
The model includes
- a time point Tp when inoculation of Vp phages is given
- a term that accounts for the extintion of a species by setting the equation to zero when the cell count goes below 1.
With this model, I get a system as:
My problem is that I would like to find the Vp and Tp that bring the blue line to zero.
I have looked into DiffEqFlux but I don’t know how to implement it.
In the example given in the link, I can change the
lotka_volterrafunction with my
SIR. I got the initial conditions
u0, the intermediary points
tspan, and the parameters
p. The problem is how do I set the
lossfunction. How do I tell
DiffEqFluxto bring u to zero using u?