Continuation of Equilibria with Bifurcation Kit - model specific issue or wrong initial conditions?

[AP9]

I am using BifurcationKit.jl to analyze the steady states of my system varying a parameter A. In the figure, the bifurcation diagram predicts a stable steady state before the first Hopf bifurcation. However, when I plot the ODE solution for A inside the flat bit before the Hopf Bifurcation it doesn’t converge to the steady state predicted. I chose as initial conditions the value of the steady state identified by the continuation algorithm and I got the solutions in the first column in the Figure, for different A. Note that the solution is never converging but if I consider instead u0 = 0 (second column), I see convergence at least for the steep bits.

I have two questions:

1. Why by choosing the predicted stable steady states as initial conditions I don’t see convergence?

2. How can I choose the initial conditions to see convergence in the flat section before the hopf bifurcation? Is there something about my model that does not behave well in BifurcationKit? Since I am only facing this problem with this model, I am attaching the code for it, even if it is not a minimal working example.

using Plots, BifurcationKit, Parameters, DifferentialEquations
const BK = BifurcationKit
const DE = DifferentialEquations

@with_kw mutable struct Params  
    α :: Float64 = 0.01
    g::Float64 = 2.
    vθ::Float64 = -35.
    vL::Float64 = -65.
    C::Float64 = 0.1 
    gL::Float64 = 0.035
    v_reset::Float64 = -50.
    vh_TC::Float64 = -70.
    vh_RE::Float64 = -60.
    vT::Float64 = 120.
    τh_m::Float64 = 20.
    τh_p::Float64 = 100.
    gT::Float64 = 0.07
    vu_GABA::Float64 = -100.
    vu_AMPA::Float64 = 0.
    τR::Float64 = 5.
    A::Float64 = 0.
end

Heaviside(x) = x >= 0 ? 1.0 : 0.0 
Θ(x) = exp(- 0.001 / (x^2)) * Heaviside(x)   #infinitiely differentiable function 
 
function FiringRate(v, p)

    @unpack α, g, vθ, vL, C, gL, v_reset, vh_TC, vh_RE, vT, τh_m, τh_p, gT, vu_GABA, vu_AMPA, τR, A = p
   
    return Θ(v - vθ) / (τR + C / gL * log(abs((v - v_reset) / (v - vθ))))

end

function System!(du, u, p, t = 0)
    
    @unpack  α, g, vθ, vL, C, gL, v_reset, vh_TC, vh_RE, vT, τh_m, τh_p, gT, vu_GABA, vu_AMPA, τR, A = p

    v_TC, v_RE, h_TC, h_RE, x_TC, x_RE, y_TC, y_RE = u
    
    W_TC_RE = 1
    W_RE_TC = 1
    IL_TC = gL * (v_TC - vL)
    IT_TC = gT * h_TC * (v_TC - vT) * Θ(v_TC - vh_TC)
    τh_TC = τh_m + (τh_p - τh_m) * Θ(vh_TC - v_TC)
    h∞_TC = Θ(vh_TC - v_TC)
    I_TC = x_TC * (v_TC - vu_GABA)
    IL_RE = gL * (v_RE - vL)
    IT_RE = gT * h_RE * (v_RE - vT) * Θ(v_RE - vh_RE)
    I_RE = x_RE * (v_RE - vu_AMPA)
    τh_RE = τh_m + (τh_p - τh_m) * Θ(vh_RE - v_RE)
    h∞_RE = Θ(vh_RE - v_RE)
   
    I_drive = A
    
    du[1] = dv_TC = (- IL_TC - IT_TC - I_TC + I_drive) / C  
    du[2] = dv_RE = (- IL_RE - IT_RE - I_RE ) / C 
    du[3] = dh_TC = (-h_TC + h∞_TC) / τh_TC
    du[4] = dh_RE = (-h_RE + h∞_RE) / τh_RE
    du[5] = dx_TC = α * (y_TC - x_TC)
    du[6] = dx_RE = α * (y_RE - x_RE)
    du[7] = dy_TC = α * (g * W_TC_RE * FiringRate.(v_RE, Ref(p)) - y_TC) 
    du[8] = dy_RE = α * (g * W_RE_TC * FiringRate.(v_TC, Ref(p)) - y_RE)
    du
   
end

Continuation of equilibria

par = Params() 
u0 = zeros(8)
prob = BifurcationProblem(System!, u0, par, (@optic _.A),
    record_from_solution = (x, p; k...) -> (v_TC = x[1], v_RE = x[2], h_TC = x[3], h_RE = x[4], x_TC = x[5], x_RE = x[6], y_TC = x[7], y_RE = x[8]),)
opts_br = ContinuationPar(p_min = 0., p_max = 30.)
br = BK.continuation(prob, PALC(tangent=Bordered()), ContinuationPar(p_min = 0., p_max = 30.); normC = norminf)
plot(br)

Solution starting at stable steady state on the bifurcation diagram

A = 5 #choose 
idx = argmin(abs.(br.param .- A))
sol = br.sol[idx]
u0 = sol[1]
#u0 = zeros(8)
par.A = A
prob = ODEProblem(System!, u0, (0., 3000.), par)
sol = DE.solve(prob)
plot(sol, idxs =(0,1), xaxis = "t", yaxis = "v_TC", title = "A = $A", legend = false)

[rveltz]

If you run your code, you see that

julia> br
 ├─ Number of points: 323
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter A starts at 0.0, ends at 30.0
 ├─ Algo: PALC
 └─ Special points:

- #  1,     hopf at A ≈ +9.15409514 ∈ (+1.05106583, +9.15409514), |δp|=8e+00, [    guess], δ = ( 2,  2), step = 128
- #  2,     hopf at A ≈ +20.92989429 ∈ (+20.92982086, +20.92989429), |δp|=7e-05, [converged], δ = (-2, -2), step = 231
- #  3, endpoint at A ≈ +30.00000000,                                                                     step = 322

The first hopf point is very strange, it belongs to the interval (1, 9)!!

If you plot the branch like this

scatter(br)

you observe the gap:

So something is wrong with the continuation. It is probably because of Θ(x) having all its derivative being 0 at 0. Anyway, you can force a maximal step size using the callback cbMaxNormAndΔp (see the help).

Using this, one gets:

opts_br = ContinuationPar(p_min = 0., p_max = 30., n_inversion = 8, dsmax = 5e-1)
br = BK.continuation(prob, 
        PALC(tangent=Bordered()), plot = true,
        opts_br; 
    callback_newton = BK.cbMaxNormAndΔp(31.1,0.1), 
    )

with

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 401
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter A starts at 0.0, ends at 24.46540196006549
 ├─ Algo: PALC
 └─ Special points:

- #  1,     hopf at A ≈ +1.05106692 ∈ (+1.05106692, +1.05106692), |δp|=9e-12, [converged], δ = ( 2,  2), step =  63
- #  2,     hopf at A ≈ +20.92991650 ∈ (+20.92982974, +20.92991650), |δp|=9e-05, [converged], δ = (-2, -2), step = 351
- #  3, endpoint at A ≈ +24.55396009,                                                                     step = 401

[rveltz]

And just before the first Hopf point, it converges to the steady state

A = 1.04 
sol = br.sol[60]
u0 = sol.x
par.A = 1.04
prob = ODEProblem(System!, u0, (0., 3000.), par)
sol = DE.solve(prob)
plot(sol, xaxis = "t", title = "A = $A", legend = false)