Producing Fold Bifurcation Diagram using BirfurcationKit

I am trying to construct a bifurcation diagram in Julia using Bifurcation.Kit package. I cannot get the graph I am expect to produce. For some context, I am modeling a system of ODEs (non linear, 6 equations, and 6 unknowns). I know that there are two alternative stable states, and therefore, I should expect a fold bifurcation. However, that is not showing up in my graph. I could really use any code debugging, insight, directions, or contacts of people who might be able to help.

I will be happy to clarify any aspect.

using Revise, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit

#Supremum Norm
norminf = x -> norm(x, Inf)

const m    = 0.250 #Mortality P and J
      s    = 0.500 #Stocking P
      f    = 1.000 #Fecundity
      c    = 0.100 #cannibalism
      kₓ   = 0.250 #Attack rate F on J in littoral
      k_y  = 0.500 #Attack rate F on J in pelagic
      r    = 0.500 #Forage fish growth rate
      b    = 0.050 #Density dependent constant forage fish
      hₓ   = 0.250 #Attack rate P on F in littoral
      h_y  = 0.500 #Attack rate P on F in pelagic
      i    = 1.000 #Stocking F
      ϵ    = 0.250 #Trophic Efficiency
      γ    = 0.500 #Area Ratio
      σ₁   = 1.000 #Dispersal
      σ₂   = 1.000 #Dispersal
      σ₃   = 1.000 #Dispersal
      σ₄   = 1.000 #Dispersal
      σ₅   = 1.000 #Dispersal 
      σ₆   = 1.000 #Dispersal 
      tmat = 5.000 #Maturation

#Defining a vector field
function VF!(dz, z, p, t)
    p₁, p₂ = p                     #Parameters
    u₁, u₂, u₃, u₄, u₅, u₆     = z #State Variables
    #Predator Equations
    dz[1] = u₅/tmat - m*u₁ - p₁*u₁ + s + ϵ*hₓ*u₁*u₃              + σ₁*γ*p₂*u₂ - σ₂*(1/γ)*p₂*u₁    #Littoral
    dz[2] = u₆/tmat - m*u₂ - p₁*u₂ + s + ϵ*h_y*u₂*u₄             - σ₁*γ*p₂*u₂ + σ₂*(1/γ)*p₂*u₁    #Pelagic

    #Foragers Equations
    dz[3] = r*u₃ - b*u₃*u₃ - hₓ*u₁*u₃ +  i + kₓ*u₃*u₅ + σ₃*γ*p₂*u₄ - σ₄*(1/γ)*p₂*u₃             #Littoral
    dz[4] = r*u₄ - b*u₄*u₄ - h_y*u₂*u₄ + i + k_y*u₄*u₆ - σ₃*γ*p₂*u₄ + σ₄*(1/γ)*p₂*u₃            #Pelagic
    #Juveniles Equations
    dz[5] = f*(u₁+γ*u₂) - u₅/tmat - m*u₅ - c*u₁*u₅ - kₓ*u₃*u₅ + σ₅*γ*p₂*u₆ - σ₆*(1/γ)*p₂*u₅       #Littoral
    dz[6] = -u₆/tmat - m*u₆ - c*u₂*u₆ - k_y*u₄*u₆ - σ₅*γ*p₂*u₆ + σ₆*(1/γ)*p₂*u₅                   #Pelagic

#Initiliazing State Variables
u₀ = zeros(Float64,6)
u₀[1] = 10 #Predator; Littoral
u₀[2] = 10 #Predator; Palagic 
u₀[3] = 10 #Forager; Littoral
u₀[4] = 10 #Forager; Palagic 
u₀[5] = 10 #Juveniles; Littoral
u₀[6] = 10 #Juveniles; Palagic 

#Initiliazing Parameters
p = zeros(Float64,2)
p[1] = 0.0 #Fishing effort
p[2] = 0.0 #Dispersal Rate

#Unclear what this step does
VF(z, p) = VF!(similar(z), z, p, 0) #Insert explanation here?????????

using OrdinaryDiffEq

#Defining Time Interval
tᵢ = 0.0 #Initial Time
t_f = 1000 #Final Time
tspan = (tᵢ,t_f)

#Defining ODE problem
prob = ODEProblem(
    VF!,   #Vector Field
    u₀,    #Initial Conditions State Variables
    tspan, #Time Interval
    p)     #Initial Conditions Parameters
sol = solve(prob,Tsit5()) #Tsit5 --> Sitff ODE Solver

#Plotting Solution
solution = solve(prob, Rodas5())
    label = "", title = ["Predator-Littoral" "Predator Littoral" "Forager-Littoral" "Forager-Pelagic" "Juvenile-littoral" "Juvenile-Pelagic"])

using BifurcationKit

opt_newton = NewtonPar(tol = 1e-11, maxIter = 10)

#Continuation Options
opts_br = ContinuationPar(
    pMin = 0.0, #minimum parameter range for p
    pMax = 1.0, #maximum paramter range for p
    newtonOptions = opt_newton,
    ds = 0.0002,
    dsmax = 0.01, #maximum arclength allowed value.
    dsmin = 0.0002,
	detectBifurcation = 3, #If set to 3, a bisection algorithm is used to locate the bifurcations points (slower).
    nInversion = 8, #number of sign inversions in bisection algorithm
    maxBisectionSteps = 25, #Maximum number of bisection steps
	nev = 6, #Number of eigenvalues to be computed. 
	maxSteps = 1000, #Maximum number of continuation steps
	θ = 0.3) # θ is a parameter in the arclength constraint. 
#Initilizaing Bifurcation Problem

#Calculating Guess
z₀ = last(sol)

#Parameters for Bifurcation
par_pp2 = (
    p₁ = p[1], #fishing effort
    p₂ = p[2]) #dispersal rate

return z₀, par_pp2

#function to record information from a solution 
recordFromSolution(x, p) = (u1 = x[1], u2 = x[2],u3 = x[3],u4 = x[4], u5 = x[5], u6 = x[6]) #Insert Explanation here

prob = BifurcationProblem(
    VF, #Vector field
    z₀, #Initial guess
    par_pp2, #parameters
    (@lens _.p₁), #It specifies which parameter axis among params is used for continuation
    recordFromSolution = recordFromSolution #Insert explanation here

#Compute the bifurcation diagram associated with the problem F(x, p) = 0 recursively.
diagram = continuation(
    prob,   #bifurcation problem
    PALC(), # continuation algorithm
    plot = true,
    verbosity = 2,
    normC = norminf
return nothing
#Unclear what this step does
VF(z, p) = VF!(similar(z), z, p, 0) #Insert explanation here?????????

BifurcationKit does not use inplace formulation for now but DifferentialEquations.jl does. You use both in your code so we need both formulations.

Your code seems fine. I use deflated continuation but got a mess, I am wondering if your model specification is right. Also, why dont you use namedtuple for the parameters (or else) instead of const?

What makes you think you have 2 equilibria?

you have many branches, only the red is stable (and I cut negative u1)

Screenshot 2023-03-24 at 15.03.08

Hi @rveltz. Thank you for taking the time to respond.

I believe that I have 2 equilibria because when p_1 = 0 and p_2 = 0.8, it seems to stabilizing on a two different equilibria respectively. (Imgur: The magic of the Internet)

My understanding is that there should be a fold bifurcation.

no necessarily. It may happen that you have two disconnected branches of solutions (see picture above).

In any case, if you know the two equilibria, you can continue them and see if the branches merge