We can use my EcoEvo package, despite the fact that this doesn't seem like an ecological model.
First, install the package (only needs to be done once):
PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]
Now, load the package, set the model and add your assumptions on parameters:
<< EcoEvo`
SetModel[{
Aux[x] -> {Equation :> a/(1 + z) - Q z},
Aux[y] -> {Equation :> Q x - q y},
Aux[z] -> {Equation :> q y - c z/(k + z)}
}]
k = 1;
Q := q;
a := c (Sqrt[c/Q] - 1);
We can find equilibria with SolveEcoEq
:
eq = SolveEcoEq[]

There are two. To get an overview of where bifurcations happen in the c
--q
plane, plot where the maximum real part of the eigenvalues at each equilibrium equal zero:
λ1[c_?NumericQ, q_?NumericQ] := Max[Re[EcoEigenvalues[eq[[1]]]]];
λ2[c_?NumericQ, q_?NumericQ] := Max[Re[EcoEigenvalues[eq[[2]]]]];
ContourPlot[{λ1[c, q] == 0, λ2[c, q] == 0}, {c, -2, 2}, {q, -2, 2},
FrameLabel -> Automatic, MaxRecursion -> 3]

To investigate more closely, let's take an arbitrary slice at q = 0.2
.
q = 0.2;
Plot[λ1[c, q], {c, -1, 2}]
Plot[λ2[c, q], {c, -1, 2}]


Seems like each equilibrium has two bifurcations in this range. We can use FindRoot
to solve for them, eyeballing the graphs for initial guesses:
bif1a = FindRoot[λ1[c, q], {c, 0.19}]
bif1b = FindRoot[λ1[c, q], {c, 1.1}]
bif2a = FindRoot[λ2[c, q], {c, 0.01}]
bif2b = FindRoot[λ2[c, q], {c, 0.5}]
(* {c -> 0.2} *)
(* {c -> 1.03904} *)
(* {c -> 0.070988} *)
(* {c -> 0.466855} *)
Now we can evaluate the eigenvalues of the jacobian matrix at the bifurcation points to see if they're Hopf bifurcations (two complex conjugate eigenvalues with zero real part, non-zero imaginary part).
c = c /. bif1a;
EcoEigenvalues[eq[[1]]]
(* {-284390., -0.199832, -0.000168003} *)
c = c /. bif1b;
EcoEigenvalues[eq[[1]]]
(* {-0.42987, 0. + 0.214415 I, 0. - 0.214415 I} *)
c = c /. bif2a;
EcoEigenvalues[eq[[2]]]
(* {-0.303951, 0. + 0.144188 I, 0. - 0.144188 I} *)
c = c /. bif2b;
EcoEigenvalues[eq[[2]]]
(* {-0.358281, 0. + 0.177922 I, 0. - 0.177922 I} *)
Nope, yep, yep, yep!
We can simulate the dynamics for select values of c
to make sure.
bif2a:
c = 0.05;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.297406, 0.0143882 + 0.148584 I, 0.0143882 - 0.148584 I} *)

c = 0.1;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.316999, -0.0155364 + 0.139724 I, -0.0155364 - 0.139724 I} *)

bif2b:
c = 0.45;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.358413, -0.00121928 + 0.176904 I, -0.00121928 - 0.176904 I} *)

c = 0.48;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 10000];
PlotDynamics[sol]
(* {-0.358171, 0.000917014 + 0.178671 I, 0.000917014 - 0.178671 I} *)

bif1b:
c = 1.0;
EcoEigenvalues[eq[[1]]]
sol = EcoSim[RuleListAdd[eq[[1]], {x -> 0.1}], 1000];
PlotDynamics[sol]
(* {-0.435545, -0.00286826 + 0.213892 I, -0.00286826 - 0.213892 I} *)

c = 1.1;
EcoEigenvalues[eq[[1]]]
sol = EcoSim[RuleListAdd[eq[[1]], {x -> 0.1}], 10000];
PlotDynamics[sol]
(* {-0.422319, 0.00401044 + 0.215014 I, 0.00401044 - 0.215014 I} *)

Interesting to note that at c = 1.1
both equilibria are unstable. If we start next to eq[[2]]
we actually end up on a different limit cycle:
sol2 = EcoSim[RuleListAdd[eq[[2]], {x -> 0.1}], 10000];
PlotDynamics[sol2]

RuleListPlot[{FinalSlice[sol, 100], FinalSlice[sol2, 100]}]

So, there's a lot going on here. BTW, I don't see anything special about q = c/4
.
q
andc
)? $\endgroup$