Stability analysis of coupled ODE's

I'd like to reproduce Panel A of Figure 7 (red curve) of this paper, which represents the steady-state solutions of four coupled ODE's, and verify this claim that the diagonal branch in the Z-shape part of the curve is the unstable branch.

The governing equations for the model are:

$$\frac{dF}{dt} = a_1 - a_2 F - a_3 F Y + a_4 C, \qquad (1)$$

$$\frac{dY}{dt} = \frac{a_5}{1 + a_6 Z} - a_7 Y - a_3 F Y + a_4 C + a_8 C, \qquad (2)$$

$$\frac{dC}{dt} = a_3 F Y - a_4 C - a_9 C - a_8 C, \qquad (3)$$

$$\frac{dZ}{dt} = \frac{a_{10} a_{11} F}{1 + a_{11} F} - a_{12} Z. \qquad (4)$$

The first step would be to calculate the equilibria:

$$0 = a_1 - a_2 F - a_3 F Y + a_4 C, \qquad (5)$$

$$0 = \frac{a_5}{1 + a_6 Z} - a_7 Y - a_3 F Y + a_4 C + a_8 C, \qquad (6)$$

$$0 = a_3 F Y - a_4 C - a_9 C - a_8 C, \qquad (7)$$

$$0 = \frac{a_{10} a_{11} F}{1 + a_{11} F} - a_{12} Z, \qquad (8)$$

which in Mathematica, we have:

Solve[{a1 - a2*F - a3*F*Y + a4*C,
a5/(1 + a6*Z) - a7*Y - a*3*F*Y + a4*C + a8*C,
a3*F*Y - a4*C - a9*C - a8*C,
a10*a11*F/(1 + a11*F) - a12*Z} == {0, 0, 0, 0}, {F, Y, C, Z}]


Mathematica doesn't return any solution for the above (it seems that it is running forever), besides this, perhaps the solutions of the above algebraic equations are not helpful at this stage, since in the figure cited above, we have a two-parameter bifurcation diagram; thus, we need further simplifications. Before going further, let's also calculate the Jacobian matrix at this stage, that is,

$$J=\begin{pmatrix} d\dot F\over dF & d\dot F\over dY & d\dot F\over dC & d\dot F\over dZ \\ d\dot Y\over dF & d\dot Y\over dY & d\dot Y\over dC & d\dot Y\over dZ \\ d\dot C\over dF & d\dot C\over dY & d\dot C\over dC & d\dot C\over dZ \\ d\dot Z\over dF & d\dot Z\over dY & d\dot Z\over dC & d\dot Z\over dZ \\ \end{pmatrix},$$

D[{a1 - a2*F - a3*F*Y + a4*C,
a5/(1 + a6*Z) - a7*Y - a*3*F*Y + a4*C + a8*C,
a3*F*Y - a4*C - a9*C - a8*C,
a10*a11*F/(1 + a11*F) - a12*Z}, {{F, Y, C, Z}}]


To simplify the above algebraic equations further, by exploiting $$(7)$$, we obtain for $$C$$: $$C = a_{13} F Y$$, where $$a_{13} = a_3 / (a_4 + a_8 + a_9)$$. Now, by adding $$(7)$$ to $$(5)$$ and $$(6)$$, we obtain:

$$a_1 = a_2 F + (a_8 + a_9) a_{13} F Y, \qquad (9)$$

$$\frac{a_5}{1 + a_6 Z} = a_7 Y + a_9 a_{13} F Y. \qquad (10)$$

By using the following transformations and definition:

$$a_{11} F \to F, \quad a_{13} Y \to Y, \quad a_6 Z \to Z, \quad a_1 / a_2 \to a_1, \quad a_{10} / a_{12} \to a_{10}, \quad a_5 / a_7 \to a_5, \quad (a_8 + a_9) / a_2 \to a_9, \quad a_{14} \equiv a_9 / a_7,$$

we can rewrite $$(8)$$, $$(9)$$, and $$(10)$$ as:

$$a_1 = F + a_9 F Y, \qquad (11)$$

$$\frac{a_5}{1 + Z} = Y + a_{14} F Y, \qquad (12)$$

$$\frac{a_{10} F}{1 + F} = Z. \qquad (13)$$

Now, by manipulating $$(11)$$, $$(12)$$, and $$(13)$$, we obtain:

$$a_9 a_{10}^2 a_5 C_2 (1 - C_2) - (1 + Z) [a_1 a_{10} - (1 + a_1) Z] [a_{10} - (1 - a_{14}) Z] = 0,$$

where $$C_2 \equiv F / (1 + F)$$; the above can be further simplified to:

$$a_9 a_{10}^2 a_5 C_2 (1 - C_2) - (1 + a_{10} C_2) [a_1 a_{10} - (1 + a_1) a_{10} C_2] [a_{10} - (1 - a_{14}) a_{10} C_2] = 0. (14)$$

Now, for the specific values of $$a_1 = 3$$, $$a_9 = 150$$, $$a_{10} = 6$$, and $$a_{14} = 500$$, we can obtain $$C_2$$ vs. $$a_5$$ by exploiting $$(14)$$ as:

eqn = 150*36*a5*C2*(1 - C2) - (1 + 6*C2) (3*6 - (1 + 3) 6*C2) (6 - (1 - 500) 6*C2) == 0;

Column[sol = Solve[eqn, C2, Reals]];

Plot[Evaluate[C2 /. sol], {a5, 0, 28},
PlotRange -> {0, 0.8},
PlotPoints -> 75,
MaxRecursion -> 5,
WorkingPrecision -> 15,
PlotLegends -> Placed[Automatic, {0.25, 0.45}],
AxesOrigin -> {0, 0},
AxesLabel -> {Subscript[a, 5], Subscript[C, 2]},
LabelStyle -> {16, GrayLevel[0]}, PlotStyle -> Thickness[0.0098]]


which results in:

thus, we have correctly reproduced the cited figure at the beginning of the post. Now, the question which remains is that, how to show the orange curve is the unstable branch? To this end, we need to calculate the above-obtained Jacobian at different equilibria and then calculate its eigenvalues or to study the Routh-Hurwitz stability criteria. Here is the place that I'm stuck; I cannot simplify the Jacobian in order to get rid of extra parameters and recast it in terms of the parameters and equations that we have obtained at the end. Any help, idea, or hint, is really appreciated!

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 and set the model.

<< EcoEvo

SetModel[{
Aux[F] -> {Equation :> βF - F - γC F Y},
Aux[Y] -> {Equation :> βY/(1 + Z) - Y - γCY F Y},
Aux[Z] -> {Equation :> βZ F/(1 + F) - Z}
}]


Note that I'm going straight to the nondimensionalized system insinuated by eqns (12)-(14) of the supplemental material:

$${dF \over dt}=\beta_F-F-\gamma_CFY \\ {dY \over dt}={\beta_Y \over 1+Z}-Y-\gamma_{C_Y}FY \\ {dZ \over dt}={\beta_Z F \over 1+F}-Z$$

Before we set parameter values, we can look at the generic Jacobian matrix:

EcoJacobian[]//MatrixForm


Now set parameter values, solve for equilibria, and plot them:

βZ = 6; γC = 150; γCY = 500; βF = 1;

eq = SolveEcoEq[SolveOpts -> {Reals}];

Plot[Evaluate[F/(1 + F) /. eq], {βY, 0, 7},
PlotRange -> {0, All}, AxesLabel -> {"βY", "F/(1+F)"}]


Because the algebra in the general case is nasty, I'm just going to pick a parameter value in the bistable region to check stability.

βY = 4.5;
eq = SolveEcoEq[]
EcoEigenvalues[eq]

(* {{F -> 0.00648905, Y -> 1.02071, Z -> 0.0386833},
{F -> 0.0759319, Y -> 0.0811313, Z -> 0.423439},
{F -> 0.579865, Y -> 0.00483027, Z -> 2.20221}} *)

(* {{-157.349, -1.39267, -0.608318},
{-51.0834, -2.64431, 0.592042},
{-291.656, -1.56241, -0.438678}} *)


From the eigenvalues, we can see that eq[[2]] is unstable and eq[[1]] and eq[[3]] are stable. If you want to see the numerical Jacobian matrices whose eigenvalues these are:

EcoJacobian[eq[[1]]] // MatrixForm
EcoJacobian[eq[[2]]] // MatrixForm
EcoJacobian[eq[[3]]] // MatrixForm


Finally, based on our suspicion that eq[[2]] is unstable throughout its range (remember, we only tested one value of $$\beta_Y$$, but it's hard to imagine how it could be otherwise), we can style the bifurcation diagram:

Clear[βY];
Plot[Evaluate[F/(1 + F) /. eq], {βY, 0, 7},
PlotRange -> {0, All}, AxesLabel -> {"βY", "F/(1+F)"},
PlotStyle -> {Black, {Black, Dashed}, Black}]
`