Failing to obtain a bifurcation plot

Specific Domains Modelling & Simulations
1

Hello everyone,

A long-time user of Julia/SciML. I have a system of ODEs and have been trying to obtain a bifurcation diagram using BifurcationKit.jl and ModelingToolkit.jl. While I am able to successfully obtain one branch, my code is failing to perform branch switching and explore other branches. It seems the eigenvalue estimated for the bifurcation point is far from zero.

The code is as follows:

#=
Section 1: Loading required pacakges and files
=#

using XLSX, CSV, DataFrames, Symbolics
using DifferentialEquations, Plots
using ModelingToolkit, BifurcationKit
using ModelingToolkit: t_nounits as t, D_nounits as D

#=
Section 2: Define the model system as per the requirement
=#

@variables T(t) I(t) N(t) P(t) E(t) M(t) R(t) A(t)
@parameters λ, β_exp, γ, fD, fL, a, μI, μL, μD, λN, rNP, αP, ϕP, θ, θN, ω, θC, αE, rEM, μE, ϕM, μX, αA, θA, T0, N0, μT, μN, m, n
# β_exp is varied for the bifurcation plot

# Dynamical equations
eqs = [
    D(T) ~ λ - (10^(β_exp))*T*I - μT*T,
	D(I) ~ (1 - fL - fD)*(10^(β_exp))*T*I - (E + R)*I - μI*I,
	D(N) ~ λN - rNP*N*I/(θN + I) - μN*N,
	D(E) ~ (ω*ϕP + αE*E)*((θC^m)/(θC^m + R^m))*I/(θ + I) - (rEM*(1 - M/ϕM) + μE)*E,
    D(M) ~ (rEM*E + rEM*R)*(1 - M/ϕM) - μX*M*(A/(1 + A))*(I/(θ + I)),
	D(R) ~ (ω*M + αE*R)*((θC^m)/(θC^m + E^m))*I/(θ + I) - (rEM*(1 - M/ϕM) + μE)*R,
    D(A) ~ αA*(I^n)/((θA)^n + I^n) - rNP*(N/ϕP)*(I/(θN + I))*A
]

# Parameters
bif_par = β_exp
plot_var = I
params = [
    λ => 352.7, β_exp => -7.0, γ => 1026.27, fD => 0.93, m => 4.0, n => 4.0,
    fL => 3.22E-7, a => 0.001, μI => 0.5, μL => 0.0033, μD => 0.069,
    λN => 10^(-4.3), rNP => 10^(0.15), αP => 0.63, ϕP => 0.5, θ => 1.0,
    θN => 0.01, ω => 10^(0.22), θC => 0.01, αE => 1.95, rEM => 10^(-0.5),
    μE => 2.0, ϕM => 1_000.0, μX => 1.2, αA => 0.0046, θA => 5.0,
    T0 => 10^(6.68), N0 => 10^(-0.3), μT => 7.37E-5, μN => 1E-4
]
@mtkbuild odesys = ODESystem(eqs, t)

# Initial guesses
u0 = ones(7); u0[1] = 10^(curr_p[26]); u0[3] = 10^(curr_p[28]);

#=
Section 3: Define and solve the bifurcation problem
=#

# Bifurcation problem
prob = BifurcationProblem(odesys, u0, params, bif_par; plot_var = plot_var, jac = true)
opts = ContinuationPar(p_min = -7.0, p_max = -5.0, n_inversion = 6, detect_bifurcation = 3, detect_fold = true, dsmax = 0.2)

# Solve
@reset opts.newton_options.tol = 1E-12
@reset opts.newton_options.max_iterations = 100
br = bifurcationdiagram(prob, PALC(), 2, ContinuationPar(opts, max_steps = 10_000), bothside = true, verbosity = 3)
# plot(br)

The following is the relevant portion of the error:

┌ Warning: The zero eigenvalue is not that small λ = 0.11281601151143743, this can alter the computation of the normal form. You can either refine the point using Newton or use a more precise bisection by increasing `n_inversion`
└ @ BifurcationKit C:\Users\vembh\.julia\packages\BifurcationKit\dbu6D\src\NormalForms.jl:82
├─ smallest eigenvalue at bifurcation = 0.11281601151143743
┌── left eigenvalues =
3-element Vector{ComplexF64}:
     0.11281601151143739 + 0.0im
   9.724595179899156e-21 + 0.0im
 -1.9227506899814033e-19 + 0.0im
├── right eigenvalue = 0.11281601151143743
└──  left eigenvalue = 0.11281601151143739 + 0.0im
┌ Warning: The bifurcating eigenvalue is not that close to Re = 0. We found 0.11281601151143739 !≈ 0.  You can perhaps increase the argument `nev`.
└ @ BifurcationKit C:\Users\vembh\.julia\packages\BifurcationKit\dbu6D\src\NormalForms.jl:36
┌── Normal form:   aδμ + b1⋅x⋅δμ + b2⋅x²/2 + b3⋅x³/6
├─── a    = -1.1804118293263817e-22
├─── b1   = 1.411061001076516
├─── b2/2 = 0.0018579145893656754
└─── b3/6 = 0.012470340030255723
...
┌ Error: Failed to compute new branch at p = -5.737703910147452
│   exception = LinearAlgebra.SingularException(8)
└ @ BifurcationKit C:\Users\vembh\.julia\packages\BifurcationKit\dbu6D\src\bifdiagram\BifurcationDiagram.jl:210

Can someone please help me out?