From T=TSE(E)BL(1T)BPP\dot T = \alpha_T S_E(E) B_L (1-T) - B_P PT=TSE(E)BL(1T)BPP:
TT=TSE(E)BL,ET=TBL(1T)SE(E),PT=BP.\begin{aligned} \partial_T \dot T &= -\alpha_T S_E(E^*) B_L,\\[4pt] \partial_E \dot T &= \alpha_T B_L (1-T^*) S_E'(E^*),\\[4pt] \partial_P \dot T &= -B_P. \end{aligned}TTETPT=TSE(E)BL,=TBL(1T)SE(E),=BP.
From E=EEE+PPTT\dot E = \mu_E - \lambda_E E + \phi_P P - \phi_T TE=EEE+PPTT:
TE=T,EE=E,PE=P.\begin{aligned} \partial_T \dot E &= -\phi_T,\\[4pt] \partial_E \dot E &= -\lambda_E,\\[4pt] \partial_P \dot E &= \phi_P. \end{aligned}TEEEPE=T,=E,=P.
From P=PP(1P/Pmax)+ESE(E)TTKK(E)P\dot P = \rho_P P(1-P/P_{\max}) + \eta_E S_E(E) - \delta_T T - \kappa_K K(E) PP=PP(1P/Pmax)+ESE(E)TTKK(E)P:
Compute derivative carefully: logistic term derivative: P[PP(1P/Pmax)]=P(12P/Pmax)\partial_P [\rho_P P(1-P/P_{\max})] = \rho_P(1 - 2P/P_{\max})P[PP(1P/Pmax)]=P(12P/Pmax).
Also P[KK(E)P]=KK(E)\partial_P[-\kappa_K K(E) P] = -\kappa_K K(E)P[KK(E)P]=KK(E) (here we treated K\kappa_KK constant; if K\kappa_KK depends on P via K(E) only E dependent so ok). And E[ESE(E)KK(E)P]=ESE(E)KPK(E)\partial_E [\eta_E S_E(E) - \kappa_K K(E) P] = \eta_E S_E'(E) - \kappa_K P \,K'(E)E[ESE(E)KK(E)P]=ESE(E)KPK(E).
Hence
TP=T,EP=ESE(E)KPK(E),PP=P(12PPmax)KK(E).\begin{aligned} \partial_T \dot P &= -\delta_T,\\[4pt] \partial_E \dot P &= \eta_E S_E'(E^*) \;-\; \kappa_K P^* K'(E^*),\\[4pt] \partial_P \dot P &= \rho_P\Big(1 - \dfrac{2P^*}{P_{\max}}\Big) \;-\; \kappa_K K(E^*). \end{aligned}TPEPPP=T,=ESE(E)KPK(E),=P(1Pmax2P)KK(E).
Collecting all:
J(X)=(TSE(E)BLTBL(1T)SE(E)BPTEPTESE(E)KPK(E)P(12P/Pmax)KK(E)).(11)\boxed{% J(X^*) = \begin{pmatrix} -\alpha_T S_E(E^*) B_L & \alpha_T B_L (1-T^*) S_E'(E^*) & -B_P\\[6pt] -\phi_T & -\lambda_E & \phi_P\\[6pt] -\delta_T & \eta_E S_E'(E^*) - \kappa_K P^* K'(E^*) & \rho_P(1 - 2P^*/P_{\max}) - \kappa_K K(E^*) \end{pmatrix}. } \tag{11}J(X)=TSE(E)BLTTTBL(1T)SE(E)EESE(E)KPK(E)BPPP(12P/Pmax)KK(E).(11)
Here SE(E)=dSEdE=E1/(E+)2S_E'(E) = \dfrac{dS_E}{dE} = \kappa\,\theta^\kappa\,E^{\kappa-1}/(E^\kappa+\theta^\kappa)^2SE(E)=dEdSE=E1/(E+)2 (or the equivalent derivative of the Hill function used). K(E)K'(E)K(E) is derivative of the cooptation function w.r.t. EEE (zero for our original KKK that depends on leadership parameters only --- if KKK depends only on \Theta and not on EEE, then K(E)=0K'(E)=0K(E)=0, simplifying the expression).
6. Bifurcation conditions and nondegeneracy
Fold (saddle-node): algebraically encoded by (9a--9b) (double root in PPP for some EE^*E). Equivalently, at (Xc,c)(X_c,\mu_c)(Xc,c) the Jacobian has a simple zero eigenvalue: detJ(Xc)=0\det J(X_c)=0detJ(Xc)=0 and rankJJJ=2. Nondegeneracy requires the quadratic nonlinearity coefficient a=12wD2F[v,v]0a = \tfrac12 w^\top D^2F[v,v] \neq 0a=21wD2F[v,v]=0 (see Section 2.D) and b=wF0b = w^\top \partial_{\mu}F \neq 0b=wF=0. For our model F=(0,1,0)\partial_{\mu}F = (0,1,0)^\topF=(0,1,0), so b=w2b= w_2b=w2 (the second component of the left nullvector).
Hopf (if it appears): occurs when the Jacobian has a simple complex conjugate pair crossing the imaginary axis. Routh--Hurwitz criteria on the cubic characteristic polynomial (or direct eigenvalue computation) determine the Hopf locus. The first Lyapunov coefficient 1\ell_11 computed from multilinear forms determines criticality.
7. Recipe to compute crit()\mu_{\mathrm{crit}}(\Theta)crit() analytically / numerically
