\section{The mechanism, its equilibria, and a determinant formula}\label{sec:model}
Let \(V=V_\kappa\) be the mass-action vector field of \eqref{eq:network} on \(\R^{3n+3}\), with
coordinates
\[
 (S_0,\dots,S_n,\;E,\;F,\;C_1,\dots,C_n,\;Y_1,\dots,Y_n),
\]
and let \(W\colon\R^{3n+3}\to\R^3\) be the linear map \eqref{eq:totals}. The totals are
conserved, so \(V\) takes values in \(\ker W\). For positive \(T=(\ET,\FT,\ST)\) the
compatibility class is \(\mathcal P_T=\{X\ge0:\ WX=T\}\).

\subsection{Pools, the stoichiometric subspace, and the boundary}
For \(1\le i\le n\) let
\begin{equation}\label{eq:pools}
 Z_i=\sum_{k=i}^n S_k+\sum_{k=i+1}^n C_k+\sum_{k=i}^n Y_k
\end{equation}
be the amount of substrate carrying at least \(i\) phosphate groups, where \(C_k\) is counted
with \(k-1\) groups and \(Y_k\) with \(k\). Binding and unbinding do not change \(Z_i\), the
kinase catalytic step at site \(i\) increases it, and the phosphatase catalytic step at site
\(i\) decreases it; hence
\begin{equation}\label{eq:poolcurrent}
 \dot Z_i=c_iC_i-\gamma_iY_i .
\end{equation}
The functions \((C_1,\dots,C_n,Y_1,\dots,Y_n,Z_1,\dots,Z_n)\) are linear coordinates on
\(\ker W\): from them and the totals one recovers \(S_n,\dots,S_1\) successively from
\eqref{eq:pools}, and then \(E,F,S_0\) from \eqref{eq:totals}. In these coordinates a binding
reaction changes exactly one complex coordinate and no \(Z_k\), and the catalytic reactions at
site \(i\) change \(Z_i\) and no other \(Z_k\). Therefore the \(6n\) reaction vectors span
\(\ker W\), and \(\ker W\) is the stoichiometric subspace; its dimension is \(3n\).

\begin{lemma}\label{lem:boundary}
Let \(T\in\R^3_{>0}\). The class \(\mathcal P_T\) is a compact convex polytope with nonempty
relative interior, the vector field is inward or tangent on each of its faces, and every
equilibrium in \(\mathcal P_T\) is strictly positive.
\end{lemma}
\begin{proof}
Every concentration is bounded by one of the totals, so \(\mathcal P_T\) is compact; it is
convex and closed by definition. An interior point is obtained by giving all complexes a small
positive concentration and distributing the remaining substrate positively among the \(S_i\).
If a species has concentration zero, every reaction that consumes it has rate zero, so its
derivative is nonnegative; this is the face condition. At an equilibrium,
\(\dot C_i=\dot Y_i=0\) gives \(C_i=p_iS_{i-1}E\) and \(Y_i=q_iS_iF\) with
\(p_i=a_i/(b_i+c_i)\) and \(q_i=\alpha_i/(\beta_i+\gamma_i)\). If \(E=0\), then all \(C_i\)
vanish and \(\ET=0\); so \(E>0\), and likewise \(F>0\). By \eqref{eq:poolcurrent},
\(c_ip_iES_{i-1}=\gamma_iq_iFS_i\), so the \(S_i\) are all positive or all zero, and the latter
would give \(\ST=0\). Then the complexes are positive as well.
\end{proof}

\subsection{Parametrization of the positive equilibria}
For given rate constants put
\begin{equation}\label{eq:pq}
 p_i=\frac{a_i}{b_i+c_i},\qquad q_i=\frac{\alpha_i}{\beta_i+\gamma_i},\qquad
 t_i=\prod_{k=1}^{i}\frac{c_kp_k}{\gamma_kq_k}\quad(t_0=1),
\end{equation}
and define three polynomials with positive coefficients,
\begin{equation}\label{eq:coeff}
\begin{gathered}
 A(u)=\sum_{i=0}^n t_i u^i,\qquad
 B(u)=\sum_{i=0}^{n-1} B_i u^i,\qquad
 D(u)=\sum_{i=0}^{n-1} D_i u^i,\\
 B_{i-1}=p_it_{i-1},\qquad D_{i-1}=q_it_i\qquad(1\le i\le n).
\end{gathered}
\end{equation}
The letter \(D\) always denotes this polynomial; the phosphatase complexes are called \(Y_i\).

\begin{proposition}\label{prop:param}
\textup{(i)} The positive equilibria of \eqref{eq:network} are exactly the points
\begin{equation}\label{eq:states}
\begin{gathered}
 S_i=t_i u^i s \quad(0\le i\le n),\qquad E=uf,\qquad F=f,\\
 C_i=B_{i-1}u^i sf,\qquad Y_i=D_{i-1}u^i sf \quad(1\le i\le n),
\end{gathered}
\end{equation}
with \((u,s,f)\in\R^3_{>0}\). Their totals are
\begin{equation}\label{eq:totalcoords}
 \ET=uf\bigl(1+sB(u)\bigr),\qquad \FT=f\bigl(1+suD(u)\bigr),\qquad
 \ST=sA(u)+usf\bigl(B(u)+D(u)\bigr).
\end{equation}
\textup{(ii)} Conversely, let \(t_0=1\) and let \(t_1,\dots,t_n\), \(B_0,\dots,B_{n-1}\),
\(D_0,\dots,D_{n-1}\) be arbitrary positive numbers. For any positive numbers
\(b_i,\beta_i,\gamma_i\) \((1\le i\le n)\), the rate constants
\begin{equation}\label{eq:ratesgeneral}
 c_i=\gamma_i\frac{D_{i-1}}{B_{i-1}},\qquad
 a_i=(b_i+c_i)\frac{B_{i-1}}{t_{i-1}},\qquad
 \alpha_i=(\beta_i+\gamma_i)\frac{D_{i-1}}{t_i}
\end{equation}
are positive and have the given numbers as their coefficients \eqref{eq:coeff}.
\end{proposition}
\begin{proof}
At an equilibrium, \(C_i=p_iS_{i-1}E\), \(Y_i=q_iS_iF\) and, by \eqref{eq:poolcurrent},
\(c_iC_i=\gamma_iY_i\), which gives \(S_i=(t_i/t_{i-1})(E/F)S_{i-1}\). With \(u=E/F\),
\(s=S_0\), \(f=F\) this is \eqref{eq:states}. Conversely, at a point of the form
\eqref{eq:states} the complex components of \(V\) vanish and \(c_iC_i=\gamma_iY_i\) for every
\(i\). Then \(a_iS_{i-1}E-b_iC_i=c_iC_i\) and \(\alpha_iS_iF-\beta_iY_i=\gamma_iY_i\), so
\(\dot S_i=-c_{i+1}C_{i+1}+c_iC_i-\gamma_iY_i+\gamma_{i+1}Y_{i+1}=0\) (terms with index \(0\)
or \(n+1\) are absent), and \(\dot E=\dot F=0\) because the complexes are stationary. Summation
gives \eqref{eq:totalcoords}. Part (ii) is a substitution: \eqref{eq:ratesgeneral} gives
\(p_i=B_{i-1}/t_{i-1}\), \(q_i=D_{i-1}/t_i\) and \(c_ip_i/(\gamma_iq_i)=t_i/t_{i-1}\).
\end{proof}

Thus the positive equilibria, and their number in each class, depend on the \(6n\) rate
constants only through the \(3n\) coefficients of \(A,B,D\). The remaining \(3n\) degrees of
freedom \(b_i,\beta_i,\gamma_i\) affect the dynamics, including stability, but not the
equilibria. We call the coefficients the \emph{geometry} and \((b_i,\beta_i,\gamma_i)\) the
\emph{kinetics}. This separation is the organizing principle of the paper.

\subsection{The scalar equation}
Fix positive totals and put \(r=\ET/\FT\),
\begin{equation}\label{eq:LM}
 L(u)=B(u)-rD(u),\qquad M(u)=B(u)-uD(u).
\end{equation}
Dividing the first two equations of \eqref{eq:totalcoords} gives \(r-u=suL(u)\). Hence a
positive equilibrium with \(u\ne r\) satisfies \((r-u)L(u)>0\), and then
\begin{equation}\label{eq:chart}
 s(u)=\frac{r-u}{uL(u)},\qquad
 f(u)=\frac{\FT}{1+s(u)uD(u)}=\frac{\FT L(u)}{M(u)},
\end{equation}
where we used \(L+(r-u)D=M\); in particular \(M(u)=L(u)\bigl(1+suD(u)\bigr)\) is nonzero and
has the sign of \(L(u)\). Conversely, for every \(u>0\) with \((r-u)L(u)>0\) the formulas
\eqref{eq:chart} define positive \(s,f\) such that \eqref{eq:states} has enzyme totals
\(\ET,\FT\). On the open set \(\mathcal U=\{u>0:(r-u)L(u)>0\}\) define the
\emph{substrate-total function}
\begin{equation}\label{eq:H}
 \Phi(u)=s(u)A(u)+u\,s(u)f(u)\bigl(B(u)+D(u)\bigr).
\end{equation}
The positive equilibria with totals \((\ET,\FT,\ST)\) and \(u\ne r\) correspond bijectively to
the solutions \(u\in\mathcal U\) of \(\Phi(u)=\ST\). This reduction is the one used in
\cite{WangSontag2008}, from which we take the following bound.

\begin{theorem}[Wang and Sontag {\cite[Theorem~3]{WangSontag2008}}]\label{thm:WS}
For every \(n\ge1\), all positive rate constants and all positive totals, the system
\eqref{eq:network} has at most \(2n-1\) positive equilibria with those totals.
\end{theorem}

\subsection{The restricted Jacobian}
Put \(y=(S_1,\dots,S_n,C_1,\dots,C_n,Y_1,\dots,Y_n)\) and \(z=(S_0,E,F)\), and order the totals
as \((\ST,\ET,\FT)\). Then \(W=(W_y\;\;I_3)\) in the coordinates \((y,z)\), so \(\ker W\) is
the graph of \(z=-W_yy\), and \(y\) is a linear coordinate system on it. Let \(\JS\) denote the
matrix of \(DV|_{\ker W}\) in the coordinates \(y\), that is, \(\JS=V_{y,y}-V_{y,z}W_y\), where
\(V_y\) denotes the \(y\)-components of \(V\). An equilibrium is \emph{nondegenerate} if
\(\det\JS\ne0\), \emph{hyperbolic} if \(\JS\) has no eigenvalue on the imaginary axis, and a
\emph{sink} if \(\JS\) is Hurwitz.

\begin{proposition}[Determinant formula]\label{prop:det}
For arbitrary positive rate constants, at any positive equilibrium with
\(u=E/F\ne r=\ET/\FT\),
\begin{equation}\label{eq:det}
 \det\JS=(-1)^{n+1}\Bigl(\prod_{i=1}^n(b_i+c_i)\,\alpha_i\,\gamma_i\Bigr)\,F^{\,n}\,u\,M(u)\,\Phi'(u),
\end{equation}
where \(\Phi\) is formed with the enzyme totals of that equilibrium. In particular the
equilibrium is nondegenerate if and only if \(\Phi'(u)\ne0\).
\end{proposition}
\begin{proof}
By the Schur complement with respect to the block \(I_3\),
\[
 \det\JS=\det\begin{pmatrix}V_{y,y}&V_{y,z}\\ W_y&I_3\end{pmatrix}.
\]
Let \(g=(g_Z,g_C,g_Y)\) consist of the \(3n\) functions \(g_{Z,i}=c_iC_i-\gamma_iY_i\),
\(g_{C,i}=\dot C_i\), \(g_{Y,i}=\dot Y_i\). By \eqref{eq:pools}--\eqref{eq:poolcurrent},
\(g=UV_y\) for a constant unitriangular matrix \(U\), since \(g_{Z,i}\) is \(\dot S_i\) plus
components of \(V_y\) that come later in the ordering of \(y\). Hence the first \(3n\) rows can
be replaced by \((g_y\;\;g_z)\) without changing the determinant. In the block form belonging
to \(y=(S;C,Y)\),
\[
 g_y=\begin{pmatrix}0&\ast\\ \ast&\Delta\end{pmatrix},\qquad
 \Delta=\diag\bigl(-(b_i+c_i);\,-(\beta_i+\gamma_i)\bigr),
\]
and the Schur complement of \(\Delta\) is the Jacobian with respect to \(S_1,\dots,S_n\) of
the reduced currents \(c_ip_iS_{i-1}E-\gamma_iq_iS_iF\). It is lower bidiagonal with diagonal
entries \(-\gamma_iq_iF\). Therefore
\[
 \det g_y=\prod_i(b_i+c_i)(\beta_i+\gamma_i)\cdot(-1)^n\prod_i\gamma_iq_iF
 =(-1)^nF^n\prod_i(b_i+c_i)\alpha_i\gamma_i\ne0 .
\]
In particular the zero set of \(g\) is locally the graph of a function \(y=\varphi(z)\), in
agreement with Proposition~\ref{prop:param}, and
\[
 \det\begin{pmatrix}g_y&g_z\\ W_y&I_3\end{pmatrix}
 =\det g_y\cdot\det\bigl(I_3-W_yg_y^{-1}g_z\bigr)
 =\det g_y\cdot\det\frac{\partial(\ST,\ET,\FT)}{\partial(S_0,E,F)},
\]
where the last matrix is the derivative of the totals along the equilibrium manifold
\(z\mapsto(\varphi(z),z)\), because \(\varphi_z=-g_y^{-1}g_z\). We compute it in the
coordinates \((u,s,f)=(E/F,S_0,F)\), for which
\(\det\partial(u,s,f)/\partial(S_0,E,F)=-1/F\). From \eqref{eq:totalcoords},
\begin{equation}\label{eq:pivot}
 \det\frac{\partial(\ET,\FT)}{\partial(s,f)}
 =\det\begin{pmatrix}ufB&u(1+sB)\\ fuD&1+suD\end{pmatrix}=uf\,(B-uD)=uf\,M(u)\ne0,
\end{equation}
so near the given equilibrium the enzyme totals determine \((s,f)\) as functions of \(u\),
namely \eqref{eq:chart}, and the Schur complement of the block \eqref{eq:pivot} in
\(\partial(\ST,\ET,\FT)/\partial(u,s,f)\) is the derivative of \(\ST\) along that curve, which
is \(\Phi'(u)\). Hence \(\det\partial(\ST,\ET,\FT)/\partial(u,s,f)=ufM(u)\Phi'(u)\), and
collecting the factors gives \eqref{eq:det}.
\end{proof}

\begin{corollary}\label{cor:parity}
At a positive equilibrium with \(u<r\) and \(\Phi'(u)\ne0\), the number of positive real
eigenvalues of \(\JS\), counted with multiplicity, is even if \(\Phi'(u)<0\) and odd if
\(\Phi'(u)>0\). In particular an equilibrium with \(u<r\) and \(\Phi'(u)>0\) is unstable.
\end{corollary}
\begin{proof}
Nonreal eigenvalues come in conjugate pairs, so \(\sign\det\JS=(-1)^{e_-}\), where \(e_-\) is the
number of negative real eigenvalues, and \(3n-e_-\) has the parity of the number \(e_+\) of
positive real eigenvalues. Thus \(\sign\det\JS=(-1)^{n}(-1)^{e_+}\). On the other hand
\(M(u)\) has the sign of \(L(u)\), hence of \(r-u\), and \eqref{eq:det} gives
\(\sign\det\JS=(-1)^{n+1}\sign\Phi'(u)\).
\end{proof}

\section{A universal ceiling for attracting equilibria}\label{sec:ceiling}
This section proves Theorem~\ref{thm:main}(a). It is independent of the construction. We use
the Brouwer degree \(\deg(g,\Omega,0)\) of a continuous map on a bounded open set without zeros
on the boundary, and the index \(\ind(g,X)\) of an isolated zero \cite{Deimling,KrasnoselskiiZabreiko}.
At a zero with invertible derivative the index is the sign of the determinant of the
derivative. Fix positive totals \(T\), use \(y\) as a coordinate system on the affine hull of
\(\mathcal P=\mathcal P_T\), and write \(g_\kappa\) for the vector field in these coordinates;
its derivative at an equilibrium is \(\JS\). We always take degrees and indices of \(-g_\kappa\),
so that a sink has index \(\sign\det(-\JS)=+1\).

\begin{lemma}\label{lem:degree}
\textup{(i)} \(\deg(-g_\kappa,\operatorname{int}\mathcal P,0)=1\) for every
\(\kappa\in\R^{6n}_{>0}\).
\textup{(ii)} The set of \(\kappa\in\R^{6n}_{>0}\) for which all equilibria in \(\mathcal P\) are
nondegenerate is dense.
\textup{(iii)} For such \(\kappa\) the equilibria are finite in number and at most \(n\) of them
have index \(+1\).
\end{lemma}
\begin{proof}
(i) Let \(y_0\in\operatorname{int}\mathcal P\) and consider
\(g^{t}(y)=(1-t)g_\kappa(y)+t(y_0-y)\), \(0\le t\le1\). On a face of \(\mathcal P\) where some
species vanishes, the corresponding component of \(g_\kappa\) is nonnegative by
Lemma~\ref{lem:boundary} and that of \(y_0-y\) is positive, so \(g^t\) has no zero there for
\(t>0\); for \(t=0\) there is none by Lemma~\ref{lem:boundary}. (The species \(S_0,E,F\) are
affine functions of \(y\), and the same argument applies to their time derivatives.) By
homotopy invariance, \(\deg(-g_\kappa,\operatorname{int}\mathcal P,0)=\deg(y-y_0,\operatorname{int}\mathcal P,0)=1\).

(ii) The vector field is linear in the rate constants:
\(V_\kappa(X)=\Gamma\,\diag\bigl(X^{\eta_1},\dots,X^{\eta_{6n}}\bigr)\kappa\), where \(\Gamma\) is
the stoichiometric matrix and \(X^{\eta_k}\) are the reactant monomials. At a positive \(X\) the
monomials are positive and \(\Gamma\) has rank \(3n=\dim\ker W\) by Section~\ref{sec:model}. Hence
\((y,\kappa)\mapsto g_\kappa(y)\) is a submersion on
\(\operatorname{int}\mathcal P\times\R^{6n}_{>0}\), and by the parametric transversality theorem
\cite[Ch.~2, \S3]{GuilleminPollack} zero is a regular value of \(g_\kappa\) on
\(\operatorname{int}\mathcal P\) for almost every \(\kappa\). By Lemma~\ref{lem:boundary} there
are no other equilibria.

(iii) Nondegenerate zeros are isolated, and an accumulation point of zeros would be a
non-isolated zero in \(\operatorname{int}\mathcal P\) or a zero on the boundary; so the zeros
are finite in number. Let \(n_\pm\) be the numbers of zeros of index \(\pm1\). Additivity of the
degree and (i) give \(n_+-n_-=1\), Theorem~\ref{thm:WS} gives \(n_++n_-\le2n-1\), and therefore
\(n_+\le n\).
\end{proof}

The degree identity in (i), for dissipative networks without boundary equilibria, is the basis
of the multistationarity criterion of \cite{CFMW2017}; see also \cite{Hofbauer1990}. Positive
index does not imply stability, so (iii) alone does not bound the number of attractors in a
degenerate system. The next lemma closes this gap.

\begin{lemma}\label{lem:index}
Let \(X_*\in\operatorname{int}\mathcal P\) be an equilibrium that is locally asymptotically
stable relative to \(\mathcal P\). Then \(X_*\) is an isolated zero of \(g_\kappa\) and
\(\ind(-g_\kappa,X_*)=+1\).
\end{lemma}
\begin{proof}
This is classical \cite{KrasnoselskiiZabreiko}; we include the short argument. Let \(\phi_t\)
be the flow and \(\mathcal B\) a closed ball around \(X_*\) inside its basin of attraction and
inside \(\operatorname{int}\mathcal P\). No point of \(\mathcal B\setminus\{X_*\}\) is an
equilibrium or lies on a periodic orbit, because its orbit converges to \(X_*\). Hence
\(I-\phi_t\) has no zero on \(\partial\mathcal B\) for any \(t>0\), and
\(\deg(I-\phi_t,\operatorname{int}\mathcal B,0)\) does not depend on \(t>0\). Stability and
attraction imply that attraction is uniform on the compact set \(\mathcal B\), so for large
\(t\) we have \(\phi_t(\mathcal B)\subset\operatorname{int}\mathcal B'\) for a smaller
concentric ball \(\mathcal B'\); the homotopy \(I-\theta(\phi_t-X_*)-X_*\),
\(0\le\theta\le1\), then has no zero on \(\partial\mathcal B\) and shows that the degree is
\(1\). On the other hand \((I-\phi_t)/t\to-g_\kappa\) uniformly on \(\partial\mathcal B\) as
\(t\downarrow0\), and \(g_\kappa\) does not vanish there, so
\(\deg(-g_\kappa,\operatorname{int}\mathcal B,0)=1\) as well.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:main}\textup{(a)}]
Suppose that for some \(\kappa\) and \(T\) there are \(n+1\) locally asymptotically stable
equilibria. They are positive by Lemma~\ref{lem:boundary}. Choose pairwise disjoint closed balls
\(\mathcal B_0,\dots,\mathcal B_n\) around them as in Lemma~\ref{lem:index}. Since
\(g_\kappa\ne0\) on the compact set \(\bigcup_k\partial\mathcal B_k\) and \(g_{\kappa'}\) depends
continuously on \(\kappa'\), we have
\(\deg(-g_{\kappa'},\operatorname{int}\mathcal B_k,0)=1\) for all \(\kappa'\) close to \(\kappa\)
and all \(k\). By Lemma~\ref{lem:degree}(ii) we may choose such a \(\kappa'\) for which all
equilibria are nondegenerate. Then each \(\mathcal B_k\) contains at least one equilibrium of
index \(+1\), so there are at least \(n+1\) of them, which contradicts
Lemma~\ref{lem:degree}(iii).
\end{proof}

The argument does not assume that a degenerate attractor persists under perturbation; it only
uses the persistence of its index. The model-specific ingredients are the exclusion of boundary
equilibria, the rank of the stoichiometric matrix, and Theorem~\ref{thm:WS}. It follows that
\(\capH\le\capC\le n\).

\section{Construction of \texorpdfstring{\(2n-1\)}{2n-1} equilibria}\label{sec:construction}

\subsection{A locus on which the scalar equation factors}
\begin{lemma}\label{lem:factor}
Suppose \(A(u)=(1+u)D(u)\) and \(\ST=\ET+\FT\), and normalize \(\FT=2\), so that \(\ET=2r\) and
\(\ST=2(r+1)\). For \(u\in\mathcal U\) put \(\ell=B(u)/D(u)-r\), \(\eta=r-u\) and
\(x=\sqrt{1+8u}\). Then
\begin{equation}\label{eq:residual}
 \Phi(u)-\ST=\frac{(1+u)\bigl(\eta(\ell+\eta)-2u\ell^2\bigr)}{u\,\ell\,(\ell+\eta)}
 =-\frac{2(1+u)}{\ell(\ell+\eta)}\Bigl(\ell-\frac{2\eta}{x-1}\Bigr)\Bigl(\ell+\frac{2\eta}{x+1}\Bigr).
\end{equation}
\end{lemma}
\begin{proof}
Since \(L=D\ell\) and \(M=D(\ell+\eta)\), \eqref{eq:chart} gives \(s=\eta/(uD\ell)\) and
\(f=2\ell/(\ell+\eta)\). Therefore
\[
 \Phi=\frac{(1+u)\eta}{u\ell}+\frac{2\eta(\ell+r+1)}{\ell+\eta},\qquad
 \Phi-2(r+1)=\frac{(1+u)\eta}{u\ell}-\frac{2\ell(r+1-\eta)}{\ell+\eta}
 =(1+u)\Bigl(\frac{\eta}{u\ell}-\frac{2\ell}{\ell+\eta}\Bigr),
\]
because \(r+1-\eta=1+u\). This is the first expression. For the second, expand the product and
use \(x^2-1=8u\):
\(-2u\bigl(\ell^2-\tfrac{4\eta}{x^2-1}\ell-\tfrac{4\eta^2}{x^2-1}\bigr)=-2u\ell^2+\eta\ell+\eta^2\).
\end{proof}

On \(\mathcal U\cap\{u<r\}\) we have \(\ell,\eta>0\), so only the first linear factor in
\eqref{eq:residual} can vanish, and the equilibrium equation is equivalent to
\begin{equation}\label{eq:interp}
 \frac{B(u)}{D(u)}=r+\frac{2(r-u)}{x-1},\qquad x=\sqrt{1+8u}.
\end{equation}
The condition \(A=(1+u)D\) says \(t_i=D_i+D_{i-1}\), and \(\ST=\ET+\FT\) says that there is as
much substrate as enzyme. They serve only to find a centre for the construction;
Theorem~\ref{thm:main}(c) removes both.

\subsection{The recurrence}
Put \(m=n-1\). Let \(x_0,\dots,x_{2m}\) be real numbers greater than \(1\), \emph{not
necessarily distinct}, and write
\begin{equation}\label{eq:ux}
 u(x)=\frac{x^2-1}{8},\qquad u_j=u(x_j),\qquad \Pi(x)=\prod_{j=0}^{2m}(x-x_j).
\end{equation}
Start with \(N=1\) and \(J=(x_0-1)/2\), and for \(k=0,\dots,m-1\), with
\((a,b)=(x_{2k+1},x_{2k+2})\), replace \(N,J\) simultaneously by
\begin{equation}\label{eq:recurrence}
\begin{aligned}
 N_{\mathrm{new}}&=\bigl[8u+(a+1)(b+1)\bigr]N+2(a+b)J,\\
 J_{\mathrm{new}}&=4(a+b)\,uN+\bigl[8u+(a-1)(b-1)\bigr]J .
\end{aligned}
\end{equation}

\begin{lemma}\label{lem:recurrence}
After \(m\) steps \(N=\sum_{i=0}^m N_iu^i\) and \(J=\sum_{i=0}^m J_iu^i\) have degree \(m\),
all coefficients \(N_i,J_i\) \((0\le i\le m)\) are positive, and
\begin{equation}\label{eq:product}
 (x-1)\,N(u(x))-2\,J(u(x))=\Pi(x).
\end{equation}
Moreover
\begin{equation}\label{eq:leading}
 N_m=8^m,\qquad J_m=8^m\,\frac{\sum_jx_j-1}{2},\qquad J_0=J(0)=\frac12\prod_{j}(x_j-1).
\end{equation}
\end{lemma}
\begin{proof}
Let \(N,J\) have positive coefficients in degrees \(0,\dots,k\) and none above. The
coefficient of \(u^i\) in \(N_{\mathrm{new}}\) is \(8N_{i-1}+(a+1)(b+1)N_i+2(a+b)J_i\) and in
\(J_{\mathrm{new}}\) it is \(4(a+b)N_{i-1}+8J_{i-1}+(a-1)(b-1)J_i\), with
\(N_{-1}=J_{-1}=N_{k+1}=J_{k+1}=0\). Since \(a,b>1\), both are positive for \(0\le i\le k+1\).
The leading coefficients obey \(N_k\mapsto8N_k\) and \(J_k\mapsto4(a+b)N_k+8J_k\), which gives
the first two formulas in \eqref{eq:leading}. For \eqref{eq:product}, the initial pair gives
\(x-x_0\), and a direct expansion using \(8u=x^2-1\) shows
\[
 (x-1)N_{\mathrm{new}}-2J_{\mathrm{new}}
 =(x-a)(x-b)\bigl[(x-1)N-2J\bigr]
\]
at \(u=u(x)\): the coefficient of \(N\) is
\((x-1)(x^2+ab+a+b)-(a+b)(x^2-1)=(x-1)(x-a)(x-b)\), and the coefficient of \(J\) is
\(2(a+b)(x-1)-2(x^2+ab-a-b)=-2(x-a)(x-b)\). Setting \(x=1\) in \eqref{eq:product} gives \(J_0\).
\end{proof}

Writing \(\Pi(x)=\Pi_{\mathrm e}(x^2)+x\,\Pi_{\mathrm o}(x^2)\) and replacing \(x\) by \(-x\) in
\eqref{eq:product} shows that, as functions of \(x^2=1+8u\),
\begin{equation}\label{eq:evenodd}
 N=\Pi_{\mathrm o},\qquad J=-\tfrac12\bigl(\Pi_{\mathrm e}+\Pi_{\mathrm o}\bigr).
\end{equation}
So \(N\) and \(J\) depend only on the multiset \(\{x_j\}\); the recurrence is a convenient way
to compute them and to see that their coefficients are positive.

\subsection{Interlacing and coefficient ratios}
\begin{lemma}[Interlacing]\label{lem:interlace}
Let \(m\ge1\). The polynomials \(N\) and \(J\) have only real, simple roots
\(\nu_m<\dots<\nu_1\) and \(\zeta_m<\dots<\zeta_1\), and
\begin{equation}\label{eq:interlace}
 \nu_m<\zeta_m<\dots<\nu_2<\zeta_2<\nu_1<\zeta_1<0 .
\end{equation}
Consequently
\begin{equation}\label{eq:pf}
 \frac{N(u)}{J(u)}=\frac{N_m}{J_m}+\sum_{k=1}^m\frac{\sigma_k}{u-\zeta_k},\qquad
 \frac{J(u)}{N(u)}=\frac{J_m}{N_m}-\sum_{k=1}^m\frac{\varsigma_k}{u-\nu_k},\qquad
 \sigma_k,\varsigma_k>0 .
\end{equation}
\end{lemma}
\begin{proof}
For \(y>0\) evaluate \eqref{eq:product} at \(x=iy\), where \(u=-(1+y^2)/8\) is real and runs
over \((-\infty,-\tfrac18)\). Each factor is
\(iy-x_j=-(x_j^2+y^2)^{1/2}e^{-i\vartheta_j}\) with \(\vartheta_j=\arctan(y/x_j)\). The number
of factors is odd, so \(\Pi(iy)=\varrho\,(-\cos\vartheta+i\sin\vartheta)\) with
\(\varrho=\prod_j(x_j^2+y^2)^{1/2}>0\) and \(\vartheta(y)=\sum_j\vartheta_j\). Comparing real and
imaginary parts in \((iy-1)N-2J=\Pi(iy)\) gives
\begin{equation}\label{eq:trig}
 N=\frac{\varrho\sin\vartheta}{y},\qquad
 J=\frac{\varrho}{2}\Bigl(\cos\vartheta-\frac{\sin\vartheta}{y}\Bigr).
\end{equation}
The function \(\vartheta\) increases strictly and continuously from \(0\) to \((2m+1)\pi/2\),
whether or not the \(x_j\) are distinct, so it takes each of the values \(\pi,2\pi,\dots,m\pi\)
exactly once, at points \(y_1<\dots<y_m\). These give \(m\) distinct roots
\(\nu_k=-(1+y_k^2)/8<-\tfrac18\) of \(N\), which are all its roots since \(\deg N=m\). By
\eqref{eq:trig}, \(J(\nu_k)=(-1)^k\varrho/2\). Hence \(J\) changes sign on each interval
\((\nu_{k+1},\nu_k)\), and since \(J(\nu_1)<0<J(0)\), also on \((\nu_1,0)\). This accounts for
all \(m\) roots of \(J\) and proves \eqref{eq:interlace}. All roots are simple, so the
partial-fraction expansions \eqref{eq:pf} exist with \(\sigma_k=N(\zeta_k)/J'(\zeta_k)\) and
\(\varsigma_k=-J(\nu_k)/N'(\nu_k)\). Exactly \(k-1\) roots of \(N\) and exactly \(k-1\) roots of \(J\)
exceed \(\zeta_k\), and both leading coefficients are positive; so \(N(\zeta_k)\) and
\(J'(\zeta_k)\) both have sign \((-1)^{k-1}\) and \(\sigma_k>0\). Similarly \(N'(\nu_k)\) has
sign \((-1)^{k-1}\) and \(J(\nu_k)\) has sign \((-1)^k\), so \(\varsigma_k>0\).
\end{proof}

\begin{lemma}[Coefficient ratios]\label{lem:ratios}
With \(J_{-1}=0\), the sequence \(J_{i-1}/J_i\) \((0\le i\le m)\) is strictly increasing, and
the sequence \(J_i/N_i\) \((0\le i\le m)\) is nondecreasing.
\end{lemma}
\begin{proof}
A polynomial with positive coefficients and only real roots satisfies Newton's inequalities
\(J_i^2>J_{i-1}J_{i+1}\) \cite[\S2.22]{HLP}; this is the first claim. For the second, let
\(m\ge1\) and \(Q_k=N/(u-\nu_k)=\sum_{i=0}^{m-1}q_iu^i\). Its roots are real and negative, so
\(q_i>0\) and, by Newton's inequalities again, \(q_{i-1}/q_i\) is nondecreasing in \(i\). With
\(q_{-1}=q_m=0\) we have \(N_i=|\nu_k|q_i+q_{i-1}\), hence
\[
 \frac{q_i}{N_i}=\frac{1}{|\nu_k|+q_{i-1}/q_i}\quad(0\le i\le m-1),\qquad \frac{q_m}{N_m}=0,
\]
which is nonincreasing in \(i\). By \eqref{eq:pf},
\(J=(J_m/N_m)N-\sum_k\varsigma_kQ_k\) with \(\varsigma_k>0\), so \(J_i/N_i\) is nondecreasing.
\end{proof}

\subsection{Conversion and coefficient positivity}
For \(r>0\) write \(N_r=N(r)\), \(J_r=J(r)\) and define
\begin{equation}\label{eq:conversion}
 Q(u)=\frac{u\,J_rN(u)-r\,N_rJ(u)}{u-r},\qquad
 D_{\raw}=4Q,\qquad
 B_{\raw}=rD_{\raw}+(4rN_r-2J_r)N-2J_rJ .
\end{equation}
The numerator of \(Q\) vanishes at \(u=r\), so \(Q\) is a polynomial, of degree \(m\) with
leading coefficient \(J_rN_m\). Comparing coefficients in \((u-r)Q=uJ_rN-rN_rJ\) gives
\begin{equation}\label{eq:qformula}
 Q_0=N_rJ_0,\qquad Q_i=N_rJ_i+\sum_{k<i}r^{k-i}\bigl(N_rJ_k-J_rN_k\bigr)\quad(1\le i\le m).
\end{equation}

\begin{lemma}\label{lem:identity}
Let \(\cs=J_r\) and \(\es=4rN_r-J_r\). Then, with \(u=u(x)\),
\begin{equation}\label{eq:interpolation}
 (x-1)\bigl(B_{\raw}(u)-rD_{\raw}(u)\bigr)+2(u-r)D_{\raw}(u)=(\cs x+\es)\,\Pi(x).
\end{equation}
\end{lemma}
\begin{proof}
We have \(B_{\raw}-rD_{\raw}=(\es-\cs)N-2\cs J\) and
\(2(u-r)D_{\raw}=8\cs uN-2(\cs+\es)J\), because \(8rN_r=2(\cs+\es)\). With \(8u=x^2-1\) the
left side becomes
\((x-1)\bigl[\es-\cs+\cs(x+1)\bigr]N-2(\cs x+\es)J=(\cs x+\es)\bigl[(x-1)N-2J\bigr]\), and
\eqref{eq:product} applies.
\end{proof}

For \(1\le k\le m\) let \(J^{(k)}(u)=J(u)/(u-\zeta_k)\). Since all \(\zeta_i\) are negative,
\(J^{(k)}\) has positive coefficients in degrees \(0,\dots,m-1\), and \(J^{(k)}(r)>0\) for
\(r>0\).

\begin{theorem}[Coefficient positivity]\label{thm:positivity}
\textup{(i)} For every \(r>0\), all coefficients of \(Q\) and \(D_{\raw}\) in degrees
\(0,\dots,m\) are positive.

\textup{(ii)} Let \(\rho=\sqrt{1+8r}\). If \(\rho>\max_jx_j\) and
\(4r+\rho\ge\sum_jx_j\), then all coefficients of \(B_{\raw}\) in degrees \(0,\dots,m\) are
positive.
\end{theorem}
\begin{proof}
For \(m=0\) we have \(N=1\), \(J=J_0\), \(Q=J_0\) and \(B_{\raw}=2(1+J_0)(2r-J_0)\). Here
\(J_0=(x_0-1)/2\), and \(\rho>x_0\) gives \(2r=(\rho^2-1)/4>(x_0-1)(x_0+1)/4>J_0\). Let
\(m\ge1\) and write \(\alpha_0=N_m/J_m\).

(i) By \eqref{eq:pf}, \(uN/J=\alpha_0u+\sum_k\sigma_k+\sum_k\sigma_k\zeta_k/(u-\zeta_k)\). Hence
\begin{align*}
 Q(u)&=J(u)J_r\,\frac{\frac{uN(u)}{J(u)}-\frac{rN_r}{J_r}}{u-r}
 =J(u)J_r\Bigl(\alpha_0+\sum_k\frac{\sigma_k|\zeta_k|}{(u-\zeta_k)(r-\zeta_k)}\Bigr)\\
 &=\alpha_0J_r\,J(u)+\sum_{k=1}^m\sigma_k|\zeta_k|\,J^{(k)}(r)\,J^{(k)}(u).
\end{align*}
Every term has nonnegative coefficients and the first has positive coefficients in all
degrees \(0,\dots,m\).

(ii) From \eqref{eq:pf}, \(N=\alpha_0J+\sum_k\sigma_kJ^{(k)}\). Substituting this and the
formula for \(Q\) into \(B_{\raw}=4rQ+4rN_rN-2J_r(N+J)\) and using
\(J_r=(r+|\zeta_k|)J^{(k)}(r)\) gives
\begin{equation}\label{eq:Bdecomp}
 B_{\raw}=\alpha_0\Bigl[4rN_r-J_r\Bigl(\sum_jx_j+1-4r\Bigr)\Bigr]J
 +\sum_{k=1}^m\sigma_k\Bigl[4rN_r-2J_r\Bigl(1-\frac{2r|\zeta_k|}{r+|\zeta_k|}\Bigr)\Bigr]J^{(k)} ,
\end{equation}
where we used \(2+2/\alpha_0=\sum_jx_j+1\), from \eqref{eq:leading}. Evaluating
\eqref{eq:product} at \(x=\rho\), where \(u=r\), gives \((\rho-1)N_r-2J_r=\Pi(\rho)>0\), that is,
\begin{equation}\label{eq:Jr}
 J_r<\tfrac12(\rho-1)N_r .
\end{equation}
Since \(4r=(\rho-1)(\rho+1)/2\), each bracket in the sum in \eqref{eq:Bdecomp} is at least
\(4rN_r-2J_r>\tfrac12(\rho-1)^2N_r>0\). If \(\sum_jx_j+1-4r\le0\), the first bracket is
positive as well. Otherwise, by \eqref{eq:Jr}, it exceeds
\[
 \tfrac12(\rho-1)N_r\Bigl[(\rho+1)-\Bigl(\sum_jx_j+1-4r\Bigr)\Bigr]
 =\tfrac12(\rho-1)N_r\Bigl(4r+\rho-\sum_jx_j\Bigr)\ \ge0 .
\]
Thus \(B_{\raw}\) is a combination of \(J,J^{(1)},\dots,J^{(m)}\) with positive coefficients,
and \(J\) has positive coefficients in all degrees \(0,\dots,m\).
\end{proof}

The hypotheses of (ii) hold for all sufficiently large \(r\). For the equally spaced roots
\(x_j=j+2\) they hold for every \(r>(4n^2-1)/8\), because then \(\rho>2n=\max_jx_j\) and
\((\rho+1)^2>4n^2+2n=2(1+\sum_jx_j)\), which is the second condition. For the coalesced roots
\(x_j=3\) used below they hold for every \(r\ge3n/2\).

\subsection{Common totals, the count, and nondegeneracy}
Fix \(r\) as in Theorem~\ref{thm:positivity}(ii); note that \(\rho>x_j\) if and only if
\(r>u_j\). Normalize and choose the totals by
\begin{equation}\label{eq:normalize}
 D=\frac{D_{\raw}}{D_{\raw}(0)},\qquad B=\frac{B_{\raw}}{D_{\raw}(0)},\qquad
 A=(1+u)D,\qquad (\ET,\FT,\ST)=(2r,\,2,\,2r+2).
\end{equation}
Then \(A(0)=D(0)=1\), all coefficients of \(A,B,D\) in the ranges required by \eqref{eq:coeff}
are positive, and Proposition~\ref{prop:param}(ii) provides rate constants for every choice of
the kinetics \((b_i,\beta_i,\gamma_i)\). By Lemma~\ref{lem:identity}, for every \(j\),
\begin{equation}\label{eq:ratio}
 \frac{B(u_j)}{D(u_j)}=r+\frac{2(r-u_j)}{x_j-1}>r .
\end{equation}
For \(x=x_j\), \(u=u_j\) define
\begin{equation}\label{eq:sf}
 s=\frac{x-1}{2uD(u)},\qquad f=\frac{4}{x+1},\qquad\text{so that}\qquad E=uf=\frac{x-1}{2},
\end{equation}
and let \(X_j\) be the point \eqref{eq:states} with these \((u,s,f)\). It is a positive
equilibrium by Proposition~\ref{prop:param}(i). Its totals follow from \eqref{eq:totalcoords}:
first \(suD=(x-1)/2\), so \(\FT=f(1+suD)=2\); next, by \eqref{eq:ratio},
\(suB=\tfrac12r(x-1)+(r-u)\), so
\(\ET=uf+f\bigl(\tfrac12r(x-1)+r-u\bigr)=\tfrac12 f\,r\,(x+1)=2r\); finally
\(sA=(1+u)(x-1)/(2u)=4(1+u)/(x+1)=(1+u)f\), that is,
\begin{equation}\label{eq:freeinventory}
 \sum_{i=0}^nS_i=E+F ,
\end{equation}
and since the bound substrate is \(\ET+\FT-E-F\), we get \(\ST=\ET+\FT=2r+2\).

\begin{proposition}[Slopes]\label{prop:slope}
For every \(j\),
\begin{equation}\label{eq:slope}
 \Phi'(u_j)=-\frac{16\,(1+u_j)\,(\cs x_j+\es)}{(r-u_j)\,(x_j+1)^2\,D_{\raw}(u_j)}\;\Pi'(x_j),
 \qquad \cs x_j+\es=4rN_r+(x_j-1)J_r>0 .
\end{equation}
In particular, if the \(x_j\) are distinct and ordered increasingly, then \(\Phi'(u_j)\ne0\) and
\(\sign\Phi'(u_j)=(-1)^{j+1}\); if \(x_j\) is a repeated root of \(\Pi\), then \(\Phi'(u_j)=0\).
\end{proposition}
\begin{proof}
By \eqref{eq:ratio}, \(u_j\) lies in \(\mathcal U\cap\{u<r\}\), where Lemma~\ref{lem:factor}
applies. Dividing \eqref{eq:interpolation} by \((x-1)D_{\raw}(u)\) gives
\[
 \Psi_1(u):=\ell-\frac{2\eta}{x-1}=\frac{(\cs x+\es)\,\Pi(x)}{(x-1)\,D_{\raw}(u)},
\]
which vanishes at \(u_j\), and \(\Phi-\ST=\Psi_2\Psi_1\) with
\(\Psi_2=-2(1+u)\bigl(\ell+2\eta/(x+1)\bigr)/\bigl(\ell(\ell+\eta)\bigr)\). Hence
\(\Phi'(u_j)=\Psi_2(u_j)\Psi_1'(u_j)\). At \(u_j\) we have \(\ell=2\eta/(x-1)\), so
\(\ell+2\eta/(x+1)=4\eta x/(x^2-1)\), \(\ell+\eta=\eta(x+1)/(x-1)\), and
\(\Psi_2(u_j)=-4(1+u)x(x-1)/\bigl(\eta(x+1)^2\bigr)\). Since \(\Pi(x_j)=0\) and \(dx/du=4/x\),
\(\Psi_1'(u_j)=4(\cs x+\es)\Pi'(x)/\bigl(x(x-1)D_{\raw}(u)\bigr)\). Multiplying gives
\eqref{eq:slope}. All factors other than \(\Pi'(x_j)=\prod_{k\ne j}(x_j-x_k)\) are positive, and
\(\Pi'(x_j)\) has sign \((-1)^{2m-j}=(-1)^j\) when the roots are distinct and ordered
increasingly.
\end{proof}

\begin{theorem}[Maximal equilibrium count]\label{thm:count}
Let \(n\ge1\), let \(x_0<\dots<x_{2n-2}\) be real numbers greater than \(1\), and let \(r\)
satisfy \(r>\max_ju_j\) and \(4r+\sqrt{1+8r}\ge\sum_jx_j\). For every choice of positive
kinetics \((b_i,\beta_i,\gamma_i)\), the rate constants \eqref{eq:ratesgeneral} formed from
\eqref{eq:normalize} are positive, and the system \eqref{eq:network} with totals
\((2r,2,2r+2)\) has exactly the \(2n-1\) positive equilibria \(X_0,\dots,X_{2n-2}\), with
\(E/F=u_j\) at \(X_j\). All of them are nondegenerate, and at the \(n-1\) equilibria with odd
index the restricted Jacobian has an odd number of positive real eigenvalues. If the \(x_j\),
\(r\) and the kinetics are rational, then so are all rate constants and all concentrations.
Consequently \(\capE=2n-1\), and exactly \(2n-1\) nondegenerate positive equilibria occur on a
nonempty open subset of \(\R^{6n+3}_{>0}\).
\end{theorem}
\begin{proof}
The points \(X_j\) are positive equilibria of one and the same system with the same totals, and
they are distinct because their ratios \(E/F=u_j\) are. Theorem~\ref{thm:WS} excludes any
others. Nondegeneracy follows from Propositions~\ref{prop:slope} and~\ref{prop:det}, the
statement about odd indices from Corollary~\ref{cor:parity}, and rationality from the
formulas. For the open set, apply the implicit function theorem to the \(3n\) equilibrium
equations in the coordinates \(y\), whose derivative is \(\JS\): each \(X_j\) continues to a
smooth family of nondegenerate positive equilibria for parameters near the constructed ones,
the continued equilibria remain distinct, and Theorem~\ref{thm:WS} again excludes others. The
perturbed parameters need not satisfy \(A=(1+u)D\) or \(\ST=\ET+\FT\).
\end{proof}

The map \(x\mapsto(x^2-1)/8\) is a bijection from \((1,\infty)\) to \((0,\infty)\), so the ratios
\(E/F\) at the \(2n-1\) equilibria can be prescribed arbitrarily. Whether the even-indexed
equilibria are stable depends on the kinetics, which the equilibrium equations do not see. The
rest of the proof of Theorem~\ref{thm:main} consists in choosing the kinetics.