How to identify oscillatory conditions in Julia ODEs?

Luigi_Marongiu · May 28, 2022, 6:29am

Dear all,
I am simulating bacterial infections and I would like to know if there is a way to spot cases when the population enters an oscillatory equilibrium, for instance, as in this case:

To note that the two theoretical bacteria “susceptible” and “resistant” (to phage infection) oscillate around certain values (around 10^7 for the resistant, a bit more than 10^5 for the susceptible).

Is it possible to have a test that defines whether the ODEs describing the consortium hit this kind of equilibrium?

The equilibrium might of course change. For instance, by simply changing the infectivity of the phage I get:

But there are also instances where this equilibrium is not reached at all, for instance:

Thank you

hamblingreg52 · May 28, 2022, 8:52am

@greg_plowman Please share the solution to this problem if possible.

Datseris · May 28, 2022, 9:22am

Yes, use AttractorsViaRecurrences from DynamicalSystems.jl to get the attractor the system ends up. Let mapper = AttractorsViaRecurrences(...). Then, id = mapper(u0) with u0 the initial condition of interst. Get A = mapper.attractors[id] the attractor the initial condition ends up, and then x = A[:, 1] is the timeseries of the first variable of your system. You can use std(x) > e, with e some threshold value like 1e-3, to distinguish periodic motion from fixed point. To distinguish periodic motion from chaotic motion you can use lyapunov_from_data(A) from DynamicalSystems.jl.

To my knowledge this is the method that imposes the least assumptions about the nature of your system. You can come up with faster methods, but that enforce several assumptions about your system on the downside.

Refs:

vettert · May 28, 2022, 11:30am

You could also consider classical stability theory for ODEs, i.e. calculate eigenvalues of your coefficient matrix at steady state (you don’t even need to solve the dynamics of the system for this); then check whether they have non-zero imaginary parts.

xtalax · May 28, 2022, 11:35am

I’d champion spectral analysis here, take a windowed fourier transform and high pass and integrate the result. If the integrated high pass energy is greater than some fraction of the total energy, you have an oscillatory state.

stevengj · May 28, 2022, 12:10pm

That defines stable fixed points, but not limit cycles. To find periodic solutions of nonlinear ODEs, I’ve seen references to something called a harmonic balance method.

Datseris · May 28, 2022, 1:45pm

That would not be able to distinguish convergence to a periodic solution, or converging to a fixed point solution via an oscillating manner, unless you integrate for a very long time, so long, that the “no-motion” state of the system takes overwhelmingly more time than transitioning to said state. But if this is so, you might as well simply compare the distance to previous state in state space until distance falls below 1e-8 or so.

rveltz · May 28, 2022, 8:51pm

you can solve this with bifurcation theory. Take a look at the tutorials of BifurcationKit.jl and the tutorials to find periodic orbits near Hopf bifurcations

Luigi_Marongiu · June 27, 2022, 8:32pm

Dear Datseris,
thank you for the reply and apologies for the delay in answering. AttractorsViaRecurrences seems a very intriguing approach. Alas, I am a biologist, so I understood about nothing of it. The manual reports a lot of applications, but I don’t understand of to apply them to my case. In the article, I understand that there is a first part where one defines a solving function and then a second part of the actual implementation. Thus:

# Part 1
function lorenz84(u, p, t)
    F, G, a, b = p
    x, y, z = u
    dx = -y^2 -z^2 -a*x + a*F
    dy = x*y - y - b*x*z + G
    dz = b*x*y + x*z - z
    return SVector(dx, dy, dz)
end

# Part 2
xg = range(-1, 3; length = 100)
yg = range(-2, 3; length = 100)
zg = range(-2, 2.5; length = 100)
grid = (xg, yg, zg)
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)
basins, attractors = basins_of_attraction(grid, lo; diffeq)

Is the lorenz84 the function under investigation? Thus, would it be based on the ODEs I used to build the figure? Or is it based on a tridimensional space needed to solve the AttractorsViaRecurrences problem?
Thanks

rveltz · June 28, 2022, 3:55pm

Did you try BifurcationKit.jl? It seems that by changing the infectivity of the phage, you can get a Hopf bifurcation from which you can compute the periodic orbits. If you post your equations, il will be easier to help you

Luigi_Marongiu · June 28, 2022, 4:37pm

I kind of forgot the values I used for the figure (the file got lost in the sync), but the point is exactly: how to find the parameters that make an oscillatory environment. I tried with:

function logistic!(du, u, p, t)
    μ = p[:μ]   # susceptible growth rate
    ν = p[:ν]   # resistant growth rate
    κ = p[:κ]   # bacterial carrying capacity
    ω = p[:ω]   # system wash-out rate
    δ = p[:δ]   # phagial infection rate
    η = p[:η]   # phagial lysis rate (inverse latency)
    β = p[:β]   # phagial burst size
    λ = p[:λ]   # phagial decay rate
    #=
    du[1] = species A sensible
    du[2] = species A infected
    du[3] = species B
    du[4] = phages
    =#
    N = u[1] + u[2] + u[3] # total bacteria
    ρ = 1 - N/κ            # logistic term
    ϡ = (δ*u[1]*u[4])      # upsampi : infected bacteria
    du[1] = (μ*u[1]*ρ) - ϡ            - (ω*u[1])
    du[2] =              ϡ - (η*u[2]) - (ω*u[2])
    du[3] = (ν*u[3]*ρ)                - (ω*u[3])
    du[4] = (η*β*u[2]) - ϡ - (λ*u[4]) - (ω*u[4])
    return(SVector(du[1], du[2], du[3], du[4]))
end
begin
    kappa = 1e9        
    mu = 0.19
    nu = 0.16
    su0   = 1e6       
    re0   = 1e4     
    omega = 0.05            
    delta = 1e-9   
    eta   = 0.25     
    lambda = 0.05        
    beta  = 50.0           
    v0    = 1e5               
    t_mx  = 1e3          
end

u0 = [su0, 0, re0, v0]  # S, I, R, V
tspan = (0.0, t_mx)
parms = (μ=mu, ν=nu, κ=kappa, δ=delta, ω=omega, η=eta, β=beta, λ=lambda)
prob = ODEProblem(logistic!, u0, tspan, parms)
soln = solve( prob, Rosenbrock23(), tstops=[t_in] )

When I try the implementation:

julia> Si = range(su0/100, su0*100; length = 100)
10000.0:1.01e6:1.0e8
julia> Ii = range(0, 0; length = 100)
StepRangeLen(0.0, 0.0, 100)
julia> Ri = range(re0/100, re0*100; length = 100)
100.0:10100.0:1.0e6
julia> Vi = range(v0/100, v0*100; length = 100)
1000.0:101000.0:1.0e7
julia> grid = (Si, Ii, Ri, Vi)
(10000.0:1.01e6:1.0e8, StepRangeLen(0.0, 0.0, 100), 100.0:10100.0:1.0e6, 1000.0:101000.0:1.0e7)
julia> u0 = [su0, 0, re0, v0]
4-element Vector{Float64}:
      1.0e6
      0.0
  10000.0
 100000.0
julia> parms = (μ=mu, ν=nu, κ=kappa, δ=delta, ω=omega, η=eta, β=beta, λ=lambda)
(μ = 0.19, ν = 0.16, κ = 1.0e9, δ = 1.0e-9, ω = 0.05, η = 0.25, β = 50.0, λ = 0.05)
julia> lo = ContinuousDynamicalSystem(logistic, u0, parms)
ERROR: UndefVarError: logistic not defined
Stacktrace:
 [1] top-level scope
   @ none:1

Yet logistic is defined at the beginning. Why is it not recognized?

Luigi_Marongiu · June 28, 2022, 4:45pm

@ rveltz Thank you, but bifurcation would be just as obscure. I need a simple implementation that will tell how much bacteria and phages to put into a test tube to obtain a sustained (that is, the concentrations variates within some ranges) microbial consortium…

rveltz · June 28, 2022, 4:47pm

well it is a shame for you because the package answers the type of questions you are interested in, namely what is the quantitative dynamics as function of a parameter.

rveltz · June 28, 2022, 9:33pm

I suggest you do

using Parameters

function logistic!(du, u, p, t)
    @unpack μ, ν, κ, ω, δ, η, β, λ, = p
    # du[1] = species A sensible
    # du[2] = species A infected
    # du[3] = species B
    # du[4] = phages
    N = u[1] + u[2] + u[3] # total bacteria
    ρ = 1 - N/κ            # logistic term
    ϡ = (δ*u[1]*u[4])      # upsampi : infected bacteria
    du[1] = (μ*u[1]*ρ) - ϡ            - (ω*u[1])
    du[2] =              ϡ - (η*u[2]) - (ω*u[2])
    du[3] = (ν*u[3]*ρ)                - (ω*u[3])
    du[4] = (η*β*u[2]) - ϡ - (λ*u[4]) - (ω*u[4])
    return du
end

u0 = [1e6, 0, 1e4, 1e5]  # S, I, R, V
tspan = (0.0, 1e3)
parms = (μ=0.19, ν=0.16, κ=1e9, δ=1e-9, ω=0.05, η=0.25, β=50.0, λ=0.05)
prob = ODEProblem(logistic!, u0, tspan, parms)
soln = solve( prob, Rodas5P() )

what parameter do you chnage to see the oscillations?

Luigi_Marongiu · June 29, 2022, 5:25am

I’ll have a better look at bifurcations, then. For the model, the most impacting parameters are μ, ν, and δ. Thanks

Luigi_Marongiu · June 29, 2022, 5:44am

The problem here was to use logistic! instead of only logistic. I obtained:

julia> Si = range(1e6/100, 1e6*100; length = 10)
10000.0:1.111e7:1.0e8

julia> Ii = range(0, 0; length = 10)
StepRangeLen(0.0, 0.0, 10)

julia> Ri = range(1e4/100, 1e4*100; length = 10)
100.0:111100.0:1.0e6

julia> Vi = range(1e5/100, 1e5*100; length = 10)
1000.0:1.111e6:1.0e7

julia> grid = (Si, Ii, Ri, Vi)
(10000.0:1.111e7:1.0e8, StepRangeLen(0.0, 0.0, 10), 100.0:111100.0:1.0e6, 1000.0:1.111e6:1.0e7)

julia> u0 = [1e6, 0, 1e4, 1e5]
4-element Vector{Float64}:
      1.0e6
      0.0
  10000.0
 100000.0

julia> parms = (μ = 0.19, ν = 0.16, κ = 1.0e9, δ = 1.0e-9, ω = 0.05, η = 0.25, β = 50.0, λ = 0.05)

julia> lo = ContinuousDynamicalSystem(𝔉_logistic!, u0, parms)
4-dimensional continuous dynamical system
 state:       [1.0e6, 0.0, 10000.0, 100000.0]
 rule f:      logistic!
 in-place?    true
 jacobian:    ForwardDiff
 parameters:  (μ = 0.19, ν = 0.16, κ = 1.0e9, δ = 1.0e-9, ω = 0.05, η = 0.25, β = 50.0, λ = 0.05)

julia> diffeq = (alg = Rosenbrock23(), reltol = 1e-9, abstol = 1e-9)
(alg = Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}(nothing, OrdinaryDiffEq.DEFAULT_PRECS), reltol = 1.0e-9, abstol = 1.0e-9)

julia> basins, attractors = basins_of_attraction(grid, lo; diffeq)
┌ Warning: The function `basins_of_attraction(grid::Tuple, ds::DynamicalSystem; ...)` is
│ deprecated in favor of the more generic
│ `basins_of_attraction(mapper::AttractorMapper, grid::Tuple`) which works for
│ any instance of `AttractorMapper`. Please use that method in the future.
│ 
│ The only reason the existing method is kept in and can be used is because the paper
│ G. Datseris and A. Wagemakers, *Effortless estimation of basins of attraction*
│ had this old call signature published.
└ @ ChaosTools ~/.julia/packages/ChaosTools/PHPDF/src/basins/basins_of_attraction.jl:13
[ Info: Automatic Δt estimation yielded Δt = 18.046615917347687
Basins of attraction: 100%|████████████████████████████████████████████████████| Time: 0:01:31
julia> basins
10×10×10×10 Array{Int16, 4}:
[:, :, 1, 1] =
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
...
julia> attractors
Dict{Int16, Dataset{4, Float64}}()

julia> fracs = basins_fractions(basins)
Dict{Int16, Float64} with 1 entry:
  -1 => 1.0

julia> for (key, att) in attractors
           u0 = att[1] # First found point of attractor
           ls = lyapunovspectrum(lo, 10000; u0)
           println("Attractor $(key) has spectrum: $(ls) and fraction: $(fracs[key])")
       end

julia>

What is the new syntax for basins_of_attraction?
What is the practical use of basin and attractors?
Thanks

rveltz · June 29, 2022, 8:35am

what are the parameters for the first 3 figures?

rveltz · June 29, 2022, 8:46am

There many methods: collocation, shooting, finite differences, fourier series. The best for precision is orthogonal collocation ala AUTO i’d say.

Luigi_Marongiu · June 29, 2022, 8:04pm

I lost the file with the original conditions and I can’t find the parameters that replicate it. But I have a more generic model where each bacterium has a phage:

function Doppel!(du, u, p, t) 
    μ = p.μ   # bacterial growth rate
    κ = p.κ   # bacterial carrying capacity
    ω = p.ω   # system wash-out rate
    δ = p.δ   # phagial infection rate
    η = p.η   # phagial lysis rate (inverse latency)
    β = p.β   # phagial burst size
    λ = p.λ   # phagial decay rate
    #=
    du[1] = species A sensible
    du[2] = species A infected
    du[3] = phage a
    du[4] = species B sensible
    du[5] = species B infected
    du[6] = phage b
    =#
    N = u[1] + u[2] + u[4] + u[5]   # total bacteria
    ρ = 1 - N/κ                     # logistic term
    ϡ = (δ[1]*u[1]*u[3])            # upsampi: infected bacteria A
    ς = (δ[2]*u[4]*u[6])            # varsigma: infected bacteria B

    du[1] = (μ[1]*u[1]*ρ)       - ϡ                     - (ω*u[1])
    du[2] =                       ϡ     - (η[1]*u[2])   - (ω*u[2])
    du[3] = (η[1]*β[1]*u[2])    - ϡ     - (λ[1]*u[3])   - (ω*u[3])

    du[4] = (μ[2]*u[4]*ρ)       - ς                     - (ω*u[4])
    du[5] =                       ς     - (η[2]*u[5])   - (ω*u[5])
    du[6] = (η[2]*β[2]*u[4])    - ς     - (λ[2]*u[5])   - (ω*u[6])
end
# parameters
kappa   = 2e7               # carrying capacity (from Weitz)
omega   = 0.05              # outflow (from Weitz)
mu      = [0.16, 0.31]      # growth rate susceptible (pathogen)
beta    = [50.0, 75.0]      # burst size           (from Weitz)
delta   = [1e-9, 1e-9]      # adsorption rate (from Weitz)
eta     = [0.25, 0.75]      # reciprocal of latency   (from Weitz)
lambda  = [0, 0]            # decay rate in hours     (from Weitz)
As0     = 1e6               # initial susceptible density A
Ai0     = 0                 # initial infected density A
a0      = 1e5               # initial phage density a
Bs0     = 1e5               # initial susceptible density B
Bi0     = 0                 # initial infected density B
b0      = 1e4               # initial phage density b
t_mx  = 5e3
# implementation
u0 = [As0, Ai0, a0, Bs0, Bi0, b0]
tspan = (0.0, t_mx)
parms = ComponentArray(μ=mu, κ=kappa, ω=omega, δ=delta, η=eta, β=beta, λ=lambda)
prob = ODEProblem(Doppel!, u0, tspan, parms)
soln = solve(prob, Rosenbrock23())
# find steady state
eq = solve(prob, Rosenbrock23(), callback = TerminateSteadyState())

In this case, the carrying capacity makes quite a difference: if I use 30 000 000, I get:

eq gives the point of steady-state, as determined in a previous post.

stevengj · June 29, 2022, 10:31pm

These methods work best if you know the period, so that you can formulate the periodic solution as a boundary-value problem (BVP), no? For an unknown period, you would then have to wrap the BVP in some kind of root-finding algorithm, and the tricky part is to efficiently combine this with the BVP solver.

As I understand it, the “harmonic-balance method” is precisely a technique for combining a Fourier-series collocation BVP solver with root finding or numerical continuation.

There are probably other techniques as well, but in my mind this combination is the key question that distinguishes the problem of finding periodic ODE solutions from ordinary periodic boundary-value problems.