* [65] A. Dabholkar, N. Drukker, and J. Gomes, “Localization in supergravity and quantum AdS4/CFT3 holography,” JHEP 10 (2014) 090, arXiv:1406.0505 [hep-th].
  * [66] A. G. Lezcano, J. Hong, J. T. Liu, and L. A. P. Zayas, “The Bethe-Ansatz approach to the $\mathcal{N}{=}4$ superconformal index at finite rank,” arXiv:2101.12233 [hep-th].
  * [67] F. Benini and G. Rizi, “Superconformal index of low-rank gauge theories via the Bethe Ansatz,” JHEP 05 (2021) 061, arXiv:2102.03638 [hep-th].
  * [68] F. Benini and E. Milan, “A Bethe Ansatz type formula for the superconformal index,” Commun. Math. Phys. 376 (2020) 1413–1440, arXiv:1811.04107 [hep-th].
  * [69] D. Cassani and L. Papini, “The BPS limit of rotating AdS black hole thermodynamics,” JHEP 09 (2019) 079, arXiv:1906.10148 [hep-th].
  * [70] N. Bobev, A. M. Charles, and V. S. Min, “Euclidean black saddles and AdS4 black holes,” JHEP 10 (2020) 073, arXiv:2006.01148 [hep-th].
  * [71] I. Biswas, D. Gaiotto, S. Lahiri, and S. Minwalla, “Supersymmetric states of $\mathcal{N}{=}4$ Yang-Mills from giant gravitons,” JHEP 12 (2007) 006, arXiv:hep-th/0606087.
  * [72] R. Arai and Y. Imamura, “Finite $N$ Corrections to the Superconformal Index of S-fold Theories,” PTEP 2019 no. 8, (2019) 083B04, arXiv:1904.09776 [hep-th].
  * [73] R. Arai, S. Fujiwara, Y. Imamura, and T. Mori, “Schur index of the $\mathcal{N}{=}4$ $U(N)$ supersymmetric Yang-Mills theory via the AdS/CFT correspondence,” Phys. Rev. D 101 no. 8, (2020) 086017, arXiv:2001.11667 [hep-th].
  * [74] E. Witten, “Constraints on Supersymmetry Breaking,” Nucl. Phys. B 202 (1982) 253.
  * [75] B. Assel, D. Cassani, L. Di Pietro, Z. Komargodski, J. Lorenzen, and D. Martelli, “The Casimir Energy in Curved Space and its Supersymmetric Counterpart,” JHEP 07 (2015) 043, arXiv:1503.05537 [hep-th].
  * [76] I. Papadimitriou, “Supercurrent anomalies in 4d SCFTs,” JHEP 07 (2017) 038, arXiv:1703.04299 [hep-th].
  * [77] O. S. An, “Anomaly-corrected supersymmetry algebra and supersymmetric holographic renormalization,” JHEP 12 (2017) 107, arXiv:1703.09607 [hep-th].
  * [78] I. Papadimitriou, “Supersymmetry anomalies in $\mathcal{N}{=}1$ conformal supergravity,” JHEP 04 (2019) 040, arXiv:1902.06717 [hep-th].
  * [79] I. Papadimitriou, “Supersymmetry anomalies in new minimal supergravity,” JHEP 09 (2019) 039, arXiv:1904.00347 [hep-th].
  * [80] S. M. Kuzenko, A. Schwimmer, and S. Theisen, “Comments on Anomalies in Supersymmetric Theories,” J. Phys. A 53 (2020) 064003, arXiv:1909.07084 [hep-th].
  * [81] G. Katsianis, I. Papadimitriou, K. Skenderis, and M. Taylor, “Supersymmetry anomaly in the superconformal Wess-Zumino model,” JHEP 02 (2021) 209, arXiv:2011.09506 [hep-th].
  * [82] A. Bzowski, G. Festuccia, and V. Procházka, “Consistency of supersymmetric ’t Hooft anomalies,” JHEP 02 (2021) 225, arXiv:2011.09978 [hep-th].
  * [83] C. Closset, L. Di Pietro, and H. Kim, “’t Hooft anomalies and the holomorphy of supersymmetric partition functions,” JHEP 08 (2019) 035, arXiv:1905.05722 [hep-th].
  * [84] B. Sundborg, “The Hagedorn transition, deconfinement and $\mathcal{N}{=}4$ SYM theory,” Nucl. Phys. B 573 (2000) 349–363, arXiv:hep-th/9908001.
  * [85] O. Aharony, J. Marsano, S. Minwalla, K. Papadodimas, and M. Van Raamsdonk, “The Hagedorn - deconfinement phase transition in weakly coupled large $N$ gauge theories,” Adv. Theor. Math. Phys. 8 (2004) 603–696, arXiv:hep-th/0310285.
  * [86] C. Closset, H. Kim, and B. Willett, “$\mathcal{N}{=}1$ supersymmetric indices and the four-dimensional $A$-model,” JHEP 08 (2017) 090, arXiv:1707.05774 [hep-th].
  * [87] M. Cvetic, H. Lu, and C. Pope, “Charged Kerr-de Sitter black holes in five dimensions,” Phys. Lett. B 598 (2004) 273–278, arXiv:hep-th/0406196.
  * [88] Z. W. Chong, M. Cvetic, H. Lu, and C. N. Pope, “General non-extremal rotating black holes in minimal five-dimensional gauged supergravity,” Phys. Rev. Lett. 95 (2005) 161301, arXiv:hep-th/0506029 [hep-th].
  * [89] M. Cvetic, M. J. Duff, P. Hoxha, J. T. Liu, H. Lu, J. X. Lu, R. Martinez-Acosta, C. N. Pope, H. Sati, and T. A. Tran, “Embedding AdS black holes in ten-dimensions and eleven-dimensions,” Nucl. Phys. B558 (1999) 96–126, arXiv:hep-th/9903214 [hep-th].
  * [90] G. Felder and A. Varchenko, “The elliptic gamma function and $SL(3,\mathbb{Z})\ltimes\mathbb{Z}^{3}$,” Adv. Math. 156 (2000) 44–76, arXiv:math/9907061 [math.QA].
  * [91] J. B. Gutowski and H. S. Reall, “Supersymmetric AdS5 black holes,” JHEP 02 (2004) 006, arXiv:hep-th/0401042 [hep-th].
  * [92] O. Madden and S. F. Ross, “On uniqueness of charged Kerr-AdS black holes in five dimensions,” Class. Quant. Grav. 22 (2005) 515–524, arXiv:hep-th/0409188.
  * [93] D. Cassani, J. Lorenzen, and D. Martelli, “Comments on supersymmetric solutions of minimal gauged supergravity in five dimensions,” Class. Quant. Grav. 33 no. 11, (2016) 115013, arXiv:1510.01380 [hep-th].
  * [94] A. Chamblin, R. Emparan, C. V. Johnson, and R. C. Myers, “Charged AdS black holes and catastrophic holography,” Phys. Rev. D 60 (1999) 064018, arXiv:hep-th/9902170.
  * [95] A. Buchel and J. T. Liu, “Gauged supergravity from type IIB string theory on $Y^{p,q}$ manifolds,” Nucl. Phys. B 771 (2007) 93–112, arXiv:hep-th/0608002.
  * [96] M. Cederwall, A. von Gussich, B. E. Nilsson, and A. Westerberg, “The Dirichlet super three-brane in ten-dimensional type IIB supergravity,” Nucl. Phys. B 490 (1997) 163–178, arXiv:hep-th/9610148.
  * [97] M. Cederwall, A. von Gussich, B. E. Nilsson, P. Sundell, and A. Westerberg, “The Dirichlet super $p$-branes in ten-dimensional type IIA and IIB supergravity,” Nucl. Phys. B 490 (1997) 179–201, arXiv:hep-th/9611159.
  * [98] E. Bergshoeff and P. Townsend, “Super D-branes,” Nucl. Phys. B 490 (1997) 145–162, arXiv:hep-th/9611173.
  * [99] E. Bergshoeff, R. Kallosh, T. Ortin, and G. Papadopoulos, “Kappa symmetry, supersymmetry and intersecting branes,” Nucl. Phys. B 502 (1997) 149–169, arXiv:hep-th/9705040.
  * [100] A. Mikhailov, “Giant gravitons from holomorphic surfaces,” JHEP 11 (2000) 027, arXiv:hep-th/0010206.
  * [101] M. Cvetic, H. Lu, and C. Pope, “Charged rotating black holes in five dimensional $U(1)^{3}$ gauged $\mathcal{N}{=}2$ supergravity,” Phys. Rev. D 70 (2004) 081502, arXiv:hep-th/0407058.
  * [102] M. Cvetic, G. Gibbons, H. Lu, and C. Pope, “Rotating black holes in gauged supergravities: Thermodynamics, supersymmetric limits, topological solitons and time machines,” arXiv:hep-th/0504080.
  * [103] S.-Q. Wu, “General Nonextremal Rotating Charged AdS Black Holes in Five-dimensional $U(1)^{3}$ Gauged Supergravity: A Simple Construction Method,” Phys. Lett. B 707 (2012) 286–291, arXiv:1108.4159 [hep-th].
  * [104] F. Benini and N. Bobev, “Two-dimensional SCFTs from wrapped branes and $c$-extremization,” JHEP 06 (2013) 005, arXiv:1302.4451 [hep-th].
  * [105] A. Sinha, J. Sonner, and N. V. Suryanarayana, “At the horizon of a supersymmetric AdS5 black hole: Isometries and half-BPS giants,” JHEP 01 (2007) 087, arXiv:hep-th/0610002.
  * [106] A. Sinha and J. Sonner, “Black hole giants,” JHEP 08 (2007) 006, arXiv:0705.0373 [hep-th].


Despite their noted potential in polynomial-system solving, there are few concrete examples of Newton-Okounkov bodies
arising from applications.
Accordingly, in this paper, we introduce a new application of Newton-Okounkov body theory to the study of chemical reaction networks, and compute several examples.
An important invariant of
a chemical reaction network is its maximum number of positive steady states
, which is realized as the maximum number of positive real roots of a parametrized polynomial system.
Here, we introduce a new upper bound on this number, namely the
`Newton-Okounkov body bound' of
a chemical reaction network.
Through explicit examples, we show that the
Newton-Okounkov body
bound of
a network gives a good upper bound on its maximum number of positive steady states.
We also compare this
Newton-Okounkov body bound
to a related upper bound, namely the mixed volume of a chemical reaction network, and find that it often achieves better bounds.

§ INTRODUCTION

A Newton-Okounkov body associated to a projective variety is a convex body
that encodes important information about the variety.
Newton-Okounkov bodies associated to projective varieties can be thought of as a vast generalization of Newton polytopes associated to projective toric varieties.
One important parallel
is that the volume of the convex body and is directly related to the degree of the associated projective variety, which then conveys the number of solutions to a general polynomial system.

In contrast to Newton polytopes of projective toric varieties,
Newton-Okounkov bodies associated to arbitrary projective varieties
require an involved construction.
Part of the contribution of this paper is a minimal, yet explicit explanation of the setting and construction of a Newton-Okounkov body.
We refer the interested reader to the foundational works <cit.> for technical details of Newton-Okounkov body theory.

With an eye towards applications, our focus is on the
of Newton-Okounkov bodies in intersection theory, specifically to compute bounds
on the number of isolated solutions to polynomial systems.
The starting data is a finite dimensional vector subspace of rational functions $V$ contained in the coordinate ring of
a $d$-dimensional complex algebraic variety $X$.
Given such a vector space $V$, Kaveh and Khovanskii <cit.> defined the (birationally-invariant) self-intersection index of $V$ to be the number of “effective” solutions in $X$ to a system of $d$ general equations $f_1=\dots=f_d=0$, where each $f_i$ is a general rational function in $V$.
Here we use “effective” solution to mean that the intersection index counts only solutions in the smooth locus of $X$, and outside the base locus (see Definition <ref>).

As an example,
the classical Kushnirenko Theorem (Proposition <ref>) computes the self-intersection index of a subspace $V$ of Laurent polynomials.
More explicitly, consider a subspace $V \subset \C[x_1^{\pm 1},\dots,x_d^{\pm 1}]$ spanned by Laurent monomials with
contained in some finite subset $A\subset \Z^d$.
The self-intersection index of $V$ is equal to the normalized Euclidean volume of the common Newton polytope $\Delta:=\conv{A}$.
In other words, a general system of $d$
polynomials in $V$
with Newton polytope $\Delta=\conv{A}$
has $d!\,\vol{\Delta}$ effective solutions (all in $\C^*$).
The Newton polytope is equal to the Newton-Okounkov body in this setting.

Generalizing from vector spaces spanned by monomials, Kaveh and Khovanskii <cit.> subsequently gave a formula for the self-intersection index for subspaces $V$ spanned by rational functions.
That is, the number of effective solutions to a general system of $d$ equations
drawn from a
subspace of rational functions $V$
is proportional to
the volume of an associated Newton-Okounkov body $\Delta$
; see Proposition <ref>.

The main challenge of Kaveh and Khovanskii's formula is the computation of an associated Newton-Okounkov body.
Such convex bodies are not necessarily polytopes, nor are the associated semigroups even finitely generated.
However, the Newton-Okounkov body will be a polytope if there exists an associated finite Khovanskii basis <cit.>.
The formalization of Khovanskii bases, which generalize sagbi bases <cit.>, is attributed to Kaveh and Manon, who described Khovanskii bases as the computational and algorithmic side of the theory of Newton-Okounkov bodies <cit.>.
This computational side of Newton-Okounkov bodies has been explored in
polynomial-system solving, such as in <cit.>.
these works focused on computing all solutions to a given polynomial system, our work is concerned with
computing bounds on the number of solutions to polynomial systems arising from applications.
Specifically, our application comes from counting steady states of a chemical reaction network.
Concretely, our focus is on
computing a
self-intersection index associated to a given chemical reaction network, and assessing the resulting bound on the maximum number of the network's steady states.

A chemical reaction network is a model of the interactions of chemical species.
classical assumption (mass-action kinetics),
are governed by an autonomous system of parametrized polynomial ordinary differential equations – see equation (<ref>).
Then, finding steady states of a chemical reaction network amounts to solving a system of (parametrized) polynomial equations.
The capacity for multiple steady states (called multistationarity) has been studied extensively in recent years, with algebraic-geometric methods at the forefront of new results <cit.>.
After deciding multistationarity, actually determining an arbitrary reaction network's maximum number of observable steady states is still an important open problem.
One strategy is to produce good upper bounds on the maximum number of observable steady states.

Accordingly, Obatake, in prior work with Shiu, Tang, and Torres, introduced the mixed volume of a chemical reaction network
and Gross and Hill introduced a related steady-state degree.
We refer the reader to <cit.> for a discussion of the several “flavors” of mixed volume theory applied to chemical reaction networks.
parameter-free numerical invariants gave good bounds on the maximum number of steady states for several
families of networks <cit.>.
However, the mixed-volume bound on a reaction network (and the related steady-state degree) is not always tight: the maximum number of observable, positive steady states can be far less than the mixed volume.
Here, through our examples, we show that the self-intersection index can give a better, tighter bound than was previously possible.

In this paper, we explain the construction of a Newton-Okounkov body and the self-intersection index of a polynomial vector space, with an eye towards applying the theory to chemical reaction networks.
In an effort to make the theory of Newton-Okounkov bodies more accessible to an applied audience, we explain background material through a motivating example (the Wnt network, see Example <ref>).
Moreover, our main contribution is one of the first concrete applications of Newton-Okounkov body theory to polynomial-system solving.
Readers interested in computing an upper bound on the number of effective solutions to a (parametrized) polynomial system may proceed to Procedure <ref>.
We emphasize that this procedure is not just for chemical reaction networks, and can be applied to any sparse parametrized polynomial system.

The outline of this paper is as follows.
Section <ref> covers the basics of chemical reaction networks.
Section <ref> explains the Newton-Okounkov body theory and self-intersection index results we will need for our application.
Section <ref> explains our procedure for explicitly computing
self-intersection index of a chemical reaction network.
In Section <ref>, we
apply our procedure to several
chemical reaction networks.
We compare the resulting “Newton-Okounkov body bound” on the maximum number of steady states with both the mixed-volume bound and the actual maximum number of positive real steady states.
The paper concludes with a discussion of the implications of this new tool for chemical reaction network theory and other applications.
We include some conjectures and avenues for future research.
Finally, Appendix <ref> collects the notation introduced in Section <ref>.

Supporting computations are archived in the following GitHub repository: <https://github.com/EliseAWalker/NOBodies-for-ChemicalRxns>

§ CHEMICAL REACTION NETWORKS AND THEIR DYNAMICAL SYSTEMS

Our notation for chemical reaction networks closely matches that of
We briefly review the basics of the mathematical models of chemical reaction networks in this section, and refer the reader to <cit.> and references therein for background.

A (chemical) reaction network $G$  (or network for short) comprises a set of $d$ species $\{X_1, X_2, \dots, X_d\}$ and a set of $m$ reactions:

\[
\alpha_{1j}X_1 +
\alpha_{2j}X_2 +  \dots +
\alpha_{dj}X_d
~ \overset{\kappa_j}{\longrightarrow} ~ %
\beta_{1j}X_1 +
\beta_{2j}X_2 +  \dots +
\beta_{dj}X_d~,
\quad \quad
{\rm for}~
j=1,2, \dots, m~,
\]

where each $\alpha_{ij}$ and $\beta_{ij}$ is a non-negative integer
called a stoichiometric coefficient and each $\kappa_j$ is a nonnegative real number called a reaction rate constant.
The stoichiometric matrix of
denoted by $N$, is the $d\times m$ matrix with
$(i, j)$-entry equal to $\beta_{ij}-\alpha_{ij}$.
Let $p=d-{\rm rank}(N)$.
The stoichiometric subspace, denoted by $S$,
is the image of $N$, that is, $S$ is the vector subspace of $\R^p$ generated by
the columns of $N$.
A conservation-law matrix of $G$, denoted by $W$, is a row-reduced $(p\times d$)-matrix whose rows form a basis of
the orthogonal complement of $S$.
If there exists a choice of $W$ for which every entry is nonnegative and each column
contains at least one nonzero entry (equivalently, if
each species occurs in at least one nonnegative conservation law), then $G$ is conservative.

We denote the concentrations of the species $X_1,X_2, \dots, X_d$ by $x_1, x_2, \dots, x_d$, respectively.
These concentrations, under the assumption of mass-action kinetics, evolve in time according to the following system of
ordinary differential equations:

\begin{equation}\label{eq:ode}
\dot{x}~=~f(x)~:=~N\cdot \begin{pmatrix}
\kappa_1 \, x_1^{\alpha_{11}}
\cdots x_d^{\alpha_{d1}} \\
\kappa_2 \, x_1^{\alpha_{12}}
\cdots x_d^{\alpha_{d2}} \\
\vdots \\
\kappa_m \, x_1^{\alpha_{1m}}
\cdots x_d^{\alpha_{dm}} \\
\end{pmatrix}~,
\end{equation}

where $x=(x_1, x_2, \dots, x_d)$,
and each $\kappa_j \in \mathbb R_{>0}$.
By considering the rate constants as a vector of parameters $\kappa=(\kappa_1, \kappa_2, \dots, \kappa_m)$, we have polynomials $f_{\kappa,i} \in \mathbb Q[\kappa,x]$, for $i=1,2, \dots, d$.
For ease of notation, we often write $f_i$ rather than  $f_{\kappa,i}$.
We refer to the dynamical system (<ref>) obtained from a chemical reaction network governed by mass-action kinetics as a chemical reaction system.

A trajectory $x(t)$ beginning at a
positive vector $x(0)=x^0 \in
\mathbb{R}^d_{> 0}$
remains, for all positive time,
in the following stoichiometric compatibility class with respect to the total-constant vector $c\coloneqq W x^0 \in {\mathbb R}^p$:

\begin{align} \label{eqn:invtPoly}
\scc_c~\coloneqq~ \{x\in {\mathbb R}_{\geq 0}^d \mid Wx=c\}~.
\end{align}

Note that when $d={\rm rank}(N)$, the network has no conservation laws. In this case, the stoichiometric subspace $S$ is $\R^d$, and the only stoichiometric compatibility class (<ref>) is $\R_{\geq 0}^d$.

§.§ Steady states

As summarized in the introduction,
one goal of
chemical reaction network theory is
to understand the fixed points,
called steady states,
the differential equations (<ref>) defining chemical reaction systems.

A steady state of (<ref>) is a
concentration vector
$x^* \in \mathbb{R}_{\geq 0}^d$ at which the
right-hand sides of the
chemical reaction system (<ref>)
vanish: $f(x^*) =0$.
We distinguish between positive steady states $x ^* \in \mathbb{R}^d_{> 0}$ and boundary steady states
$x^*\in {\mathbb R}_{\geq 0}^d\backslash {\mathbb R}_{>0}^d$.
In particular, a boundary steady state $x^*\not\in (\C^*)^d$.
A network exhibits multistationarity
if there exists a positive rate-constant vector $\kappa \in \mathbb{R}^m_{>0}$ such that there exist two or more positive steady states of (<ref>) in some stoichiometric compatibility class (<ref>).

To analyze steady states within a stoichiometric compatibility class, we will
use conservation laws
in place of linearly dependent steady-state equations, as follows.
Let $I = \{i_1 < i_2< \dots < i_p\}$ denote the indices of the first nonzero coordinate of the rows of conservation-law matrix $W$.
Consider the function $\F: {\mathbb R}_{\geq 0}^d\rightarrow {\mathbb R}^d$ defined by

\begin{equation}\label{consys}
f_{c,\kappa,i} =\F(x)_i :=
\begin{cases}
f_{i}(x)&~\text{if}~i\not\in I,\\
(Wx-c)_k &~\text{if}~i~=~i_k\in I .
\end{cases}
\end{equation}

We call system (<ref>), the system augmented by conservation laws or just the augmented system. By construction, positive roots of the system of polynomial equations $\F=0$ are precisely the positive steady states of (<ref>) in the stoichiometric compatibility class (<ref>) defined by the total-constant vector $c$.
When the network has no conservation laws, then the augmented system is just the chemical reaction system (<ref>).

We conclude this section with an example of a chemical reaction network, the Wnt network (Figure <ref>), which will serve as the motivating example throughout this paper.

[Motivating Example]
The Wnt network was introduced in <cit.> as a shuttle model of the Wnt signalling pathway, and it was systematically studied using algebraic-geometric methods in <cit.>.
Depicted in Figure <ref>, the network comprises 19 chemical species (whose concentrations are denoted by the 19 variables $x_1,\dots,x_{19}$) participating in 31 reactions, labeled by $k_1,\dots,k_{31}$.
The network is not conservative (the network contains inflow and outflow reactions), however, the network has 5 conserved quantities $c_1,\dots,c_5$, and
it has
5 linear conservation laws.

$\ce{$x_1$ <-->[$k_1$][$k_2$] $x_2$}$                  $\ce{$x_2$ <-->[$k_{26}$][$k_{27}$] $x_3$}$              $\ce{$x_5$ <-->[$k_{28}$][$k_{29}$] $x_7$}$              $\ce{$x_{10}$ <-->[$k_{30}$][$k_{31}$] $x_{11}$}$

$\ce{$x_{10}$ <-->[$k_{13}$][$k_{12}$] \emptyset}$     $\ce{$x_{11}$ ->[$k_{23}$] \emptyset}$     2c$\ce{$x_{11}$ + x_{12} <-->[$k_{24}$][$k_{25}$] $x_{13}$}$

2c$\ce{$x_2$ + $x_4$ <-->[$k_3$][$k_4$] $x_{14}$ ->[$k_5$] $x_2$ + $x_5$}$      2c$\ce{$x_5$ + $x_8$ <-->[$k_6$][$k_7$] $x_{16}$ ->[$k_8$] $x_4$ + $x_8$}$

2c$\ce{$x_3$ + $x_6$ <-->[$k_{14}$][$k_{15}$] $x_{15}$ ->[$k_{16}$] $x_3$ + $x_7$}$      2c$\ce{$x_7$ + $x_9$ <-->[$k_{17}$][$k_{18}$] $x_{17}$ ->[$k_{19}$] $x_6$ + $x_9$}$

2c$\ce{$x_4$ + $x_{10}$ <-->[$k_9$][$k_{10}$] $x_{18}$ ->[$k_{11}$] $x_4$}$      2c$\ce{$x_6$ + $x_{11}$ <-->[$k_{20}$][$k_{21}$] $x_{19}$ ->[$k_{22}$] $x_6$}$

The Wnt network, a formalism of the shuttle model of the Wnt signalling pathway <cit.>.

Let $x(t) = (x_1(t),x_2(t),\dots,x_d(t))\in \Rnn^d$
denote the vector of concentrations of the species as a differentiable function of time $t$.
Let $\kappa:=(k_1,\dots,k_{31})$ denote the vector of rate constants, and $c:=(c_1,\dots,c_5)$ denote the vector of conserved quantities.
In the chemical reaction system (<ref>) of the Wnt network arising from mass-action kinetics, the differential equations $\dot{x}_1,\dot{x}_4,\dot{x}_8,\dot{x}_9,\dot{x}_{12}$ are linear combinations of the remaining $\dot{x}_i$.
This implies the 5 linear conservation equations $f_{c,\kappa}(x)_{1},~f_{c,\kappa}(x)_{4},~f_{c,\kappa}(x)_{8},~f_{c,\kappa}(x)_{9},~f_{c,\kappa}(x)_{12}$ seen in (<ref>).
The augmented system (<ref>) of the Wnt network arising from mass-action kinetics is:

\begin{equation}\label{eq:wnt-aug-sys}
\resizebox{.78\textwidth}{!}{$%
\begin{cases}
f_{2}:=f_{c,\kappa}(x)_2=-k_3 x_2 x_4+k_1 x_1+(-k_2-k_{26}) x_2+k_{27} x_3+(k_4+k_5) x_{14}\\
f_{3}:=f_{c,\kappa}(x)_3=-k_{14} x_3 x_6+k_{26} x_2-k_{27} x_3+(k_{15}+k_{16}) x_{15}\\
f_{5}:=f_{c,\kappa}(x)_5=-k_6 x_5 x_8-k_{28} x_5+k_{29} x_7+k_5 x_{14}+k_7 x_{16}\\
f_{6}:=f_{c,\kappa}(x)_6=-k_{14} x_3 x_6-k_{20} x_6 x_{11}+k_{15} x_{15}+k_{19} x_{17}+(k_{21}+k_{22}) x_{19}\\
f_{7}:=f_{c,\kappa}(x)_7=-k_{17} x_7 x_9+k_{28} x_5-k_{29} x_7+k_{16} x_{15}+k_{18} x_{17}\\
f_{10}:=f_{c,\kappa}(x)_{10}=-k_9 x_4 x_{10}+(-k_{13}-k_{30}) x_{10}+k_{31} x_{11}+k_{10} x_{18}+k_{12}\\
f_{11}:=f_{c,\kappa}(x)_{11}=-k_{20} x_6 x_{11}-k_{24} x_{11} x_{12}+k_{30} x_{10}+(-k_{23}-k_{31}) x_{11}+k_{25} x_{13}+k_{21} x_{19}\\
f_{13}:=f_{c,\kappa}(x)_{13}=k_{24} x_{11} x_{12}-k_{25} x_{13}\\
f_{14}:=f_{c,\kappa}(x)_{14}=k_3 x_2 x_4+(-k_4-k_5) x_{14}\\
f_{15}:=f_{c,\kappa}(x)_{15}=k_{14} x_3 x_6+(-k_{15}-k_{16}) x_{15}\\
f_{16}:=f_{c,\kappa}(x)_{16}=k_6 x_5 x_8+(-k_7-k_8) x_{16}\\
f_{17}:=f_{c,\kappa}(x)_{17}=k_{17} x_7 x_9+(-k_{18}-k_{19}) x_{17}\\
f_{18}:=f_{c,\kappa}(x)_{18}=k_9 x_4 x_{10}+(-k_{10}-k_{11}) x_{18}\\
f_{19}:=f_{c,\kappa}(x)_{19}=k_{20} x_6 x_{11}+(-k_{21}-k_{22}) x_{19}\\
\end{cases}$%
\end{equation}

§.§ Bounds on the number of steady states

As is true for many polynomial systems arising from applications,
chemical reaction systems are structured and sparse.
Consequently, the study of chemical reaction networks benefits from the study of sparse polynomial systems.
We conclude this section by recalling a few fundamental upper bounds on the number of isolated solutions to sparse polynomial systems, and their applications to reaction networks.

a polynomial $f = a_1 x^{\sigma_1} + a_2 x^{\sigma_2} + \dots + a_{\ell} x^{\sigma_{\ell}} \in \mathbb{R}[x_1,x_2,\dots, x_d]$, where
exponent vectors $\sigma_i \in \mathbb{Z}^d$ are distinct and
$a_i \neq 0$ for all $i$,
the Newton polytope of $f$
is the convex hull of its exponent vectors:
${\rm Newt}(f) \coloneqq {\rm conv} \{\sigma_1, \sigma_2,\dots ,\sigma_{\ell}\} \subseteq \mathbb{R}^d.$
The (finite) set of exponent vectors $\{\sigma_1,\dots,\sigma_{\ell}\}$ is called the support of $f$.
The Kushnirenko Theorem gives
one classical bound on the maximum number of nonzero isolated complex solutions to a polynomial system where each polynomial has the same Newton polytope:

Consider $d$ real polynomials $g_1, g_2, \dots, g_d \in \mathbb{R}[x_1, x_2, \dots, x_d]$ with common support, and let $G:=\newt{g_i}$ for all $i$.
Then the number of isolated solutions in $(\C^*)^d$, counted with multiplicity, of the system
$g_1(x) = g_2(x) = \cdots = g_d(x) = 0$ is at most $d!~\vol{G}$.

In order to state a generalization of this result, we will need the following definition.

Let $Q_1, Q_2, \dots,Q_d\subseteq \R^d$ be polytopes.
The normalized volume of the Minkowski sum $\lambda_1Q_1+ \lambda_2 Q_2+ \dots+\lambda_d Q_d$ is a homogeneous polynomial of degree $d$ in nonnegative variables $\lambda_1,\lambda_2,\dots,\lambda_d$. In this polynomial, the coefficient
of $\lambda_1 \lambda_2\cdots\lambda_d$,
denoted by $\mv{Q_1, Q_2, \dots,Q_d}$, is the mixed volume of $Q_1,Q_2,...,Q_d$.

The Bernstein-Khovanskii-Kushnirenko theorem uses the mixed volume of the Newton polytopes of a polynomial system to count its solutions in $(\C^*)^d$.

Consider $d$ real polynomials $g_1, g_2, \dots, g_d \in \mathbb{R}[x_1, x_2, \dots, x_d] $.
Then the number of isolated solutions in $(\C^*)^d$, counted with multiplicity, of the system
$g_1(x) = g_2(x) = \cdots = g_d(x) = 0$ is at most $\mv{\newt{g_1}, \dots , \newt{g_d}}$.

The Kushnirenko Theorem and its generalization, the
Bernstein Theorem, connect convex geometry with algebraic geometry: the volumes of polytopes are used to count the number of solutions in $(\C^*)^d$ to a generic polynomial system
drawn from the
polynomial vector space spanned by the monomials with exponent vectors in the relevant polytopes.
Proofs of these bounds leverage rich theory from algebraic geometry.
Recently, Kaveh and Khovanskii <cit.> realized the Kushnirenko and Bernstein Theorems as special cases of their more general result, as we will see in Section <ref>.

Motivated by Proposition <ref>, <cit.> introduced the definition of the mixed volume of a chemical reaction network.

Let $G$ be a network with $d$ species,
$m$ reactions, and
a $p \times d$ conservation-law matrix $W$,
which results in the system augmented by conservation laws
$f_{c,\kappa}$, as in (<ref>).
Let $c^*\in \mathbb{R}^p_{\neq 0}$, and let
$\kappa^* \in \mathbb{R}^m_{>0}$ be generic.
Let $Q_1, Q_2, \dots,Q_d\subset \R^d$ be the Newton polytopes of $f_{c^*,\kappa^*,1},f_{c^*,\kappa^*,2}, \dots, f_{c^*,\kappa^*,d}$, respectively. The mixed volume of $G$ (with respect to $W$)
is the mixed volume of $Q_1,Q_2,\dots,Q_d$.

Every positive steady state is a steady state over $\mathbb{C}^*$.
Also, the mixed volume pertains to polynomial systems with the same supports
as the augmented system $f_{c,\kappa}=0$ (but without constraining the coefficients to come from a reaction network).  We obtain, therefore, the bounds in the following result:

For every reaction network, the following inequalities hold among the
maximum number of positive steady states,
the maximum number of steady states over $\mathbb{C}^*$, and
the mixed volume of the network (with respect to any conservation-law matrix):

\[
\mathrm{max~\#~of~positive~steady~states}
~\leq ~
\mathrm{max~\#~of~steady~states~over~} \mathbb{C}^*
~\leq ~
{\rm mixed~volume}~.
\]

Note that some networks have boundary steady states,
but boundary steady states are not counted by mixed volume.

As an example, we now compute the mixed volume of our motivating example (also computed in <cit.>).
(Other related mixed volumes were also computed in <cit.>.)

[Motivating Example, continued]
Recall the Wnt network from Example <ref>.
Let $e_i$ denote the $i$-th standard basis vector of $\R^{19}$.
The Newton polytope of $f_1$ is $\conv\{e_1,e_2,e_{14},e_{15}\}$, the Newton polytope of $f_2$ is $\conv\{e_2+e_4,e_1,e_2,e_3,e_{14}\}$, and
we can compute the remaining Newton polytopes $\newt{f_3},\dots,\newt{f_{19}}$ analogously.
By Definition <ref>, the mixed volume of $G$ is 56, and by Proposition <ref>, we conclude that $G$ has at most 56 positive steady states (for a generic choice of $\kappa,c$).
Supporting computations are archived in the repository <https://github.com/EliseAWalker/NOBodies-for-ChemicalRxns>

network     max number of positive steady states     mixed volume     max number of solutions over $\C^*$
fully irreversible ERK     1 <cit.>     3     3 (attained over $\R^*$ <cit.>)

full ERK     3 (conjecture) <cit.>     7     7 (attained over $\C^*$ <cit.>)

Edelstein    3     3    3 <cit.>

Wnt     3 (conjecture) <cit.>     56 (Example <ref>)     9 (attained over $\R$ <cit.>)

Selected results on mixed volume of chemical reaction networks and related bounds on the maximum number of positive steady states.

Table <ref> catalogs the maximum number of steady states, the mixed volume, and the computed maximum number of nonzero complex steady states for selected chemical reaction networks.
For some networks, like the fully irreversible ERK network <cit.> and the Edelstein network <cit.>, the mixed volume gives a good bound on the maximum number of observable positive steady states.
For others, like the Wnt network, the mixed volume greatly exceeds the maximum number of observable steady states, motivating the need for a much sharper bound.
using Gröbner-basis methods,
<cit.> showed that the Wnt model has a maximum number of 9 complex steady states.
Moreover, <cit.> numerically provided a witness of 3 positive steady states, and thus conjectured that the maximum number of positive steady states cannot exceed 3.
We conclude that the mixed-volume bound overcounts the actual number of observable steady states by at least $45$.
Can we do better?
In this paper, we introduce a complementary approach to those in <cit.>, adding a new algebraic-geometric tool for studying chemical reaction networks.
In the next section, we will show that for the Wnt network, the resulting bound is tighter than the mixed-volume bound.

§ NEWTON-OKOUNKOV BODY BOUND BACKGROUND

In <cit.>,
Kaveh and Khovanskii investigated a generalization of the mixed-volume bound. That is, instead of using the monomial-support system to give a bound on the number of solutions, they instead considered systems with polynomial support. The bound they gave is the birationally invariant intersection index. They furthermore gave a formula for the self-intersection index using Newton-Okounkov bodies. We will consequently refer to this bound as the Newton-Okounkov body bound. This bound is complicated to define and compute. This section concerns itself with the necessary definitions and concepts.

We begin with a brief overview of
the setup.
Let $X$ be a $d$-dimensional algebraic variety and let $V\subset \C(X)$ be a $n$-dimensional vector space with finite Khovanskii basis $\mathcal{B}=\{b_0,\dots,b_n\}$.
The rational Kodaira map $\phi_V:X\dashrightarrow\mathbb{P}(V')\simeq \mathbb{P}^n$, defined by $x\mapsto [v_0(x),\dots,v_n(x)]$ parametrizes a variety as follows.
The Zariski closure of the image of $\phi_V$ is a projective variety; call it $\kodproj:=\overline{\phi_V(X)}$.
The coordinate ring of $\kodproj$ is
a graded ring
Given a finite Khovanskii basis for $\rv$, one may compute a Newton-Okounkov body of $\rv$.
Convex bodies – Newton polytopes and Newton-Okounkov bodies – associated to projective varieties give us a beautiful way to compute
the self-intersection index of $V$.
Next, we define the objects needed to construct a Newton-Okounkov body associated to
a vector space $V$ of functions defined on a projective variety.

§.§ Polynomial vector space, graded algebra

A (nonempty, finite-dimensional) subspace $V\subset \C(X)$ gives rise to a graded algebra $\rv$ as follows.
Let $s$ be a formal grading variable, and let $V^k \subset \C(X)$ denote the subspace generated by all $k$-fold products of elements in $V$, for $k\in \N$ (here $V^0$ is just the base field $\C$).
Any nonzero element $f$ in the $\C$-algebra generated by $Vs$
can be decomposed into a sum of $s$-homogeneous components: $f=f_ks^k+\dots+f_1s+f_0$, where $f_k\neq 0$ and $f_j\in V^j$ for each $j$.
Decomposing $f\in\rv$ in this way, we see that $\rv \subset \C(X)[s]$, and $\rv$ is graded by the highest degree of $s$.
Let $\rv$$:=\bigoplus_{k\geq 0} V^ks^k$ denote this graded algebra associated to $V$.

[Motivating Example, continued]
In this example, we associate a finite-dimensional vector space $V$ to the Wnt network.
We also compute a finite Khovanskii basis of the corresponding algebra $\rv$.
Consider the following vector subspace $V$ of $\C(X)$ spanned by 29 generators $v_0,\dots,v_{28}$:

\begin{align*}
V = \spn_{\C}\{
&-k_3 x_2 x_4+ (-k_2-k_{26}) x_2+ (k_4+k_5) x_{14},~
k_1 x_1+ k_{27} x_3,~
-k_{14} x_3 x_6 -k_{27} x_3 +( k_{15}+ k_{16}) x_{15},\\
&k_{26} x_2,~
-k_6 x_5 x_8 -k_{28} x_5+k_7 x_{16},~
k_{29} x_7+k_5 x_{14},~
-k_{14} x_3 x_6 -k_{20} x_6 x_{11}+ k_{15} x_{15},~\\
&k_{19} x_{17}+( k_{21}+ k_{22}) x_{19},~
-k_{17} x_7 x_9-k_{29} x_7+ k_{18} x_{17},\\
&k_{28} x_5+ k_{16} x_{15},~
-k_9 x_4 x_{10}+(-k_{13}-k_{30}) x_{10}+ k_{10} x_{18},~
k_{31} x_{11}+ k_{12},\\
&-k_{24} x_{11} x_{12} + (-k_{23}-k_{31}) x_{11} + k_{25} x_{13},~
-k_{20} x_6 x_{11}+ k_{21} x_{19}+k_{30} x_{10},~
k_{24} x_{11} x_{12}-k_{25} x_{13},~\\
&k_3 x_2 x_4+(-k_4-k_5) x_{14},~
k_{14} x_3 x_6+(-k_{15}-k_{16}) x_{15},~
k_6 x_5 x_8+(-k_7-k_8) x_{16},\\
&k_{17} x_7 x_9+(-k_{18}-k_{19}) x_{17},~
k_9 x_4 x_{10}+(-k_{10}-k_{11}) x_{18},~
k_{20} x_6 x_{11}+(-k_{21}-k_{22}) x_{19},~\\
\}
\spn_{\C}\{v_0,\dots,v_{28}\}
\end{align*}

Note that $V$ is chosen such that the defining polynomials $f_i$ of the Wnt network's augmented system (<ref>) are $\C$-linear combinations of the basis polynomials $v_j$.
Let $\rv=\bigoplus_{k\geq 0} V^k s^k$ be the graded
algebra associated to $V$.

§.§ Degree of Kodaira map

Given a vector space $W$ over a field $F$ its dual space $W'$ is the $F$-vector space
of all linear functionals $W\rightarrow F$.
The projectivization of an $F$-vector space $W$ is $\pp(W)$ $:= W'/\sim$ where $\sim$ is the equivalence relation $u \sim v$ if $u=\lambda v$ for $u,v\in W'$ and $\lambda\in F^*$.

Let $V \subset \C(X)$ be nonzero finite dimensional vector subspace.
Suppose $\{v_0,\ldots,v_n\}$ is a vector space
basis for $V$.
Then for $x\in X$ such that $v_i(x)$ is defined for all $i$ and $v_i(x)\neq 0$ for some $i$,
the rational Kodaira map $\varphi_V:X\dashrightarrow \pp(V')$ of the subspace $V$ is given by:

\[\varphi_V(x)=[v_0(x):\dots:v_n(x)]~.\]

The degree $\text{deg}~\varphi_V$ of the Kodaira map $\varphi_V$
is the maximum
number of isolated points in the inverse image of any point:
The homogeneous coordinate ring of
the projective variety $\overline{\phi_V(X)}$
is isomorphic to the graded algebra $\rv$ <cit.>.

[Motivating Example, continued]
Recall that a vector space basis for $V$ given in Example <ref> is $\{v_0,\dots,v_{28}\}$.
The Kodaira map $\phi_V$ is defined by $x\mapsto [v_0(x): \dots: v_{28}(x)]$.
Using elimination techniques, we compute the degree $\deg(\phi_V)$ of the Kodaira map to be $1$.

§.§ Value semigroup

A valuation is a map that translates subalgebra computations into computations in an associated semigroup.
We give the following definition for a valuation on a $\C$-algebra, $\C(X)$.

Let $\succ$ be a total order on $\Z^d$ that restricts to a well-ordering on $\Z^d_{\geq 0}$.
A $\Z^d$-valuation on $\C(X)$ is a
surjective group homomorphism
$\nu:\C(X) \rightarrow \Z^d \cup \{\infty\}$
that satisfies the following properties for all $f,g\in \C(X)^*$ and $c\in\C^*$:

* $\nu(0)=\infty$, where $\infty \succ \alpha$, and $\alpha+\infty=\infty$ for all $\alpha \in \Z^d$,

* $\nu(c)=(0,\dots,0)$, and

* $\nu(f+g)\succeq \text{min}\{\nu(f),\nu(g)\}$.

Throughout this paper, we will restrict ourselves to a specific type of valuation induced by a monomial term order called a Gröbner valuation.

Fix a total ordering on $\Z^d$.
Define a Gröbner valuation $\nu$ on $\C(X)$ as follows.
For $f\in \C[X]$, let $cx_1^{a_1}\dots x_d^{a_d}$ be the term in $f$ with the smallest exponent $(a_1,\dots,a_d)$ with respect to the total ordering on $\Z^d$.
Define $\nu(f)=(a_1,\dots,a_d)$.
Extend $\nu$ to $\C(X)$ by defining $\nu(f/g)=\nu(f)-\nu(g)$ for any $f,g\in \C[X]$ with $g\neq 0$.
A Gröbner valuation $\nu:\C(X)\rightarrow \Z^d$ is a
$\Z^d$-valuation with one-dimensional leaves on $\C(X)$ <cit.>.

The graded reverse lexicographic term order (grevlex)
on the set of monomials in a (multivariate) polynomial ring $k[{\bf x}]$ is defined by: $x^A > x^B$ if either degree($x^A$) $>$ degree($x^B$) or degree($x^A$) = degree($x^B$) and the last nonzero entry of the vector of integers $A-B$ is negative.

Every monomial term order, including grvlex, is extendable to a Gröbner valuation where the smallest exponent in the valuation (in the total order $\succ$ on $\Z^d$) corresponds to the exponent of the leading term of the monomial term order.

We extend a $\Z^d$-valuation $\nu$ on $\C(X)$ to a $\Z^d\oplus \N$-valuation
$\widehat{\nu}:\rv\rightarrow \Z^d\oplus \N$ on the graded algebra $\rv$ as follows:
for $f\in \rv^*$ with degree $k$ as a polynomial in $\C(X)[s]$,

\[\widehat{\nu}(f) = \widehat{\nu}(f_ks^k+\dots+f_1s+f_0) = (\nu(f_k), k)~.\]

The image of $\rv^*$ under $\widehat{\nu}$ is defined to be the value semigroup $S$
of $(V,\nu)$:
that is, $S:= S(V,\nu) = \{\widehat{\nu}(f)~|~f\in \rv^*\}$.

We remark that although $S$ is more precisely a monoid ($0\in S$ since $\nu(c)=0$ for any $c\in \C^*$), we follow the convention to refer to $S$ as a semigroup.
The rank of a valuation $\widehat{\nu}$ is the rank of the sublattice of $\Z^d\oplus \N$ generated by the value semigroup $S(V,\nu)$.
In all cases we consider, the valuation has full rank, that is, $\text{rank}(\widehat{\nu})=d+1=\dim(S)$.

§.§ Khovanskii basis

A Khovanskii basis is a special generating set of an algebra, first formalized in <cit.>.

Let $\nu$ be a valuation.
subset $\mathcal{B}\subset \rv$ is a Khovanskii basis for $(V,\nu)$
if the set $\{\nu(b)~|~b\in\mathcal{B}\}$ generates $S(V,\nu)$ as a semigroup.

It follows from the definition that a Khovanskii basis $\mathcal{B}$ generates $\rv$ as a $\C$-algebra, and $\mathcal{B}\cap Vs$ is a vector space basis for $Vs$.

We will only work with
graded polynomial
algebras $\rv$
that have a finite Khovanskii
basis with respect to some Gröbner valuation induced by a monomial term order.
Under this assumption (namely, a finite Khovanskii basis exists),
<cit.> gives an algorithm for computing a finite Khovanskii basis $\mathcal{B}=\{b_1,\dots,b_n\}$ from a finite set of $k$-algebra generators for $\rv$.
The recent Macaulay2 package SubalgebraBases.m2 <cit.> implements <cit.>'s algorithm;
we use SubalgebraBases.m2 to compute Khovanskii bases in our examples (Section <ref>).

[Motivating Example, continued]
Let $\nu$ be the Gröbner valuation on $\C(X)$ induced by the graded reverse lexicographic (grevlex) monomial term order on $\C[x_1,\dots,x_{19}]$.
We extend the Gröbner valuation on $\C(X)$ to $\C(X)[s]$, so that $\widehat{\nu}(f) = (\nu(f), \deg_s(f))$.
The vector space basis $v_0,\dots,v_{28}$ induces a vector space basis $v_0s,\dots,v_{28}s$ for $Vs$.
Using the vector space basis $v_0s,\ldots,v_{28}s$ for $Vs$ as input to SubalgebraBases.m2, we compute the following Khovanskii basis $\mathcal{B}$ (with respect to $\nu$) for $\rv$:

\begin{align}\label{eq:wnt-khovanskii-basis}
\Big\{    &s,\,{x}_{19}s,\,{x}_{18}s,\,{x}_{17}s,\,{x}_{16}s,\,{x}_{15}s,\,{x}_{14}s,\,{x}_{12}s+{x}_{13}s,\,{x}_{11}s,\,{x}_{10}s,\,{x}_{9}s,\,{x}_{8}s,\,{x}_{7}s,\,{x}_{5}s,\,{x}_{4}s+{x}_{6}s,\,{x}_{3}s,\,
\notag \\
\,{x}_{6}{x}_{11}s,\,{x}_{4}{x}_{10}s,\,{x}_{7}{x}_{9}s,\,{x}_{5}{x}_{8}s,\,{x}_{3}{x}_{6}s,\,{x}_{2}{x}_{4}s,\,{x}_{11}{x}_{13}s^{2}+\frac{{k}_{25}}{{k}_{24}}\,{x}_{13}s^{2},\,
\notag \\
\notag \\
\begin{split}
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad +\frac{2\,{k}_{25}}{{k}_{24}}\,{x}_{2}{x}_{6}^{2}{x}_{11}{x}_{13}^{2}s^{5}+\frac{{k}_{25}^{2}}{{k}_{24}^{2}}\,{x}_{2}{x}_{4}{x}_{6}{x}_{13}^{2}s^{5}+\frac{{k}_{25}^{2}}{{k}_{24}^{2}}\,{x}_{2}{x}_{6}^{2}{x}_{13}^{2}s^{5},\,\end{split}\notag \\
\begin{split}  &{x}_{4}{x}_{6}{x}_{10}{x}_{11}^{2}{x}_{13}^{2}s^{5}+{x}_{6}^{2}{x}_{10}{x}_{11}^{2}{x}_{13}^{2}s^{5}+\frac{2\,{k}_{25}}{{k}_{24}}\,{x}_{4}{x}_{6}{x}_{10}{x}_{11}{x}_{13}^{2}s^{5}\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad +\frac{2\,{k}_{25}}{{k}_{24}}\,{x}_{6}^{2}{x}_{10}{x}_{11}{x}_{13}^{2}s^{5}+\frac{{k}_{25}^{2}}{{k}_{24}^{2}}\,{x}_{4}{x}_{6}{x}_{10}{x}_{13}^{2}s^{5}+\frac{{k}_{25}^{2}}{{k}_{24}^{2}}\,{x}_{6}^{2}{x}_{10}{x}_{13}^{2}s^{5} \Big\}\end{split}
\end{align}

The Khovanskii basis has 32 generators, which we denote $b_0,\dots,b_{31}$.
For each polynomial, the leading monomial (with respect to grevlex) of each basis vector is the first monomial listed.
By definition, the semigroup $S(V,\nu)$ is generated by the valuations $\widehat{\nu}(b_i)$ of the Khovanskii basis elements.
Let $e_i$ denote the $i$-th standard basis vector for $\R^{19}$.
Then, for example,
})\,{x}_{13}s\right)=(e_{11}+e_{12},1)\in S(V,\nu)$
since the Khovanskii basis element $b_{19} = {x}_{11}{x}_{12}s-({k}_{25}/{k}_{24
})\,{x}_{13}s$ has leading term ${x}_{11}{x}_{12}s$ under grevlex.
The Gröbner valuation of the remaining $b_i$ are computed analogously.
The valuations of the 32 Khovanskii basis elements are the columns of the following matrix $A$:

\[
\resizebox{.8\textwidth}{!}{
\]

By definition, since $\mathcal{B}$ is a Khovanskii basis for $(V,\nu)$, the semigroup $S(V,\nu)\subset \Z^{19}\bigoplus \N$ is generated by the columns of the matrix $A$.

Deciding whether a finitely generated subalgebra has a finite Khovanskii basis is still an open question, in general.
We refer the interested reader to <cit.> for examples of
intricacies of computing Khovanskii bases.

§.§ Index of graded algebra

In order to define the index associated to $\rv$, we will need a couple more definitions.
Let $G(S)\subset \R^{d+1}$ denote the group generated by the semigroup $S$, so
$G(S)~=~\left\{\sum_i k_i a_i~|~a_i\in S, k_i\in \Z\right\}$.
Let $\pi:\R^{d+1}\rightarrow \R$ denote projection onto the $(d+1)$-st coordinate.
Let $G_0(S)\subset G(S)$ denote the subgroup $\pi^{-1}(0)\cap G(S)$ of $\Z^d\oplus \{0\}$.
The index of $\rv$, $\text{ind}(\rv)$,
is the index of the subgroup
$G_0(S)$ in $\Z^d \bigoplus \{0\}$.

Consider $\rv$ where $V=\spn\{x_1(x_1^2+x_2^2 -2x_1),~x_1(5-4x_2),~x_2(x_1^2+x_2^2 -2x_1),~x_2(5-4x_2)\}$. Using $s$ as a grading variable, this is equivalent to considering the algebra generated by $\{x_1(x_1^2+x_2^2 -2x_1)s,~x_1(5-4x_2)s,~x_2(x_1^2+x_2^2 -2x_1)s,~x_2(5-4x_2)s\}$ over $\C$.

Using the Grobner valuation induced by grevlex
as the valuation $\widehat{\nu}$,
one finds that the the algebra generators are also a finite Khovanskii basis for $\rv$ with respect to this valuation.
Thus, the value semigroup $S$ is generated by the points: $(0,2,1), (1,1,1),(2,1,1), (3,0,1)$.
The group generated by $S$ is $\Z^2 \oplus \Z$.
By taking integer combinations of the generators  (specifically, pairwise sums of $\pm(0,2,1), \pm(1,1,1), \pm(2,1,1), \pm(3,0,1)$), then one finds that $G_0(R)$ contains the points $(1,-1,0), (-1,0,0), (1,0,0), \text{ and } (-1,1,0)$, which is enough to generate all of $\Z^2 \oplus \{0\}$.
Thus $G_0(R) = \Z^2 \oplus \{0\}$, and so $\ind(R) = [\Z^2 \oplus \{0\}: G_0(R)] = 1$.

[Motivating Example, continued]
Recall that the semigroup $S(V,\nu)$ is generated by the images $\widehat{\nu}(b_i)$ of the Khovanskii basis elements $b_0,\dots,b_{31}\in\mathcal{B}$, computed in Example <ref>.
Also, $G(S)$ is precisely the group generated by these $\widehat{\nu}(b_i)$.
Here, $G(S)=\Z^{19}\bigoplus \Z$, meaning that the valuations generate the entire integer lattice as a group.
Consequently, $G_0(S)=\Z^{19}\bigoplus \{0\}$.
From this, we conclude that $\text{ind}(\rv)=[\Z^d \times \{0\}:G_0(S)]=1$.

§.§ Newton-Okounkov body

The Newton-Okounkov cone of $(V,\nu)$
is defined to be the closure of the convex hull of the value semigroup $S$
in $\R^{d+1}$.

Finally, we define our main computational object: a Newton-Okounkov body associated to $(V, \nu)$.

Let $(V,\nu)$ be a subspace equipped with a valuation $\nu$ as defined in Section <ref>. The Newton-Okounkov body $\Delta(V,\nu)$ associated to $(V,\nu)$ is defined to be the base of the Newton-Okounkov cone $\text{cone}(V,\nu)$, that is,
$\Delta(V,\nu):=\text{cone}(V,\nu)\cap (\R^d\times \{1\}).$

\[\Delta(V,\nu) = \overline{\conv\left(\bigcup\limits_{k>0}^{} \{\nu(f)/k~|f\in V^k s^k\}\right)}\subset \R^d~.\]

We conclude this section by returning to our motivating example and illustrating the definitions and notation introduced here.

[Motivating Example, continued]
In Example <ref>, we determined a Khovanskii basis (<ref>) associated to the Wnt network, and from this basis we constructed the associated semigroup $S(V,\nu)$.
The cone $\cone(V)$ is the full-dimensional cone in $\R^{20}$ spanned by the columns $a^{(1)},\ldots,a^{(32)}\in \Z^{20}$ of $A$.
Writing the column as $a^{(i)}$ as $\left(a_1^{(i)}, a_2^{(i)}, \dots, a_{19}^{(i)}, a_{20}^{(i)}\right)$, the generators of the Newton-Okounkov body $\nobv$ are obtained as follows:
divide each of the first $19$ components
by the last component
What remains is $\widetilde{a^{(i)}} := \left(a_{1}^{(i)}/a_{20}^{(i)},~a_2^{(i)}/a_{20}^{(i)},\dots,a_{19}^{(i)}/a_{20}^{(i)}\right) \in \R^{19}$, for $i=1,\dots,32$.
The Newton-Okounkov body for Wnt is then the convex hull of the $\widetilde{a^{(i)}}$ for $i=1,\dots,32$.
It is a $19$-dimensional polytope sitting in $\R^{19}$.

In Proposition <ref>, we will see how to use the Newton-Okounkov body associated to $(V,\nu)$ to bound the number of isolated solutions to a general polynomial system in $V$.
Kaveh and Khovanskii <cit.> showed that, analogous to the setting of a toric variety and its corresponding Newton polytope, asymptotic behavior of the Hilbert function of a graded algebra associated to an arbitrary projective variety can be understood from combinatorial data of its Newton-Okounkov body.
As we will see in the next section, a consequence of this is that
the volume of the Newton-Okounkov body associated to $(V,\nu)$ is an invariant that bounds the number of effective solutions to a general polynomial system drawn from $V$.

§.§ The self-intersection index: a bound on isolated solutions

Kushnirenko's theorem implies that the number of solutions (in ${(\C^*)}^d$) to a general polynomial system with support in a finite set $A\subset \Z^d$ is an invariant of the polynomial vector space $V$ generated by the monomials with exponents in $A$.
In other words,
the number of solutions to a general system of polynomials drawn from $V$ depends on the subspace $V$, and has less on the particular polynomial system drawn from it.
This generalizes this to more general vector spaces $V$: <cit.> showed that
the number of solutions (in some distinguished subset of an algebraic variety) of a system of rational functions in a rational space is an invariant of the subspace $V$, and is called the self-intersection index of V (Definition <ref>).
Moreover, the actual number of solutions is computed using the constructions from the previous section (Proposition <ref>).

We carefully explain the results from <cit.> now.
As in the previous section,
$X$ is an irreducible $d$-dimensional complex algebraic variety with function field $\mathbb{C}(X)$,
and $V$ is a nonzero finite-dimensional subspace of rational functions on $X$.

Let $(V_1,\dots,V_d)$ be a $d$-tuple of nonzero finite-dimensional subspaces of rational functions on $X$.
Let ${\bf V}:= V_1 \times \dots \times V_d$.
Take $\uv \subset X$ to be the set of all nonsingular points in $X$ at which every function in each $V_i$ is regular.
Let $\zv\subset \uv$ be the collection of all common zeros from each $V_i$:

\[\zv :=\bigcup_{i=1}^d\{x\in U_{\bf V}~|~f(x)=0~\text{~for~all~}f\in V_i\}.~\]

We refer to $\uv$ as the smooth locus
of $X$ with respect to
${\bf V}$ and to $\zv$ as the base locus of ${\bf V}$.

We will only consider vector spaces $V$ spanned by a finite set of polynomial functions.
Consequently, every function in $V$ is regular, and $\uv$ is just the set of nonsingular points of $X$.
We will be able to count effective solutions to general polynomial systems.

Let $(f_1,\dots,f_d)\in V_1\times\dots\times V_d$ be a $d$-tuple of rational functions.
An effective solution of $f_1=\dots=f_d=0$ is a nondegenerate $x$ in the set $\uv\setminus\zv$ such that $f_1(x)=\dots=f_d(x)=0$.

Kaveh and Khovanskii <cit.> showed that any $d$-tuple $(V_1,\dots,V_d)$ of rational functions has a proper algebraic subvariety $\textbf{R} \subset \textbf{V}:=V_1\times \dots \times V_d$ such that for any  $(f_1,\dots,f_d)\in \textbf{V}\setminus\textbf{R}$:
the number of solutions to the system $f_1=\dots=f_d=0$ in the set $U_{\textbf{V}}\setminus Z_{\textbf{V}}$ is independent of the choice of general $(f_1,\dots,f_d)$,
each solution $a\in U_{\textbf{V}}\setminus Z_{\textbf{V}}$ of the system $f_1=\dots=f_d=0$ is nondegenerate.
We call such a $d$-tuple $(f_1,\dots,f_d)$ of rational functions a general system drawn from $V$.

Let $(V_1,\dots,V_d)$, where $V_i\subset \C(X)$, be a $d$-tuple of subspaces of rational functions on $X$.
Let $(f_1,\dots,f_d)\in V_1\times \dots \times V_d$ be a general system.
The (birationally-invariant) intersection index $[V_1,\dots,V_d]$ is defined to be the number of
effective solutions to  $f_1=\dots=f_d=0$.

In this paper, we will only consider the case where the $V_i$ are all equal, so the $d$-tuple from Definition <ref> is $(V,\dots,V)$ for some subspace $V$ of rational functions on $X$.
In this case, $[V,\dots,V]$ is called the (self-)intersection index of $V$.
We emphasize that the intersection index of $V$ counts effective solutions to a system $f_1=\dots=f_d=0$,
where each $f_i$ is a general element
drawn from $V$.

We now introduce the key result central to our application: the self-intersection index of a subspace of rational functions is proportional to the volume of its associated Newton-Okounkov body.

Let $X$ be a $d$-dimensional irreducible complex algebraic variety and
let $V$ be a nonzero finite-dimensional vector subspace of rational functions on $X$ with Kodaira map $\varphi_V$.
If $\varphi_V$ has finite mapping degree, then

\[[V,\dots,V]=\dfrac{d!~\deg~\varphi_V}{\ind~\rv} \vold{\Delta(\rv)}~.\]

Proposition <ref> says that if $\{f_1=\dots=f_d=0\}$ is a system of polynomial equations, where each $f_i$ is general in $V$,
then the number of effective solutions to $\{f_1=\dots=f_d=0\}$
can be computed explicitly in terms of the algebraic-geometric objects introduced in this section.

[Motivating Example, continued]
In Example <ref>,
we constructed the 19-dimensional Newton-Okounkov body of the Wnt network, a 19-dimensional polytope that is the convex hull of the $\widetilde{a^{(1)}}, \dots, \widetilde{a^{(32)}}$.
So, $d=19$, and also
the normalized volume $19!\vnob$ of this convex body is 32.
In Example <ref>, we computed $\text{ind}(\rv)=1$.
In Example <ref>, we showed that $\deg{\phi_V}=1$.
So, by Proposition <ref>, the self-intersection index of $V$ is 32, from which we conclude  the number of effective solutions to the polynomial system (<ref>) is at most 32.

As described in this section, a Newton-Okounkov body is associated to a
polynomial vector space $V$ and a valuation $\nu$ (more specifically, to a finite Khovanskii basis $\mathcal{B}$ for $(V,\nu)$) via a graded ring $\rv$.
Through our motivating example, we showed how to apply the setup to a chemical reaction network to obtain a Newton-Okounkov body associated to a chemical reaction network, whose volume gives a bound on the maximum number of steady states of the network.
In the next section, we state this process more explicitly
(Procedure <ref>).
While the motivation in this paper is polynomial systems arising from chemical reaction networks, we emphasize that the procedure applies more generally to any sparse polynomial system one wishes to solve.

§ NEWTON-OKOUNKOV BODIES OF CHEMICAL REACTION NETWORKS

Given a chemical reaction network, this section describes how to compute an associated Newton-Okounkov body bound. This first involves finding a vector space of functions that includes the chemical reaction system. After
choosing such a
vector space $V$, one computes the self-intersection of $V$.
The corresponding Newton-Okounkov body bound for the associated network is defined as this self-intersection index.

§.§ Procedure to compute the self-intersection index of a chemical reaction network

We give a procedure for computing the self-intersection index of a chemical reaction system defined by polynomials $\{f_1,\ldots,f_d\}$, where $f_i\in \C[x_1,\dots,x_d]$.
Specifically, we construct
a vector space $V$ (Section <ref>),
the Kodaira map $\varphi_V$ (Section <ref>),
a valuation $\nu$ (Section <ref>),
a Khovanskii basis $\mathcal{B}$ (Section <ref>),
the index of $\rv$ (Section <ref>),
the Newton-Okounkov body $\nobv$ (Section <ref>),
and finally apply Proposition <ref>.

Input: A chemical reaction system defined by polynomials $\{f_1,\ldots,f_d\}$, where $f_i\in \C[x_1,\dots,x_d]$.

Output: The self-intersection index for the chemical reaction system, i.e., an upper bound on the number of its positive steady states for any generic choice of rate constants.

*  Fix a monomial term order $<$.

*  Choose a generating set $\{v_1,\dots,v_n\}$ for a polynomial vector space $V$ such that each $f_i$ is
in the complex span of the $v_j$, i.e. $f_i\in V:=\spn_{\C}\{v_1,\ldots,v_n\}$ for $i = 1, \dots, d$.

*  Compute $\deg(\varphi_V)$ according to Section <ref>.

*  Compute
a finite
Khovanskii basis $\mathcal{B}$ for the subalgebra $\rv$ associated to $V=\spn_{\C}\{v_1,\ldots,v_n\}$ according to Section <ref>.
If a finite Khovanskii basis does not exist, return to Step 0, and choose a different monomial term order $<$.

*  Compute $\ind(\rv)$ according to Section <ref>.

*  Compute the volume $\vnob$ of the Newton-Okounkov body $\nobv:=\Delta(V,<,\mathcal{B})$
according to Section <ref>.

*  Output the birationally invariant self-intersection index $[V,\dots,V]$, according to Proposition <ref>.

§.§ Practical Considerations

One strategy for choosing a generating set for input to SubalgebraBases.m2 is to look for commonalities among the defining polynomials that we wish to solve.
For instance, when the defining polynomials of the system admit common factors (c.f., a product decomposition as in <cit.>) or common summands, we choose these shared polynomials in the list of vector space generators.
We apply this strategy in Example <ref> and Example <ref>.
Step <ref> is the most challenging step of Procedure <ref> and requires some trial-and-error by the user.

The bound $[V,\dots,V]$ computes all effective solutions (Definition <ref>) to any general system $g_1=0,\dots,g_d=0$, where $g_i$ is general in $V$.
Recall that this means that all solutions counted are outside $Z_{\textbf{V}}$, the base locus of $V$ as in Definition <ref>.
Note in particular that this means there could be additional real positive steady states in the base locus.
For each example we consider, the base locus is easily computed
and does not contain any positive real solutions.

For our examples, in Step 0 of Procedure <ref>, we will always select the monomial term order grevlex first.
The leading term under grevlex yields the smallest exponent under the Gröbner valuation (Definition <ref>).

A finite Khovanskii basis is not guaranteed.
However, if the semigroup $S(V,\nu)$ is finitely generated, then there is a finite Khovanskii basis for $(V,\nu)$ <cit.>.
If a finite Khovanskii basis is not achieved for a chosen valuation, one strategy is to choose a different basis
as input to  SubalgebraBases.m2.
Another strategy is to try a different valuation, like the Gröbner valuation induced by a different monomial term order.
For Example <ref>, the same vector space for $V$ basis together with the Gröbner valuation corresponding to the lex monomial order instead produced an infinite sagbi basis.

For small examples, the degree of the Kodaira map can be computed by hand.
We verify the computation using the Macaulay2 package Cremona.m2 <cit.>.
Note that Cremona.m2 requires a homogeneous map as input: in our computations, we multiply each $v_i$ by a power of a homogenizing variable.

§ EXAMPLES

In this section, we explicitly compute Newton-Okounkov bodies of select chemical reaction networks.
We also calculate the mixed volume bound for each reaction network, as in <cit.>.
The supporting computations for each example are archived at: <https://github.com/EliseAWalker/NOBodies-for-ChemicalRxns>.
Recall (Remark <ref>) that for each example, we choose the Gröbner valuation induced by the grevlex monomial order when applying Step <ref>.

[Wnt, Motivating Example <ref>, concluded]
As shown in Example <ref>, the Newton-Okounkov body bound for the Wnt network is 32, which is tighter than the mixed-volume bound of 56.
Furthermore, there are no positive real solutions in the base locus (indeed, the base locus is empty), and hence we expect no more than 32 positive steady states for Wnt, for generic parameters.

Newton-Okounkov bodies for Examples <ref>, <ref>, and <ref>, respectively.

[KST <cit.>]
Consider the following system:

\begin{equation}\label{eq:poly-sys-2circ-ex-non-gen-coeff}
\begin{aligned}
f_1 &= k_1(x^2+y^2) -k_2 x + k_3\\
f_2 &= k_1(x^2+y^2) -k_4 y + k_5 ~.
\end{aligned}
\end{equation}

where each $k_i\in\R_{> 0}$.
This polynomial system can be
obtained from a chemical reaction network with $f_{\ell}=\dot{x}_{\ell}$ since each defining polynomial is of the shape $f_{\ell} = p_{\ell} - x_{\ell} q_{\ell}$ for real polynomials $p_{\ell},~q_{\ell}$ with nonnegative coefficients <cit.>.
We note here that some information is known about the
in particular, the coefficients of $(x^2+y^2)$ in
$f_1$ and $f_2$
are equal.

In this small example, we can explicitly compute the maximum number of complex solutions.
Geometrically, the system defines the intersection of two circles, so the system has at most 2 isolated complex solutions.
Moreover, we produce a witness for two positive, real solutions by taking $k_1=1,~k_2=8,~k_3=7,~k_4=6,~k_5=3$.
In other words, the maximum number of steady states for the associated chemical reaction system is 2.
Next we apply Procedure <ref>, and show that the associated self-intersection index gives a tight bound.

Input: The chemical reaction system $\{f_1,f_2\}\subset\C[x,y]$.

Step <ref>: Let $V=\spn_{\C}\{1,x,y,(x^2+y^2)\}$.
Here, since the support of $f_1$ and $f_2$ share the summand $(x^2+y^2)$, we choose the polynomial instead of the individual monomials $x^2$ and $y^2$ as generators of the vector space.

Step <ref>: The Kodaira map of $V$ is defined by $(x,y)\mapsto [s:sz:ys:(x^2+y^2)s]$, and it has degree 1.

Step <ref>: Let $Vs$ be generated by $\{s,xs,ys,(x^2+y^2)s\}$.
We verify that the generating set at height 1 is in fact a Khovanskii basis
for the graded algebra $\rv$ associated to $V$: so, $\mathcal{B}=\{s,xs,ys,(x^2+y^2)s\}$.

Step <ref>: The valuations of the elements of $\mathcal{B}$ are $(0,0,1), (1,0,1), (0,1,1),$ and $(2,0,1)$.
They generate the value semigroup $S(V,\nu)$.
The subgroup $G_0(\rv) = \Z^2\oplus \{0\}$.
The index of the lattice $\ind~\rv$, that is, the index of $G_0(\rv)$ in $\Z^2\oplus \{0\}$, is 1.

Step <ref>: Then $\nobv$ is the triangle $\conv\{(0,0),(0,1),(2,0)\}$ in $\R^2$.

Step <ref>: The self-intersection index $[V,V]$ is $2!\vnob = 2$.

Output: This chemical reaction system has a maximum of 2 positive steady states (outside the base locus).
Note that for this example, the base locus, which consists of all solutions to $\{1=0, x=0, y=0, (x^2+y^2)=0\}$, is empty.

A general system drawn from $V$ has the form

\[
\begin{cases}
\kappa_1(x^2+y^2) + \kappa_2 y + \kappa_3 x + \kappa_4 = 0\\
\kappa_5(x^2+y^2) + \kappa_6 y + \kappa_7 x + \kappa_8 = 0~.
\end{cases}
\]

with $\kappa_i\in \C$. Thus, the mixed volume of a general system in $V$ is 4.
However, any such generic system has exactly 2 solutions, since, the system defines the intersection of two circles.
In other words, the maximum number of complex-number steady states is equal to the Newton-Okounkov body bound.
This Newton-Okounkov body bound is also tighter than the mixed volume bound of 4 for the original system (<ref>).

[PM <cit.>]
The following system of polynomials comes from a chemical reaction network and was studied as a network with absolute concentration robustness in <cit.>:

\begin{equation}\label{eq:pm}
\begin{aligned}
f_1(x) &= x_1[(x_1-1)^2 + (x_2-2)^2] \\
f_2(x) &= x_2[(x_1-1)^2 + (x_2-2)^2]
\end{aligned}
\end{equation}

It is straightforward to see that this system has exactly one positive steady state $(x_1^*,x_2^*)=(1,2)$.
In what follows, we will ignore the fact that this system is simple to solve by inspection and highlight how to apply our tools from Section <ref> to this toy example.
In particular, we
emphasize that the degree of the Kodaira map is not always 1 and needs to be computed explicitly in order to apply Proposition <ref>.

Input: the chemical reaction system $\{f_1,f_2\}\subset\C[x_1,x_2]$.

Step <ref>: This system belongs to the vector space $V = \spn_{\C}\{x_1(x_1^2+x_2^2 -2x_1),~x_1(5-4x_2),~x_2(x_1^2+x_2^2 -2x_1),~x_2(5-4x_2)\}$.
In particular, $f_1$ is equal to the sum of the first two basis elements of $V$, and $f_2$ is the sum of the last two basis elements.

Step <ref>: The Kodaira map of $V$ is defined by $(x,y)\mapsto [x_1(x_1^2+x_2^2 -2x_1)s: x_1(5-4x_2)s: x_2(x_1^2+x_2^2 -2x_1)s: x_2(5-4x_2)s]$, and it has degree 2.

Step <ref>: The generating set induces the Khovanskii basis $\mathcal{B}=\{x_1(x_1^2+x_2^2 -2x_1)s,~x_1(5-4x_2)s,~x_2(x_1^2+x_2^2 -2x_1)s,~x_2(5-4x_2)s\}$
for the graded algebra $\rv$ associated to $V$.

Step <ref>:
In Example <ref>, we computed that $\ind(\rv)=[\Z^2\oplus\{0\}:G_0(\rv)]=1$.

Step <ref>: The Newton-Okounkov body $\nobv$ has vertices $(3,0),(1,1),(2,1),(0,2)$ and volume 1.

Step <ref>: The self-intersection index $[V,V]$ is $d!\cdot \deg~\varphi_L\cdot\vnob /~\ind~\rv= 2!\cdot 2\cdot 1 / 1= 4$.

Output: The maximum number of positive steady states (outside the base locus) is 4.

Moreover, for this example, no positive solutions are missed by the self-intersection bound.
Indeed, the base locus of $V$ contains three points: $Z_{{\bf V}}= \{ (1\pm.75i, 1.25), (0,0)\}$.
While four is an over-calculation for the single known positive steady-state, the steady-state $(1,2)$ actually has multiplicity four, and so the calculation for the number of roots (counted with multiplicity) is tight.

[BCY <cit.>]
Networks with infinitely many positive steady states in a compatibility class were investigated in <cit.>.
These networks were the first such “weakly reversible” networks discovered, and <cit.> showed that these networks
have infinitely many positive steady states defined by a curve of steady states arising from a common factor of the individual equations in the chemical reaction system.
However, in general, the number of isolated (complex) solutions to these networks is small.

[] (0) at (-4,0) $\ce{0}$;
[] (x) at (-1.5,0) $\ce{X}$;

[] (y) at (-4,2) $\ce{Y}$;
[] (xy) at (-1.5,2) $\ce{X + Y}$;

[->,> = latex]  (0) – node[anchor=north] $1$ (x);
[->,> = latex]  (xy) – node[anchor=south] $5$ (y);

[->,> = latex]  (xy) – node[anchor=east] $1$ (x);

[->,> = latex]  (y) – node[anchor=east] $1$ (0);

[] (2x) at (5,0) $\ce{2X}$;
[] (3x) at (7.5,0) $\ce{3X}$;

[] (2xy) at (5,2) $\ce{2X + Y}$;
[] (3xy) at (7.5,2) $\ce{3X + Y}$;

[->,> = latex]  (2x) – node[anchor=north] $1$ (3x);
[->,> = latex]  (3xy) – node[anchor=south] $1$ (2xy);

[->,> = latex]  (3x) – node[anchor=east] $1$ (3xy);

[->,> = latex]  (2xy) – node[anchor=east] $5$ (2x);

[] (2y) at (0.25,0) $\ce{2Y}$;
[] (x2y) at (2.75,0) $\ce{X + 2Y}$;

[] (3y) at (0.25,2) $\ce{3Y}$;
[] (x3y) at (2.75,2) $\ce{X + 3Y}$;

[->,> = latex]  (3y) – node[anchor=east] $1$ (2y);
[->,> = latex]  (x3y) – node[anchor=south] $1$ (3y);

[->,> = latex]  (x2y) – node[anchor=east] $5$ (x3y);

[->,> = latex]  (2y) – node[anchor=north] $1$ (x2y);

[] (2x2y) at (10,0) $\ce{2X + 2Y}$;
[] (3x2y) at (13,0) $\ce{3X + 2Y}$;

[] (2x3y) at (10,2) $\ce{2X + 3Y}$;
[] (3x3y) at (13,2) $\ce{3X + 3Y}$;

[->,> = latex]  (2x2y) – node[anchor=north] $5$ (3x2y);
[->,> = latex]  (3x3y) – node[anchor=south] $1$ (2x3y);

[->,> = latex]  (3x2y) – node[anchor=east] $1$ (3x3y);

[->,> = latex]  (2x3y) – node[anchor=east] $1$ (2x2y);

A chemical reaction network with infinitely many positive steady states, but only one isolated positive steady state.

Consider the chemical reaction network in Figure <ref> from <cit.>.
The corresponding chemical reaction system is given by:

\[
\begin{cases}
\dot{x} = (x^2 y^2 + x^2 + y^2 + 1 - 5xy)[1 - xy]\\
\dot{y} = (x^2 y^2 + x^2 + y^2 + 1 - 5xy)[x - y]~.
\end{cases}
\]

The set of positive steady states consists of the non-isolated points $(x,y)$ on the curve $x^2 y^2 + x^2 + y^2 + 1 - 5xy = 0$, and the only isolated solution that is positive is the steady state $(1,1)$ <cit.>.
Notice that there is one more complex (non-positive) solution to $\dot{x}=\dot{y}=0$, namely $(-1,-1)$.

Step <ref>: Let $p_1:=x^2 y^2 + x^2 + y^2 + 1 - 5xy$ (the common factor of the system's two defining polynomials), and let $p_2:=1$, $p_3:=xy$, $p_4:=x$, $p_5:=y$ (these are the monomials in the remaining factors).
Notice that $\dot{x} = p_1(p_2 - p_3)$ and $\dot{y} = p_1(p_4 - p_5)$.
Take $V=\spn_{\C}\{p_1p_2,~p_1p_3,~p_1p_4,~p_1p_5\}$, so $\dot{x}, \dot{y} \in V$.

Step <ref>: The Kodaira map of $V$, defined by $(x,y)\mapsto [p_1p_2s:~p_1p_3s:~p_1p_4s:~p_1p_5s]$, has degree 1.

Step <ref>: The generating set induces a Khovanskii basis $\mathcal{B}=\{p_1p_2s,~p_1p_3s,~p_1p_4s,~p_1p_5s\}$ for the associated graded algebra $\rv$.

Step <ref>: The value semigroup is generated by the valuations of the Khovanskii basis elements, namely: $(2,2,1)$, $(2,3,1)$, $(3,2,1)$, and $(3,3,1)$.
It is straightforward to see that $(0,1,0) = (2,3,1) - (2,2,1)$ and $(1,0,0)=(3,2,1)-(2,2,1)$ are in $G_0(S)$ and generate $\Z^2\oplus \{0\}$.
So, $\ind(\rv)$ is 1:

Step <ref>: Then $\nobv$ is the square $\conv\{(2,2),~(2,3),~(3,2),~(3,3)\}$ in $\R^d$, with volume 1.

Step <ref>: The self-intersection index $[V,\dots,V]$ is $d!\deg(\varphi_V) \vnob /~\ind{\rv}= 2!\cdot 1 \cdot 1 / 1 = 2$.

Output: The maximum number of positive steady states (outside the base locus) is at most 2.

The base locus consists only of purely complex solutions, hence no positive solutions are missed by the Newton-Okounkov body bound.

In contrast, by Definition <ref>, the mixed volume of the chemical reaction system is 18.
We already saw that the maximum number of isolated positive steady states is 1.
This is a great example of the volume of the self-intersection index achieving a tighter bound than the mixed volume: the mixed-volume overcount is $18 - 1 = 17$, but the analogous Newton-Okounkov overcount is $2 - 1 = 1$.

Table <ref> summarizes the bounds for each of the networks considered in this section. We observe that the Newton-Okounkov bound is always better (or at least as good) as the mixed-volume bound.

nickname,     example
max number      max number      mixed-volume     Newton-Okounkov

[source]         steady states      $\C$-steady-states      bound      body bound

Wnt <cit.>     Example <ref>     $\geq 3$     9     56     32

KST <cit.>     Example <ref>
2     4     4     2

PM <cit.>     Example <ref>
1     4     4     4

BCY <cit.>     Example <ref>
1     2     18     2

The various root counts associated to the chemical reaction systems in Section <ref>. The Newton-Okounkov body bound matches or outperforms the mixed-volume bound for every example.

§ DISCUSSION

The self-intersection index of a finite dimensional subspace of rational functions $V$ bounds the number of isolated
solutions to a general polynomial system in the linear span of a Khovanskii basis for the associated graded algebra $\rv$.
However, in practice, very few examples of Newton-Okounkov bodies had
been computed.
In this paper,
we defined the first notion of a Newton-Okounkov body of a chemical reaction system, and gave a procedure for computing this body and the resulting intersection index  (Procedure <ref>).
Applying this procedure on
concrete examples, we showed that the volume of the intersection index of a chemical reaction network
the mixed volume bound (Section <ref>).
These examples provide a proof-of-concept for our new method.
In summary, the Newton-Okounkov body is a promising new tool for achieving our ultimate goal of tight bounds on the number of steady states of an arbitrary reaction network.
We conclude with some final remarks and ideas for next steps.

§.§ When can a Newton-Okounkov body improve bounds?

There are two important situations when a Newton-Okounkov body achieves better upper bounds on the maximum number of steady states of chemical reaction networks than those discussed in Section <ref>.

First, we recall that the Kushnirenko Theorem and the Bernstein Theorem take into account only the monomial support and are independent of precise rate constants in the chemical reaction system (<ref>) (provided that the rate constants are sufficiently general).
However, in certain applications or for specific networks, the rate constants are likely to be non-generic.
For instance, we may
a priori know some relationship among the rate constants, such as equality.
On the other hand, we may have specific values for some or all of the rate constants.
For example, the chemical reaction system from Example <ref> has rate constants $\kappa =
(k,k,5k,k,k,5k,k,k,k,k,k,k,k,k,5k,k,5k,k,k,k)$, for $k\in\R_{\geq 0}$.
Recall that for this example, the Newton-Okounkov body bound greatly outperformed the mixed-volume bound.
We anticipate the theory of Newton-Okounkov bodies for chemical reaction networks will improve bounds on the maximum number of steady states in similar cases where we can utilize extra information about the rate constants.

Second, the intersection index of a subspace $V$ counts solutions to a polynomial system drawn from $V$ outside an algebraic variety
contains the base locus and the set of poles of $V$.
The mixed volume on the other hand counts all isolated solutions in $(\C^*)^d$.
Hence, using the theory of Newton-Okounkov bodies, we can restrict to counting solutions in some specified Zariski-open subset of $(\C^*)^d$, for example.
With an appropriately chosen basis for $V$,
the theory of Newton-Okounkov bodies for chemical reaction networks can bound steady states outside the specified base locus.

§.§ Non-conservative networks

The examples in Section <ref> were all pathological in the sense that none of the networks were conservative.
Recall that for a conservative network, each species occurs in at least one nonnegative conservation law.
Based on our investigations, we conjecture the following.

For every conservative reaction network, the mixed-volume bound is equal to the Newton-Okounkov bound.
In other words, conservative networks are “Bernstein-general”.

§.§ Non-polynomial chemical reaction systems

We remark that the theory developed in Section <ref>  supports systems of the form $p_1/q_1=0,\dots,p_d/q_d=0$, where $p_i,q_i$ are polynomials in $\C[x_1,\dots,x_d]$.
That is, the self-intersection index is defined for any nonzero finite dimensional vector space $V$ of rational functions on a projective variety $X$, and hence
bounds the number of effective solutions to systems of rational equations.
We considered here chemical reaction systems arising from mass-action kinetics.
Other kinetic assumptions, or other steady-state approximations, produce alternative chemical reaction systems that are not polynomials.

How tight are Newton-Okounkov body bounds associated to rational steady-state approximations, e.g., Michaelis–Menten kinetics, or associated to positive rational steady-state parametrizations (c.f., <cit.>)?

§.§ Best bounds

Procedure <ref> detailed a method for computing one possible Newton-Okounkov body for a chemical reaction system.
However, the resulting Newton-Okounkov body is not well-defined.
There are a number of choices, including choice of valuation/monomial term order and choice of a vector space basis, made in Procedure <ref> that can potentially yield different associated Newton-Okounkov bodies.
For instance, the defining polynomial system in (<ref>) may in principle give different vector spaces $V$ depending on the order of the variables.
We would like one well-defined definition of a Newton-Okounkov body, or at least some detailed guiding principles for how best to define one to achieve the tightest bound.
More explicitly, a next research step is to remove the “trial-and-error” in applying the Procedure <ref>; see Remark <ref>).

Under which hypotheses on a basis for $V$ does the associated Newton-Okounkov body bound give the tightest bound on the maximum number of steady states of a chemical reaction system?

In this work, we showed how Newton-Okounkov body theory can and should be used in practice to obtain good bounds on the number of solutions to polynomial systems arising from applications.
Moving forward, we are hopeful that Procedure <ref> will be an important computational tool for applied mathematicians.
A researcher interested in counting solutions to their own polynomial system can make simple modifications of our code and subsequently calculate the associated self-intersection index, thereby obtaining better bounds than was previously possible.

§.§ Acknowledgements

Ideas leading to this research emerged during the Newton-Okounkov Bodies Working Seminar organized by Frank Sottile and Emanuele Ventura at Texas A&M University in Spring 2018.
When this research was initiated,
NO was supported by an American Fellowship from the AAUW.
We thank
Michael Burr,
Heather Harrington,
Anne Shiu,
Aleksandra Sobieska,
Frank Sottile,
Ruth Williams
for helpful discussions.

§ NOTATION

We collect the notation used throughout the paper.

* $\N$, the natural numbers, including 0

* $e_i\in \R^d$, the $i$-th standard basis vector

* $F^*:=F\setminus \{0\}$, the set of nonzero elements of $F$

* $X$, an irreducible $d$-dimensional complex algebraic variety

* $\C(X)$, function field of $X$

* $V\subset \C(X)$, a nonzero finite-dimensional complex vector space of rational functions on $X$

* $s$, a formal grading variable

* $V^k\subset \C(X)$, the subspace spanned by all $k$-fold products of elements in $V$

$\rv:=\displaystyle \bigoplus_{k\geq 0} V^ks^k$, the graded subalgebra of $\C(X)[s]$ associated to $V$

* $\nu$, a valuation on $\C(X)$

* $\widehat{\nu}$, a “lifted” valuation on $\rv$

* $S:=S(V,\nu)=\widehat{\nu}(\rv)\subset \Z^d\oplus \Z$, the value semigroup associated to $(V,\nu)$

* $\mathcal{B}\subset \rv$, a finite Khovanskii basis for $\rv$

* $G(\rv)=\Z S(V,\nu)$, the group generated by $S(V,\nu)$

* $\pi : \R^{d+1} \rightarrow \R$ given by projection onto the last coordinate, i.e., $(x_1, \dots, x_d, x_{d+1})\mapsto x_{d+1}$

* $G_0(\rv) := \pi^{-1}(0) \cap G(\rv)$, the subgroup of $\Z^n\bigoplus \{0\}$ generated by $S(V,\nu)$

* $\text{ind}(\rv)$, the index of $G_0(\rv)$ as a subgroup of $\Z^d \times \{0\}$

* $\Delta(V,\nu,
\mathcal{B})$, the Newton-Okounkov body associated to $V$, given valuation $\nu$ and Khovanskii basis $\mathcal{B}$