Unable to reproduce a simple bifurcation diagram

[vembha]

Hello,

New to BifurcationKit.jl here! I have a simple ODE-based model from an article (Fig. 2A in https://www.pnas.org/doi/10.1073/pnas.1419162112), which is known to have three bifurcation points when varied over a parameter m. The figure is attached:

The parameter is varied between 0 and 1, the dashed branches are unstable, and the solid branches are stable.

However, I am unable to reproduce the bifurcation diagram. The PALC() algorithm fails to identify the bifurcation points. It does not give out any errors. The code is below:

using BifurcationKit, Plots

# ODE system
function ConwayPerelson_2015(u, p)
    (λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ) = p
    T, L, I, V, E = u # The variable V is plotted in the bifurcation diagram
    [
        λ - d*T - β*V*T,
        αL*β*V*T + (ρ - a - dL)*L,
        (1 - αL)*β*V*T - δ*I + a*L - m*E*I,
        pr*I - c*V,
        λE + bE*E*I/(KB + I) - dE*E*I/(KD + I) - μ*E
    ]
end

# Parameters
params = (
    λ = 1E4, d = 0.01, β = 1.5E-8, δ = 1, pr = 2_000, c = 23,
    a = 0.001, dL = 0.004, ρ = 0.0045, αL = 1E-6,
    λE = 1, bE = 1, KB = 0.1, dE = 2, KD = 5, μ = 2, m = 0.42
)

# Initial guess
u0 = zeros(5)

# Bifurcation problem
prob = BifurcationProblem(ConwayPerelson_2015, u0, params, (@optic _.m);
record_from_solution = (x, p; k...) -> (V = x[4]))
opts = ContinuationPar(p_min = 0.0, p_max = 1.0)

# Solve
br = continuation(prob, PALC(), opts, bothside = true)

Can someone please help me out? Is there a way to ensure that the algorithm does not miss out on identifying the bifurcation points?

[rveltz]

You made a tiny mistake (semicolon missing) in the vector field

 (;λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ) = p

You should also make the computation of bifurcation points more precise opts = ContinuationPar(p_min = 0.0, p_max = 1.0, n_inversion=6)

[vembha]

Thank you very much for the prompt reply, and I am sorry for the silly miss of the semicolon.

The algorithm works, but it is giving me only one of the bifurcation points.

FYI, the estimated bifurcation point is highlighted as red II in the figure in the question. But there are two more. This is despite modifying the opts to detect fold bifurcations.

opts = ContinuationPar(p_min = 0.0, p_max = 1.0, n_inversion = 6,
        detect_bifurcation = 3, detect_fold = true, dsmin = 1E-4)
br = continuation(prob, PALC(), opts, bothside = true, normC = norminf)

Any idea how to ensure the algorithm finds all of them?

[rveltz]

diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max_steps = 10000))

gives the warning

┌ Error: Failure to converge with given tolerance = 1.0e-12.
│ Step = 9741
│  You can decrease the tolerance or pass a different norm using the argument `normC`.
│  We reached the smallest value [dsmin] valid for ds, namely 0.0001.
│  Stopping continuation at continuation step 9741.

but plot(diagram) yields:

The warning basically says that the newton tolerance are too high, hence:

@reset opts.newton_options.tol=1e-10
diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max_steps = 150000))
plot(diagram)

or:

plot(diagram, applytoY = x->log(1+abs(x)))

[vembha]

Thank you very much! I now see where I fell short, and the package is super cool . I was able to recapitulate the figure:

Two quick questions, however:

1. A non-existent hopf point is being predicted, very close to the point where a simple fold bifurcation is known to exist (see the StackTrace below). So, I am guessing it is an error of the algorithm. Are there any settings that I could change to avoid such incorrect estimations?

[Bifurcation diagram]
 ┌─ From 1-th bifurcation point.
 ├─ Children number: 0
 └─ Root (recursion level 2)
      ┌─ Curve type: EquilibriumCont from Transcritical bifurcation point.
      ├─ Number of points: 150001
      ├─ Type of vectors: Vector{Float64}
      ├─ Parameter m starts at 0.6087198085395702, ends at 0.22451213423608637
      ├─ Algo: PALC
      └─ Special points:

- #  1,       bp at m ≈ +0.40983773 ∈ (+0.40983773, +0.40983773), |δp|=2e-11, [converged], δ = ( 1,  0), step =  53
- #  2,       bp at m ≈ +0.83909847 ∈ (+0.83909847, +0.83909847), |δp|=2e-15, [    guess], δ = ( 1,  0), step = 7161
- #  3,     hopf at m ≈ +0.83853493 ∈ (+0.83853493, +0.83853493), |δp|=2e-09, [converged], δ = (-2, -2), step = 8138
- #  4, endpoint at m ≈ +0.22450773,                                                                     step = 150001

2. Is there a way to force the top stable branch to go all the way to m=0.0? I figured increasing max_steps does that, but it seems I would need a large value for it, probably because the top branch is being explored very slowly.

[rveltz]

@vembha

Is there a way to force the top stable branch to go all the way to m=0.0? I figured increasing max_steps does that, but it seems I would need a large value for it, probably because the top branch is being explored very slowly.

You can use a larger dsmax:

diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max_steps = 150000, dsmax = 0.7))

The issue here is that the range of values of ConwayPerelson_2015 is pretty large. I would be better numerically to rescale the model (log or other).

A non-existent hopf point is being predicted, very close to the point where a simple fold bifurcation is known to exist (see the StackTrace below).

You can see

julia> diagram[2]
[Bifurcation diagram]
 ┌─ From 1-th bifurcation point.
 ├─ Children number: 0
 └─ Root (recursion level 2)
      ┌─ Curve type: EquilibriumCont from Transcritical bifurcation point.
      ├─ Number of points: 150001
      ├─ Type of vectors: Vector{Float64}
      ├─ Parameter m starts at 0.6087056695593298, ends at 0.7526129563574361
      ├─ Algo: MoorePenrose
      └─ Special points:

- #  1,       bp at m ≈ +0.40983773 ∈ (+0.40983773, +0.40983773), |δp|=5e-12, [converged], δ = ( 1,  0), step = 225
- #  2,       bp at m ≈ +0.83909847 ∈ (+0.83909847, +0.83909847), |δp|=9e-16, [converged], δ = ( 1,  0), step = 35762
- #  3,     hopf at m ≈ +0.83853493 ∈ (+0.83853493, +0.83853493), |δp|=2e-09, [converged], δ = (-2, -2), step = 40644
- #  4, endpoint at m ≈ +0.75261202,  

It seems the Hopf bifurcation is not a fluke. See

julia> eigenvals(diagram[2].γ, 40644)
5-element Vector{ComplexF64}:
 -2.0506719532705528e-10 - 0.005389356755654721im
 -2.0506719532705528e-10 + 0.005389356755654721im
   -0.000499605123267488 + 0.0im
     -2.8970614512243027 + 0.0im
     -24.290528399220722 + 0.0im

julia> eigenvals(diagram[2].γ, 40645)
5-element Vector{ComplexF64}:
 -1.2226521003841756e-7 - 0.005389880502352536im
 -1.2226521003841756e-7 + 0.005389880502352536im
 -0.0004996051977692838 + 0.0im
    -2.8970638500412025 + 0.0im
     -24.29052805667206 + 0.0im

julia> eigenvals(diagram[2].γ, 40643)
5-element Vector{ComplexF64}:
  8.657536849168388e-8 - 0.005388984365665424im
  8.657536849168388e-8 + 0.005388984365665424im
 -0.000499605070280656 + 0.0im
   -2.8970597457475535 + 0.0im
    -24.29052864275146 + 0.0im

[rveltz]

You could use MTK to have analytical jacobian if you are unsure

[vembha]

Yes. Thank you!

However, despite using MTK for using the symbolic jacobian, the hopf is being predicted. Need to understand why … Thanks again for the help!

@variables T(t) L(t) I(t) V(t) E(t)
@parameters λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ
eqs = [
    D(T) ~ λ - d*T - β*V*T,
    D(L) ~ αL*β*V*T + (ρ - a - dL)*L,
    D(I) ~ (1 - αL)*β*V*T - δ*I + a*L - m*E*I,
    D(V) ~ pr*I - c*V,
    D(E) ~ λE + bE*E*I/(KB + I) - dE*E*I/(KD + I) - μ*E
]
@mtkbuild odesys = ODESystem(eqs, t)

bif_par = m
plot_var = V
params = [
    λ => 1E4, d => 0.01, β => 1.5E-8, δ => 1, pr => 2_000, c => 23,
    a => 0.001, dL => 0.004, ρ => 0.0045, αL => 1E-6,
    λE => 1, bE => 1, KB => 0.1, dE => 2, KD => 5, μ => 2, m => 0.0
]
u0 = zeros(5)

prob = BifurcationProblem(odesys, u0, params, bif_par;
        plot_var = plot_var, jac = true)

[rveltz]

Cool, maybe there is a small paper there if they did not see that.