Hello,

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[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
# 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[1]
r_i0 = 0.0 # initial infected population u[2]
r_v0 = 1000.0 # initial infectious agent population u[3]
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[3] += Vp
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(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_volterra`

function 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`loss`

function. How do I tell`DiffEqFlux`

to bring u[1] to zero using u[3]?

Thank you