How to bring differential equations to wanted value with DiffEqFlux?

Luigi_Marongiu · January 29, 2021, 7:15am

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:

treatment_81920×1440 110 KB

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

ChrisRackauckas · January 29, 2021, 8:05am

What do you mean “using u[3]”? You can do a loss function to find what parameters cause u[1] to go to zero at a time T, but the question you’re asking doesn’t seem to make sense?

Luigi_Marongiu · January 29, 2021, 8:17am

If the parameters are μ, κ, φ, ω, η, β = p, then no, I am not interested in that. I would like to know how many phages = Vp = u[3] I need to add and at what time point Tp in b1. As you can see in the other post (of which the present post is, ideally, a more focused continuation) there is an amount of phage = u[3] that collapses u[1]. As you suggested, how do I train DiffEqFlux to use u[3] as such a parameter? and how do I tell DiffEqFlux to take into account also the time? Thanks

ChrisRackauckas · January 29, 2021, 8:21am

Make Vp a parameter.

parms = [mu, kappa, phi, omega, eta, beta, Vp]
affect1!(integrator) = integrator.u[3] += integrator.p[7]

soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp], sensealg=ForwardDiffSensitivity())

should be all you need.

Luigi_Marongiu · January 29, 2021, 8:22am

Thanks, I’ll try…

Luigi_Marongiu · January 30, 2021, 11:02am

Something weird happened. Until the day before yesterday, when I was introducing 200 million phages, the system went into a cyclic phase and then reached an equilibrium:

Yesterday I ran this code:

parms = [mu, kappa, phi, omega, eta, beta, Vp]
# modification
condition1(u, t, integrator) = t==Tϕ
affect1!(integrator) = integrator.u[3] += integrator.p[7]
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], sensealg=ForwardDiffSensitivity())

and, for the same amount I got the extinction:

Since the amount inserted is the same, I thought the solver itself has changed. So, I closed Julia’s core and started a new session but the results are the same.
What happened?
Can the solver be reversed as it was?
I think the question I have (found the amount of phage that extinguishes the bacterium) cannot be found mathematically, but from the other post (obtained with the solver as it was) it can be seen that there must be an amount that triggers such extinction.
So, I’ll rephrase the question: is there a computational savvy way to run the sir function in a loop providing different Vp?
Thanks

ChrisRackauckas · January 30, 2021, 11:35am

You can change versions at any time. Are you sure it’s a version thing?

https://diffeq.sciml.ai/stable/features/ensemble/

Luigi_Marongiu · January 30, 2021, 12:56pm

I don’t think is a solver’s version problem, for I did not install a new version of DifferentialEquations: I just added BoundaryValueDiffEq, Flux, Optim, DiffEqFlux, DiffEqSensitivity. But how can it be that now the solver brings the system to zero even when using a fresh session but before it did not? I assumed that the code I used changed the solver in a stable manner (if possible)…

ChrisRackauckas · January 30, 2021, 6:46pm

repro?

Luigi_Marongiu · January 31, 2021, 7:50am

Hello, I don’t know what “repro” means but I tried again. If I give the value Vp = 200 million phages, I now get extinction, whereas before that was not the case. Moreover, If I use up to Vp/21 I still get extinction but I don’t with amounts below Vp/22. Standing that I don’t understand why I was getting a different result before running DiffEqFlux, it is still evident that the minimal amount of phages is not Vp but another value that DiffEqFlux did not help finding. I’ll try the ensemble simulations you’ve suggested. Thanks

Luigi_Marongiu · February 10, 2021, 8:26pm

Hello, how do i set the code for the parallel ensemble simulation? I can see that:
EnsembleProblem(prob::DEProblem; output_func = (sol, I),-> (sol, false), prob_func = (prob, i, repeat) -> (prob), reduction = (u, data, I) -> (append!(u, data), false) where prob is the one I have set in the first post. Is this verbatim OK? and how do I tell what parameter to check and in what range? The example in the link gives the code function prob_func(prob, I, repeat) @. [...] but I don’t understand what it does.

ChrisRackauckas · February 11, 2021, 10:04am

You just say how to change the base problem to new ones. So, remake it with new parameters.

Luigi_Marongiu · February 23, 2021, 6:35am

Hello, could you please write an example of how should I set the function? I tried to focalize the problem with this post but I got no answer… Thank you

ChrisRackauckas · February 23, 2021, 8:45am

https://diffeq.sciml.ai/stable/features/ensemble/#Random-Initial-Conditions

Luigi_Marongiu · February 23, 2021, 10:32am

I have seen that but I don’t know how to use it…

ChrisRackauckas · February 23, 2021, 12:55pm

What’s the exact question here? There’s plenty of examples already, so without a simple exact question it’s hard to do much more. It’s literally just a function that spits out the problem to solve for trajectory i.

Luigi_Marongiu · February 26, 2021, 6:53am

I have been looking at example 2 of the vignette since it covers a Lotka-Volterra system. This example has a main function that resembles the one I called growth in my second post. The problem is that the example than sets another function:

function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
end

What is the use of this function in my case?
I just want to test over u[3] in growth. How do I write that down?
shall I do:

function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
  du[3] = p[5]*u[3]
end

or simply

function g(du,u,p,t)
  du[3] = p[1]*u[3]
end

What about the indices?

ChrisRackauckas · February 26, 2021, 10:15am

What do you even mean by this? What are you trying to do mathematically? Do you mean you want to sample over the parameters of the noise? Then yes, make there be a parameter p that is used in the diffusion equation.

That’s completely dependent on the model. Is that the same parameter that is used in the drift equation f? Then use the same index. Is it a different one? Then use a different index. It’s literally just the same p that you pass in. There’s nothing fancy there.

Luigi_Marongiu · February 26, 2021, 10:41am

I just need to pass different values to u[3] at the selected step. I got:

condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += i

I just want to test several combinations of Tp and i.

ChrisRackauckas · February 26, 2021, 12:00pm

condition1(u, t, integrator) = t==integrator.p[4]
affect1!(integrator) = integrator.u[3] += integrator.p[5]

and follow the examples.