BifurcationKit.jl uses numerical continuation to produce the bifurcation diagram and numerical continuation can only find connected branches of solutions. In this case, your run of BifurcationKit finds an initial branch of unstable steady-states which doesn’t undergo any bifurcations. It might find other branches if you vary the initial conditions. That said, the code you’ve used will only find steady-state solutions though (equilibria) - you need a slightly different code to find/track periodic orbits. (BifurcationKit cannot track chaotic attractors directly.)
The code I used to check that the branch of solutions computed is unstable is
using NonlinearSolve
using ForwardDiff
prob = NonlinearProblem(HR, zeros(6), par) # root finding
sol = solve(prob, NewtonRaphson()) # find the steady-state
J = ForwardDiff.jacobian(u -> HR(u, par), sol.u) # find it's linearised stability
isstable = all(real.(eigen(J).values) .< 0) # check all eigenvalues have negative real part
Given the plots you’ve made of the stable orbits you see from simulation, you might be best to do a brute-force bifurcation diagram, at least to start with. That is, use Monte-Carlo simulation over the parameter range you are interested in and plot some summary statistics against the bifurcation parameter.