How to bring a system of differential equations to zero with DifferentialEquations?

Specific Domains Biology, Health, and Medicine
1

I have a population model described by this system of differential equations:

dN/dt = μN(1-(N+I)/κ) -φNP -ωN
dI/dt = φNP -ηI -ωI
dP/dt = βηI -φNP -ωP

with

κ = 2.2 ⋅ 10^7
μ = 0.47
φ = 10^-9
β = 50
η = 1.0
ω = 0.05

This system reaches equilibrium after a phase of oscillation. I would like to find the initial value of P that brings N to zero.
To find it, I have equated dN/dt to zero and solved for P, obtaining:
(μ/φ) * (1-(N+I/κ)) - ω/φ
N+I is the initial number of susceptible individuals because the infection simply shifts individuals from one compartment to another. I tried it out with this system where I have introduced an extinction clause (if N goes below one, I reset everything):

function 𝔅!(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
# extra parms
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(𝔅!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])

The susceptible N does indeed go down after the infection at t=1000, but does not go extinct:

treatment_1
treatment_11920×1440 117 KB

How can I find the value of P (in terms of concentration and time) that extinguish N? Is there some procedure in DifferentialEquations?
Thank you

1 Like
2

If I am not missing something, this can be solved as a Boundary Value Problem:
https://diffeq.sciml.ai/stable/types/bvp_types/

2 Likes
3

If you solve for dN/dt = 0, that only gives you an equilibrium for N, but not necessarily driving N\rightarrow 0. It sounds like you want to find a parameter (e.g. initial conditions, or one of the other constants) that will drive N to zero. I do agree with @AndiMD’s suggestion, and would add DiffEqFlux.jl as a great set of tools for this.

Also, I didn’t look closely at your equations, but I would guess that the way to drive N (susceptible) down is to infect everyone, and that could be accomplished with large P(0). If so, you could try simulating with increasing P(0), and in a few manual attempts you might find a suitable value without resorting to optimization.

4

It’s obvious by inspection that every value of P gives dN/dt = 0 when N = 0. I think you must have divided by N to get that solution. That doesn’t work when N = 0.

As a first step, I’d solve for lots of initial values of P, and plot the final value of N against the initial value of P. That will show you where to start looking. Bisection search might do the rest.

As a sanity check, solve for an initial value of P where N only just becomes zero, then another value where it only just fails to become zero, and see what’s different.

5

Thank you. Since I used 10^8 phages against 50 000 bacteria, all of them should have been infected. Instead, some bacteria shifted the system back. Also, when N=0%, I=100% so N+I is never zero; but even if we get N+I=0, then the logistic term becomes 1; moreover, N+I/K=50000/10^7, which cannot be 1 (perhaps I should increase N+I to 10^6-10^7 to have a more variable logistic term; like this is always about 0…).
The only parameters I can change are the input of phages (and bacteria) and the time of administration, for all the others are inherent biological factors. I would like to use an optimization method to obtain an objective and fast approach to determine these values for any given system. I’ll have a look into BVP. Thank you all!

6

You can use a callback to force the extinction when it gets small enough.

7

I think I did it with cb2 and cb3…

1 Like
8

So you’re saying, the callback triggered, it went to zero, but it came back?

9

In the other simulations I have done, no: once it went down it stayed there:

Scenario_8
Scenario_81920×1440 165 KB

With this simulation, the problem is that the blue line goes down too fast and the blue goes extinct. In the current problem, the blue line does not reach zero… I need to find the right time and amount to dispense the phage.

10

How would I set a Boundary Value Problem? What would be the main function? (du[1], I guess) And what the implicit function? (du[3], I guess). Thank you

11

EDIT2: see Chris’ comment for the correct function bc!(…)
EDIT: The following does not work, maybe the BVP solver is not the best way to solve this or I dont use it correctly.
I am not 100% sure the following is correct, but you can try this:

using BoundaryValueDiffEq
using OrdinaryDiffEq

function modelFun!(du, u, p, t) # base infection model with logistic term
    κ = 2.2E7
    μ = 0.47
    φ = 10^-9
    β = 50
    η = 1.0
    ω = 0.05
    #=
    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

"""Residual calculation. Assumes u[1] > 0"""
function bc!(residual, u, p, t)
    residual[1] = (u[1][1] - N0)/N0 # keep the starting value constant
    residual[2] = (u[1][2] - I0)/I0 # keep the starting value constant
    # residual[3] = (u[1][3] - P0)/P0 # keep the starting value constant
    residual[1] += u[end][1]/N0 # u[end][1] should be zero at the final timestep
end

N0 = 1.0*10^7
I0 = 1E5
P0 = 2*10^8

pGuess = 1.0
tspan = (0.0,1500.0)
bvp0 = BVProblem(modelFun!, bc!, [N0,I0,P0], tspan)
sol0 = solve(bvp0, Shooting(Vern7()), dtmin=0.01, dt=0.5)
12

Is the true solution actually zero at that point?

I just really don’t get the question. You can’t just say “and now the solution be zero!” and the solution shall becometh zero. When you define the problem, it has a solution. You can change the problem with callbacks, because that’s now a hybrid equation. But the question, if I’m reading it correctly, is how to bring a value to zero when the solution to the ODE does not have that value as zero. That requires that you define an intervention, because otherwise that is simply not a behavior of that actual solution. The solver gives you the solution to what you write down.

So I don’t think the BVP will necessarily work either, since there is possibly no true solution to the ODE that is truly zero at that endpoint. Your confusion with the bc! definition reflects this misunderstanding:

A BVP on 3 states only has 3 degrees of freedom. If you define 3 initial conditions, you’ve already defined the value at the end. You can’t change it to zero: see a standard elementary differential equation book on this. This is just mathematically not possible. You can ask the question where you say, the end value has to be zero, and thus you have to remove a degree of freedom elsewhere, which is essentially finding what value of one of the initial conditions will give you a zero at the end. But notice that what you want to define is 4 conditions on 3 equations, which is undetermined and thus there is no guarantee of a solution, and indeed in this case it looks like mathematically there is no solution. So a numerical method that solves the equations correctly will not work here, since there is no solution to find!

1 Like
13

As I understand it, we are looking for an initial value P0 to obtain N(end)==0. Assuming it can be 0 for the fixed starting values I0 and N0.

A BVP on 3 states only has 3 degrees of freedom. If you define 3 initial conditions, you’ve already defined the value at the end

Makes sense, I thought those were just initial parameters as well.

14

Oh okay, then the right BVP would be:

function bc!(residual, u, p, t)
    residual[1] = (u[1][1] - N0)/N0 # keep the starting value constant
    residual[2] = (u[1][2] - I0)/I0 # keep the starting value constant
    residual[3] = u[end][1]/N0 # u[end][1] should be zero at the final timestep
end
1 Like
15

Sorry, I am no mathematician, so most of this is beyond my comprehension. But I have found this – the problem is to distinguish the bacteria (blue line) by adding phage (green line) in what is known as phage therapy. I calculated the minimal amount to add by the formula Vp = (mu/phi) * (1-(r_s0/kappa)) - omega/phi (which I obtained by equating dN/dt to zero end solving for P). I don’t get extinction using such amount, but If I add more and more phage, I achieve extinction:

treatment_8
treatment_81920×1440 110 KB

treatment_1
treatment_11920×1440 113 KB

treatment_half9
treatment_half91920×1440 86.2 KB

So my question is: how do I calculate the right amount of phage and its time of delivery?
Perhaps the boundary problem is not the right tool, but the fact that if I add one billion phages I do reach extinction tells me there must be a solution. But Instead of randomly selecting values I would like to find the right amount straight away, even because the conditions will change in future models. Same thing for the time of administration. Here I subjectively chosen t=1000, but there should be a better time (specifically when bacteria are at a certain density, even if in this context the density is constant).
I tried the suggestions but:

julia> function bc!(residual, u, p, t)
           residual[1] = (u[1][1] - N0)/N0 # keep the starting value constant
           residual[2] = (u[1][2] - I0)/I0 # keep the starting value constant
           residual[3] = u[end][1]/N0 # u[end][1] should be zero at the final timestep
       end
bc! (generic function with 1 method)
julia> N0 = 1.0*10^7
1.0e7
julia> I0 = 1E5
100000.0
julia> P0 = 2*10^8
200000000
julia> pGuess = 1.0
1.0
julia> tspan = (0.0,1500.0)
(0.0, 1500.0)
julia> bvp0 = BVProblem(modelFun!, bc!, [N0,I0,P0], tspan)
BVProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 1500.0)
u0: [1.0e7, 100000.0, 2.0e8]
julia> sol0 = solve(bvp0, Shooting(Vern7()), dtmin=0.01, dt=0.5)
ERROR: MethodError: no method matching iterate(::DiffEqBase.NullParameters)
Closest candidates are:
  iterate(::Base.EnvDict) at env.jl:119
  iterate(::Base.EnvDict, ::Any) at env.jl:119
  iterate(::JuliaInterpreter.ExprSplitter) at /home/gigiux/.julia/packages/JuliaInterpreter/LBEYA/src/construct.jl:535
  ...
Stacktrace:
 [1] indexed_iterate(::DiffEqBase.NullParameters, ::Int64) at ./tuple.jl:84
 [2] 𝔅!(::Array{Float64,1}, ::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64) at /home/gigiux/Documents/model/util/funCollect.jl:7
 [3] (::ODEFunction{true,typeof(modelFun!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Float64,1}, ::Vararg{Any,N} where N) at /home/gigiux/.julia/packages/DiffEqBase/3iigH/src/diffeqfunction.jl:248
{etc}

Can I reach such a solution with boundary problems?
Thank you

16

Just do a DiffEqFlux style optimization.

1 Like