Issue with BifurcationKit.continuation and arbitrary ODE system - MethodErrors (Convert, ambiguous, Incorrect matrix dimensions...)

Hi all! I’m trying to set up a pipeline for creating an ODE system based on a network of interactions. As a simple example, the corresponding system for the following network would be as such:
\frac{dA}{dt}=-A+\Phi(FB) \\ \frac{dB}{dt}=-B+\Phi(A)
Where A and B are processes that depend on and amplify each other, and F is some parameter (that I have been trying to bifurcate with). \Phi(x) is a sigmoid function, in my case I use:

Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))

This is my code for generating my custom ODE equations of motion:

# create equations of motion
function create_viability_model(odenetworkmodel::ODENetworkModel)
    #should be in place
    return function model!(du, u, p, t=0)
        Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))
        for i in range(1, odenetworkmodel.num_nodes)
            du[i] = -u[i]
            input = 1.0

            for j in range(1, odenetworkmodel.num_nodes)
                if odenetworkmodel.adj_matrix[i,j] == 1
                    input *= u[j]
            for k in range(1, length(odenetworkmodel.parameters))
                if odenetworkmodel.par_adj_matrix[k,i] != 0
                    input *= p[k]

            du[i] += Φ(input)

        return du

ODENetworkModel is a struct I created that holds the information that the user inputs to make such a model. If odenetworkmodel.adj_matrix[i,j] == 1, then node i depends on node j, so we multiply the input by the current value of node j (I’ve been instructed that if A depends on C, D, E, then the product of the inputs is used for the sigmoid function, i.e. \frac{dA}{dt}=-A+\Phi(CDE). The parameter adjacency matrix has a row for each parameter, so for the kth parameter, if the ith column == 1, then that parameter is being fed into node i.

struct ODENetworkModel

I have been able to use interactive_trajectory from DynamicalSystems.jl to visualise the trajectories for when there are 2 nodes (such as the example above), and even find_attractors(ds) from Attractors.jl seems to be doing fine and finding attractors that make sense for the parameters and context, BUT I can’t seem to get BifurcationKit.continuation to work! I’ve been wrestling with some dimension errors, and am currently stuck on a conversion error:

function create_bifurcation_diagram(problem)
    p = problem.parameters
    u₀ = problem.u₀

    model! = create_viability_model(problem)

    ds = ContinuousDynamicalSystem(model!, u₀, p)
    # initial_conditions = vec(collect.(Iterators.product(0:0.1:2, 0:0.1:2)))
    attractors = find_attractors(ds)
    equilibrium_guess = vec(attractors[1])

!line 76!    bfprob = BifurcationProblem(model!, equilibrium_guess, vec(p), (@lens _[problem.p_index]), record_from_solution = recordFromSolution) 

    # continuation options
	opts_br = ContinuationPar(p_min = problem.p_start, p_max = problem.p_end)

    br = BifurcationKit.continuation(bfprob, PALC(), opts_br; bothside = true)

    # scene = Plots.plot(diagram; code = (), title="Viability Bifurcation Diagram", ylim = (0.0,1.0), color=:blue, ylabel = "A")

    # display(scene)

function main()
    # don't load from file- does that fix it?
    model = ODENetworkModel(2, [0 1; 1 0], [1 0; 0 0], [0.4, 0.9], 1, 0.1, 0.9, [0.4, 0.5])


Mapping initial conditions to attractors: 100%|███████████████████████████████████████████████████████████████| Time: 0:00:00        
ERROR: MethodError: Cannot `convert` an object of type Float64 to an object of type SVector{2, Float64}

Closest candidates are:
  convert(::Type{SVector{N, T}}, ::CartesianIndex{N}) where {N, T}
   @ StaticArrays C:\Users\smahn\.julia\packages\StaticArrays\yXGNL\src\SVector.jl:8
  convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
   @ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
  convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
   @ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59

 [1] create_bifurcation_diagram(problem::ODENetworkModel)
   @ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:76
 [2] main()
   @ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:93
 [3] top-level scope
   @ REPL[5]:1

Does anyone have any tips for when you’re defining an arbitrarily large system like this? When I explicitly write out the equations of motion, I can get BifurcationKit.continuation running smoothly, but not with this method.

Many thanks for your help!!