\section{Introduction}
Post-translational modification (PTM) of proteins (for instance, by
phosphorylation or methylation)
plays a key role 
in signal transduction, with disruptions having
implications for human health~\cite{cohen-role,Anas}. 
Much research 
therefore has focused on determining how networks of PTM's
process and encode information~\cite{escape,SK}. 
The goal, on the theory side, is to understand the dynamical systems
arising from such networks.  Indeed, qualitative properties,
such as the capacity for bistability or oscillations, 
  may indicate that the underlying biochemical mechanism acts as a
  biological switch or clock \cite{tyson-albert}. 
Bistability~\cite{bistab-cascade,mapk-bistable} and oscillations~\cite{yeast-mapk-oscillations,oscillations-mapk-cancer}, have been observed in cellular signaling networks, such as 
mitogen-activated protein kinase (MAPK) cascades, in which PTM networks play a key role.
 
  
Although PTM networks are often large, with 
parameter values that are difficult to estimate, 
a growing body of research has achieved success in 
determining whether a given network admits, for instance, bistability or oscillations, and even, in some cases, obtaining a mathematical description of
the parameter regions that give rise to such properties.

The aim of this Perspective is to highlight these recent mathematical developments.


The basic building block of a PTM network is formed by the following two sets of reactions:
  \begin{equation}
    \label{eq:basic_network}
    S + M \leftrightarrows S M \to S^* + M   \quad \text{ and }  \quad  S^*+U \leftrightarrows S^* U \to S + U   ~,
  \end{equation}
in which a substrate $S$ forms a complex

 
  with a modifier $M$, and then $M$ and the modified substrate~$S^*$ 
dissociate 
  (left); this modification can be undone by another modifier $U$
  (right). 
  A PTM network is adapted from Eq.~\ref{eq:basic_network}
  by adding known biological interactions: for instance, the complexes $S M
  $ and $S^* U$ might undergo more modifications, there
  might be many substrates and modifiers, and a (modified)
  substrate of one reaction might be a modifier of another. 


Arguably the most detailed and accurate description of the dynamics
of a PTM network is achieved at the level of mass-action. 
Here we focus on  systems in which the  spatial distribution of
concentrations can be neglected
  and where models therefore take the form of ordinary differential
  equations. Most models in the literature are of this type. 
  (When spatial effects may not be neglected  \cite{Kholodenko2006,KHOLODENKO901,Kholodenko2010}, partial differential equations
  must be considered and
the connection between network structure and qualitative properties
has not yet been studied to the same extent; for some recent results,
see \cite{BifTheo-020,CasianMethodLines,diff-001}.)
  
  Even if spatial effects are neglected, the resulting dynamical
  models (systems of ODEs) 
contain many variables and
  parameters  and are therefore difficult to analyze.
  Biophysicists and biochemists have long developed tools to deal with
  this situation, in particular the King-Altman method and the
  Michaelis-Menten reduction under the quasi-steady-state assumption
  proposed by Briggs and Haldane (see 
  \cite[Chapters~2\&4]{cornish2012fundamentals} and the references
  therein). 
  These methods guide the systematic elimination of
  variables and thus reduce the dimension of a
  system. 
It is important to note that the steady states of these reduced systems are in one-to-one
  correspondence with those of the original system \cite{feliu-wiuf-crn},
  while this need not be the case for the time-dependent trajectories of
  the respective ODE systems. 

In recent years,
  the abundance of PTM networks in biology and the accompanying
  need to analyze large networks has stimulated further research on
  the systematic elimination of variables. Much of this research
  originated at the interface of Chemical Reaction Network Theory,
  mathematical biology, and algebraic geometry. 
 Some of the results obtained thereby are a natural extension of the 
King-Altman method and to some extent the Michaelis-Menten reduction. 


 In Section~3, we will summarize
 the most consequential such elimination results for PTM
  networks
(after introducing the mathematical
  formalism used to study PTM networks in Section~2).

  In Section~4 we will turn to a specific class of PTM networks, those
  describing the multisite phosphorylation/dephosphorylation of a single substrate.
Much has been discovered in recent years about 
the steady states and the dynamics of these systems --
  more than for general PTM networks. 
What has been revealed is a close connection between a network's structure and its qualitative dynamics.  Consider, for instance, 
  the sequential, double phosphorylation (also called dual-site phosphorylation) of a protein 
in which the enzymatic mechanism is 
fully {\em processive} versus those with fully {\em distributive} mechanism (see Figure~\ref{fig:proc-dist}). 
 Processive systems are globally stable -- convergent to a unique steady state -- while 
 distributive systems may be bistable.  
These properties -- and those of several related networks -- are summarized later in Figure~\ref{fig:multisite-aims}.

\begin{figure}[hbt]
\begin{center}
	\begin{tikzpicture}[scale=1.7]
\draw (0,3) rectangle (1, 3.4);
\draw (2,3) rectangle (3, 3.4);
\draw (4,3) rectangle (5, 3.4);
\draw (6,3) rectangle (7, 3.4);
\draw (0,2.3) rectangle (1, 2.7);
\draw (2,2.6) rectangle (3, 3);
\draw (4,2.6) rectangle (5, 3);
\draw (6,2.3) rectangle (7, 2.7);
\draw (4,3.5) circle (0.1);	
\draw (6,3.5) circle (0.1);	
\draw (6.3,3.5) circle (0.1);	
    	\node[] at (0.5, 3.2) {Substrate};
    	\node[] at (2.5, 3.2) {Substrate};
    	\node[] at (4.5, 3.2) {Substrate};
    	\node[] at (6.5, 3.2) {Substrate};
    	\node[] at (0.5, 2.85) {$+$};
    	\node[] at (0.5, 2.5) {Kinase};
    	\node[] at (2.5, 2.8) {Kinase};
    	\node[] at (4.5, 2.8) {Kinase};
    	\node[] at (6.5, 2.85) {$+$};
    	\node[] at (6.5, 2.5) {Kinase};
    	\node[] at (4, 3.5) {P};
    	\node[] at (6, 3.5) {P};
    	\node[] at (6.3, 3.5) {P};
	 \draw[->] (1.2, 3.25) -- (1.8, 3.25);	
	 \draw[<-] (1.2, 3.15) -- (1.8, 3.15);	
	 \draw[->] (3.2, 3.2) -- (3.8, 3.2);	
	 \draw[->] (5.2, 3.2) -- (5.8, 3.2);
\draw (-0.2,3.7) rectangle (7.2,2);
\node[] at (7.9, 3.6) {{\sc Processive}};
\node[] at (8.5, 1.85) {{\sc Distributive}};
\draw (0,1) rectangle (1, 1.4);
\draw (2,1) rectangle (3, 1.4);
\draw (4,1) rectangle (5, 1.4);
\draw (6,1) rectangle (7, 1.4);
\draw (8,1) rectangle (9, 1.4);
\draw (0,0.3) rectangle (1, 0.7);
\draw (2,0.6) rectangle (3, 1);
\draw (4,0.3) rectangle (5, 0.7);
\draw (6,0.6) rectangle (7, 1);
\draw (8,0.3) rectangle (9, 0.7);
    	\node[] at (0.5, 1.2) {Substrate};
    	\node[] at (2.5, 1.2) {Substrate};
    	\node[] at (4.5, 1.2) {Substrate};
    	\node[] at (6.5, 1.2) {Substrate};
    	\node[] at (0.5, 0.85) {$+$};
    	\node[] at (0.5, 0.5) {Kinase};
    	\node[] at (8.5, 1.2) {Substrate};
    	\node[] at (8.5, 0.85) {$+$};
    	\node[] at (8.5, 0.5) {Kinase};
    	\node[] at (2.5, 0.8) {Kinase};
    	\node[] at (4.5, 0.85) {$+$};
    	\node[] at (4.5, 0.5) {Kinase};
    	\node[] at (6.5, 0.8) {Kinase};
\draw (4,1.5) circle (0.1);	
\draw (6,1.5) circle (0.1);	
\draw (8,1.5) circle (0.1);	
\draw (8.3,1.5) circle (0.1);	
    	\node[] at (4, 1.5) {P};
    	\node[] at (6, 1.5) {P};
    	\node[] at (8, 1.5) {P};
    	\node[] at (8.3, 1.5) {P};	
	 \draw[->] (1.2, 1.25) -- (1.8, 1.25);	
	 \draw[<-] (1.2, 1.15) -- (1.8, 1.15);	
	 \draw[->] (3.2, 1.2) -- (3.8, 1.2);	
	 \draw[->] (5.2, 1.25) -- (5.8, 1.25);	
	 \draw[<-] (5.2, 1.15) -- (5.8, 1.15);	
	 \draw[->] (7.2, 1.2) -- (7.8, 1.2);
\draw (-0.2,1.7) rectangle (9.2,0);
    	\end{tikzpicture}
\end{center}
\caption{Dual-site phosphorylation may be {\em processive} (top) or {\em distributive} (bottom) or somewhere in between; cf.\ \cite[Figure 2B]{salazar}.}
\label{fig:proc-dist}
\end{figure}

In this Perspective, and especially at the close, we will highlight important future directions.















 \subsection{Related review articles}
Excellent reviews on various aspects of PTM networks
have been published (which we highlight below).  Our goals therefore are complementary to what earlier reviews accomplished and, to some extent, pick up where they left off.  Indeed, there is currently no review summarizing recent progress on elimination for PTM networks or the dynamics of
 multisite phosphorylation networks involving a single substrate.

As mentioned earlier, we are interested in the question of how PTM networks process information.  
Theoretical aspects of this question, and experimental approaches also, were reviewed
by Prabakaran {\em et al.}~\cite{escape}.  
Our Perspective, in contrast, is limited to the theory side, and specifically to the 
question of how a PTM network constrains the resulting dynamics (which in turn constrain the capacity for information encoding).

Another related review is that of Salazar and H\"ofer~\cite{salazar}, which describes 
the many biological functions of (multisite) phosphorylation (the best-studied PTM), a thorough list of possible mechanisms, and the resulting kinetic models.  
Their review,
to abridge its title, goes ``from molecular mechanisms to kinetic models'' --
whereas, here we go from kinetic models to dynamics.  
  To be clear, Salazar and H\"ofer did examine dynamics (including
  ultrasensitivity and bistability but not oscillations), but not for all
  the systems described here. 

Finally, much of the interest in phosphorylation networks is due to their role in MAPK cascades.  The dynamics of MAPK cascades was reviewed by Hell and Rendall in~\cite{hell-2015}.


\section{Background} \label{sec:background}
{\em Post-translational modification}
is a biochemical alteration of a protein that occurs after its mRNA has been translated into a sequence of amino acids.
The most common PTM is {\em phosphorylation}, the enzyme-mediated addition of a phosphate group to a protein substrate, which, to quote from~\cite{25}, ``can alter its behavior in almost every conceivable way.''
  The basic phosphorylation/dephosphorylation mechanism is a relabeled version of the building block in Eq.~\ref{eq:basic_network}:
  \begin{equation}
    \label{eq:basic_phospho}
    S_0 + K \leftrightarrows S_0 K \to S_1+K \quad \text{ and }  \quad S_1 + F \leftrightarrows S_1 F \to S_0+F~,
  \end{equation}
  where a {\em kinase} and a {\em phosphatase} enzyme ($K$ and $F$) act as
  modifiers,
the non-phosphorylated (unmodified) protein is
  denoted by $S_0$, and the phosphorylated (modified) protein is $S_1$. 
The mechanism in Eq.~\ref{eq:basic_phospho} also describes other PTMs, such as acetylation, the
addition of an acetyl group.  Although some PTMs are non-binary (more
than one modification is possible, even at a single
site~\cite{escape}), we will focus on PTMs built on
the mechanism in Eq.~\ref{eq:basic_phospho} and use, for ease,
the language of
phosphorylation.

So far we have discussed phosphorylation on a single site.  
Many substrates, however, have more than one site at which phosphate
groups can be attached 
(in which case we continue to use subscripts to denote the number of phosphorylations, 
e.g., $S_2$ is a doubly phosphorylated protein).



Such multisite phosphorylation
may be 
  {\em sequential} or {\em random}, or somewhere in between, 
and {\em processive} or {\em distributive}, or somewhere in
between~\cite{Guna_threshold, PM, salazar}. 
  In {\em sequential} phosphorylation
  the last phosphate group to be attached is the first to be removed, so,
  for example, 
  there is a unique once-phosphorylated substrate ($S_1$). This Perspective focuses solely on sequential phosphorylation. 
(In {\em random} phosphorylation, phosphate groups can be added in any order.)
{\em Processive} phosphorylation, as depicted in Figure~\ref{fig:proc-dist},
requires an enzyme and 
substrate to bind only once in order to add 
  or remove 
all phosphate groups.
{\em Distributive} phosphorylation, in contrast, requires multiple bindings
to add 
  or remove 
all phosphate groups, because 
   at most one group is added with every binding of enzyme
  and protein. 



\begin{example}[Mixed-mechanism network] \label{ex:mixed}
The following is a {\em mixed-mechanism} dual-site network (of Figure~\ref{fig:multisite-aims}(E)), in which phosphorylation is processive and dephosphorylation is distributive~\cite{SK}:
\begin{equation}
  \label{eq:mixed-network}
  \begin{split}
  \xymatrix{
   
    S_0 + K  \ar @<.4ex> @{-^>} [r] ^-{k_1}
    &\ar @{-^>} [l] ^-{k_{2}} 
    S_0 K \ar 
   
    [r] ^-{k_3} 
    & S_1 K \ar
    [r] ^-{k_4}
    & S_2 + K\\
   
    S_2 + F \ar @<.4ex> @{-^>} [r] ^-{k_{5}}
    &\ar @{-^>} [l] ^-{k_6} S_{2} F 
    \ar
    [r] ^-{{k_7}} 
   
   
   
   
   
   
    &S_1 + F
    \ar @<.4ex> @{-^>} [r] ^-{k_{8}}
    &\ar @{-^>} [l] ^-{k_9}
    S_{1} F \ar
    [r]^-{{k_{10}}} & S_0 + F
  }
  \end{split}
\end{equation}


Under the commonly used assumption that the (cellular) volume is constant, 
the network in Eq.~\ref{eq:mixed-network} gives rise, via mass-action kinetics,
to the following ODEs:
\begin{align} 
    \notag
 \dot x_1 &~=~     -k_1 x_1 x_2+k_2 x_3+k_{10} x_9 \\
    \notag
 \dot x_2 &~=~      -k_1 x_1 x_2+k_2 x_3+k_4 x_4 \\
    \notag
 \dot x_3 &~=~      k_1 x_1 x_2-(k_2+k_3) x_3  \\
    \notag
 \dot x_4 &~=~      k_3 x_3- k_4 x_4 \\
    \label{eq:OEs-mixed}
 \dot x_5 &~=~      k_4 x_4-k_5 x_5 x_6+k_6 x_7 \\
    \notag
 \dot x_6 &~=~      -k_5 x_5 x_6-k_8 x_8 x_6+(k_6+k_7) x_7+(k_9+k_{10}) x_9 \\
    \notag
 \dot x_7 &~=~      k_5 x_5 x_6-(k_6+k_7) x_7 \\
    \notag
 \dot x_8 &~=~      k_7 x_7-k_8 x_6 x_8+k_9 x_9 \\
    \notag
 \dot x_9 &~=~      k_8 x_6 x_8-(k_9+k_{10}) x_9~,
\end{align}
where $x_1,x_2, \ldots, x_9$ denote, respectively, the concentrations of the species
$S_0$, $K$, $S_0K$, $S_1 K$, $S_2$, $F$, $S_2 F$, $S_1$, $S_1 F$.  The total amounts of free and bound enzyme or substrate, denoted by 
$(K_{\mbox{tot}}, F_{\mbox{tot}}, S_{\mbox{tot}} )\in \mathbb{R}_{>0}^3$,
remain constant as the dynamical system in Eq.~\ref{eq:OEs-mixed} progresses, so we obtain the following {\em conservation laws}:
\begin{align}
\label{eqn:conservation}
K_{\mbox{tot}} ~=~ x_2 + x_3 + x_4~, ~ 
F_{\mbox{tot}} ~=~ x_6+x_7+x_9  ~, ~
S_{\mbox{tot}} ~=~ x_1+x_3+x_4+x_5+x_7+x_8+x_9  ~.
\end{align}
Thus, each trajectory $x(t)$ of Eq.~\ref{eq:OEs-mixed} is confined to
an {\em invariant set} defined by some $(K_{\mbox{tot}}, F_{\mbox{tot}}, S_{\mbox{tot}} )$: 
\begin{align} \label{eq:invt-set}
\{x \in \mathbb{R}_{\geq 0}^9 \mid 
    \text{ the conservation equations $\eqref{eqn:conservation}$ hold} \}~.
\end{align}
\end{example}

One focus of this Perspective is on the dynamics arising from 
  multisite phosphorylation networks such as the mixed-mechanism
  network of Example~\ref{ex:mixed} and those displayed in
  Figure~\ref{fig:multisite-aims}. We report what is known about
  their qualitative dynamics: we describe, whether or not a network
  is {\em bistable}, i,e. whether or not it has the capacity for two
  or more steady states, stable to small perturbations, within some
  invariant set in Eq.~\ref{eq:invt-set}.
(The biological significance of bistability has been described, for
instance, in \cite{thomson2009unlimited}). If not, is the network {\em
  globally stable}, with a unique steady state in each invariant set to which
all trajectories converge? Or, does the network admit {\em
  oscillations}, i.e., periodic orbits?
To answer such questions, one key technique is elimination, which we discuss next.





\begin{figure}[hb]
\begin{center}
	\begin{tikzpicture}[scale=1]
    	\node[above left] at (0,0.25) {(E)};
    	\node[above left] at (0,3.25) {(A)};
    	\node[above left] at (0,1.75) {(C)};
    	\node[above right] at (15,1.75) {(D)};
    	\node[above right] at (15,3.25) {(B)};
	\draw (0,0) rectangle (6.5,1.2);
    	\node[above] at (3.2,0.5) {Mixed-mechanism dual-site};
    	\node[above] at (3.2,0) {(unique steady state, oscillations)};
	\draw (0,1.5) rectangle (6.5,2.7);
    	\node[above] at (3.2,2) {ERK-mechanism dual-site};
    	\node[above] at (3.2,1.5) {(bistable, oscillations)};
	\draw (0,3) rectangle (6.5,4.2);
    	\node[above] at (3.2,3.5) {Fully {distributive} dual-site};
    	\node[above] at (3.2,3) {(bistable)}
	\draw (10,1.5) rectangle (15,2.7);
    	\node[above] at (12.5,2) {Fully {processive} dual-site};
    	\node[above] at (12.5,1.5) {(globally stable)};
	\draw (10,3) rectangle (15,4.2);
    	\node[above] at (12.5,3.5) {Michaelis-Menten dual-site};
    	\node[above] at (12.5,3) {(bistable, no oscillations)};
    	\node[above] at (8.3,3.2) {{\em \dots reduces to \dots}};
    	\node[above] at (8.3,2) {{\em \dots reduces to \dots}};
 	\end{tikzpicture}
\end{center}
\caption{Dual-site phosphorylation/dephosphorylation networks, also called {\em futile cycles}: five types and their properties. A network reduces to another if the latter approximates the former if certain conditions on the parameters hold.  See \S \ref{sec:a}--\S \ref{sec:e}.} 
\label{fig:multisite-aims}
\end{figure}




\section{Elimination and the role of intermediate complexes} \label{sec:intermediates}
We begin this section with the well-known reversible Michaelis-Menten mechanism
  (cf.\ \cite{cornish2012fundamentals}):\\
  \begin{minipage}[t]{0.4\linewidth}
    \underline{Reaction Network}
    \begin{align*}
      &\xymatrix{
        S_{0}+ K \ar @<.4ex> @{-^>} [r] ^-{k_1}
      &\ar @{-^>} [l] ^-{k_2}  S_{0}K \ar @<.4ex> @{-^>} [r] ^-{k_3}
      &\ar @{-^>} [l] ^-{k_4}   S_{1}+K 
        }\\
      \intertext{\underline{Conservation laws}}
      & x_1 + x_3 + x_4 = S_{\mbox{tot}}  \\
      & x_2 + x_3 = K_{\mbox{tot}}
    \end{align*}
  \end{minipage}
  \begin{minipage}[t]{0.55\linewidth}
    \underline{ODEs (using $x_1$ for $S_0$, $x_2$ for $K$, $x_3$ for
      $S_0K$, and $x_4$ for $S_1$):}
    \begin{align*}
      \dot x_1 &= -k_{1} x_{1} x_{2}+k_{2} x_{3} \\
      \dot x_2 &= -k_{1} x_{1} x_{2}+(k_{2} +k_{3}) x_{3}-k_{4}
                 x_{2} x_{4} \\ 
      \dot x_3 &=  \phantom{-}k_{1} x_{1} x_{2}-(k_{2}+k_{3}) x_{3}+k_{4} x_{2}
                 x_{4} \\
      \dot x_4 &= \phantom{-}k_{3} x_{3}-k_{4} x_{2} x_{4}
    \end{align*}
  \end{minipage}\\

The total amount of substrate $S_{\mbox{tot}}$ is
often much
larger than
that of the enzyme $K_{\mbox{tot}}$
  (Briggs-Haldane assumption \cite{cornish2012fundamentals}), and one can
  therefore assume $\dot x_3 \approx 0$. It is then standard practice
  to use the equation $\dot x_3 = 0$
and the
 relation $ x_2 + x_3 = K_{\mbox{tot}} $

 to \emph{eliminate}  the variables $x_2$ and $x_3$.
That is, one solves for $x_2$ and $x_3$ in terms
  of $x_1$, $x_4$, $K_{\mbox{tot}}$, and the rate constants (parameters): 

  \begin{align} \label{eq:param-mm}
    x_2 &= \frac{K_{\mbox{tot}} (k_{2}+k_{3})}{k_{2}+k_{3}+k_{1} x_{1}+k_{4} x_{4}}, &
    x_3 &=\frac{K_{\mbox{tot}} (k_{1} x_{1}+k_{4}x_{4})}{k_{2}+k_{3}+k_{1}
          x_{1}+k_{4} x_{4}}~, \\
    \intertext{and substitutes these expressions into
    $\dot x_1$ and
    $\dot x_4$ to obtain the ODEs of the 
   
    reduced
    system:}
    \dot x_1 &= \frac{K_{\mbox{tot}} (-k_{1} k_{3} x_{1}+k_{2} k_{4}
               x_{4})}{k_{2}+k_{3}+k_{1} x_{1}+k_{4} x_{4}} = -\dot x_4~. \notag
  \end{align}%
Mathematically, Eq.~\ref{eq:param-mm} is a {\em parametrization} of $x_2$ and $x_3$ in terms of $x_1$, $x_4$, $K_{\mbox{tot}}$,
  and rate constants.


One can apply similar reasoning to PTM networks, removing
intermediate complexes like $S_0K$ to obtain a reduced system with fewer variables. 
Although this process can 
become arduous for large networks, the well-known King-Altman method
provides a systematic way to eliminate intermediates
\cite{cornish2012fundamentals}.
However, this method
relies at least implicitly on the
  ability to identify the enzyme (to guide the choice of
  variables to solve for). Furthermore, it does not explain the
  positivity of the parametrization. The first requirement, in
  particular, limits the applicability of the method as
 it is, frequently the case (for
  example, in cellular signaling) that the
  substrate of one modifier acts as a modifier for a second
  substrate. Such cascades cannot be handled by the King-Altman
  method. 

Fortunately, both concerns were resolved recently.
First, positivity in the King-Altman framework 
was explained
  by  Thomson and Gunawardena  \cite{TG}, who initiated the theory underlying elimination in PTM networks (see also~\cite{linear-framework}). 
Subsequently, Feliu and Wiuf extended this to address the second concern: how to eliminate 
in cascades~\cite{feliu-signaling} and other general networks~\cite{feliu-wiuf-crn,
feliu-wiuf-ptm} (see also \cite{saez-gph}).
These works describe procedures to systematically eliminate
  intermediates and some substrates.  
In fact, the resulting parametrizations arise from spanning
  trees of certain graphs -- much like King-Altman!
 
  


\subsection{Elimination and steady-state parametrizations}
  For the Michaelis-Menten mechanism above, elimination
  was guided by biophysical insight as well as the aim to simplify the
  system and understand the dynamics (of the simplified
  system). If one were interested only in the positive steady states
  of the original ODEs, or in the number of steady states, 
then one might pursue a different approach which also uses a parametrization.  
Specifically, we set to zero the right-hand sides of the differential equations, and then 
solve for each $x_i$ in terms of just $x_2$ and $x_4$,
and see that the set of all steady states is the 2-dimensional image of the map

$ \mathbb{R}^2_{>0}  \to \mathbb{R}^{4}_{>0} $ given by:  
  \begin{align} \label{eq:param}
    (x_2, x_4) ~\mapsto ~ \left( \frac{k_{2} k_{4} }{k_{1} k_{3}} x_{4},~ x_2, ~
    \frac{k_{4} }{k_{3}} x_{2} x_{4},~ x_4\right)~.
  \end{align}  
This {\em steady-state parametrization} is particularly simple, as it
  consists of the \emph{monomials} $x_2$, $x_4$, and $x_2x_4$ (or their scalar multiples). It is
  consequently called a \emph{monomial parametrization}.


The existence of a positive parametrization is guaranteed for PTM networks
  \cite{feliu-wiuf-ptm,messi,TG}. 
 

  However, such a parametrization need not be monomial.

As for the number of steady states, this corresponds mathematically to 
counting the number of times the image of the parametrization intersects some invariant set. 
If at least one such intersection (for some choice of rate constants) has at least two elements (i.e., at least two positive steady states in some invariant set), then we say that the network is {\em multistationary}.  Deciding whether a network is multistationary is in general difficult~\cite{mss-review}, 
but this can be accomplished easily when the parametrization is monomial
  \cite{signs,TSS}. Such an argument, for instance, was used to prove
  that 
  a reversible version of
the following {\em fully processive, $n$-site phosphorylation network} (for which the $n=2$ case is the network of Figure~\ref{fig:multisite-aims}(D)) is {\em not} multistationary~\cite{ConradiShiu}:
\begin{equation}
\label{eq:n-site-processive}
\begin{split}
  \xymatrix{
    S_0 + K  \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_0 K \ar [r] ^-{}
    &S_1 K \ar [r] ^-{}
    &\hdots \ar  [r] ^-{}
    & S_{n-1} K \ar [r] 
    ^-{}
    & S_n + K\\
    S_n + F \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_{n} F \ar  [r] ^-{} 
    ^-{} 
    &\hdots \ar [r]
    ^-{}
    &S_2 F \ar  [r]
    ^-{}
    &S_1 F \ar [r]
    ^-{}
    &S_0 + F
  }
\end{split}
\end{equation}
Another approach to precluding multistationarity in this network 
is described below in Section~\ref{sec:remove-intermediates}.
  

  For large (not necessarily PTM) networks  it can be challenging to obtain a
  parametrization, as there is (currently) no
  procedure apart from computing spanning trees 
  (unless the network has more structure~\cite{messi}). 
 
 
 
  Guided by
  experience, we suggest the following approach: if $p$ is the number of conservation laws, then we
  expect to parametrize the steady states in
  terms of $p$ variables (there is no guarantee though). Therefore, using a computer algebra system like {\tt Mathematica} or 
  {\tt Maple}, solve for all but $p$ variables, eliminating all
  intermediates (and some substrates)~\cite{feliu-wiuf-ptm,TG}. If the result is not
  satisfactory, try another subset of variables.

Finally, we can express a steady-state parametrization in terms of 
biochemically meaningful parameters, namely, 
the catalytic constants of the enzyme and the $K_m$ values.
For instance, for the above parametrization (\ref{eq:param-mm}),
  there are two $K_m$ values: $K_{m,1} = \frac{k_2+k_3}{k_1}$ for the
  substrate $S_0$ and $K_{m,2} = \frac{k_2+k_3}{k_4}$ for the
  substrate $S_1$.  
We solve these $K_m$ equations for $k_1$ and $k_4$ and substitute these
expressions into the parametrization (\ref{eq:param-mm}), which becomes:
\begin{align*}
    (x_2, x_4) ~\mapsto ~ \left(  \frac{k_2 K_{m,1}}{k_3 K_{m,2}} x_{4},~ x_2, ~
   \left(1+\frac{k_{2}}{k_{3}}\right) \frac{1}{K_{m,2}} x_{2} x_{4},~ x_4\right)~.
  \end{align*}  
Such an approach was used to analyze how catalytic constants enable the emergence of 
bistability in the fully distributive dual-site phosphorylation system (the network of Figure \ref{fig:multisite-aims}(A))~\cite{a6maya}. 


\subsection{Intermediates and steady states} \label{sec:remove-intermediates}
Going beyond parametrizations, Feliu and Wiuf also
  analyzed elimination at the level of the reaction network, and its effect on steady states and their stability. Here
removal of intermediates yields a reduced
  reaction network, 

  and they proved
  the following:
 
 if the reduced network admits at least $N\geq 1$ (locally stable) steady states (for
  some parameter values), then the original network also has at 
least $N$ (locally stable) steady states (for some parameters) \cite{FeliuWiuf}. 
(See~\cite{BP-inher,ME_entrapped,mss-review} for other results that ``lift'' multistationarity/bistability from small networks to large.)




As an
example, we revisit the fully processive network in Eq.~\ref{eq:n-site-processive}.
By iteratively removing intermediates $S_iK$ and $S_iF$  \cite{FeliuWiuf}, we
obtain the reduced network:



\begin{equation}
  \label{eq:1-site}
  \begin{split}
  \xymatrix{
    S_0 + K  \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_0 K \ar [r] ^-{}
    & S_n + K\\
    S_n + F \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_{n} F \ar  [r] ^-{} 
    ^-{} 
    &S_0 + F
  }
\end{split}
\end{equation}
The network in Eq.~\ref{eq:1-site} is mathematically equivalent to the
single-site phosphorylation network (the $n=1$ case of the 
network in Eq.~\ref{eq:n-site-processive}), which is known to be globally
stable: each invariant set has a unique steady state, this steady
state is positive, and it is the global attractor of the invariant
set~\cite{AS}.  
Thus, as explained above,
we conclude that the network in Eq.~\ref{eq:n-site-processive}
  has at least one (locally stable) steady state.
  
In fact, we can say more about the network in
Eq.~\ref{eq:n-site-processive}: because it reduces to the globally stable network in Eq.~\ref{eq:1-site}, it too is globally stable~\cite{ConradiShiu,EithunShiu}.  We proved this via results on how the dynamics are affected when intermediates are incorporated into a network. We turn to this topic next.









\subsection{Intermediates and dynamics}
As mentioned above, scientists have long used approximations of reaction systems
 (e.g., Michaelis-Menten), in which steady states correspond exactly to those of the original system --  but trajectories need not have such a correspondence. 
In the case of a reaction system with intermediates, however, the trajectories are close to each other:
Cappelletti and Wiuf proved that the trajectories can be uniformly approximated on compact time intervals by the trajectories of the reduced reaction system obtained by removing intermediates~ \cite{approx-elim}.  

Similarly, global stability of a network obtained by removing intermediates, or by making some reactions irreversible, can imply that the original system also is globally stable.  
Marcondes de Freitas {\em et al.} proved such a result, using monotone systems theory~\cite{Angeli2010}, for removing intermediates~\cite{Freitas1,Freitas2}; and Ali Al-Radwahi and Angeli proved such a result for making reactions irreversible~\cite{ali-angeli}.  










\subsection{Future direction: intermediates and oscillations}
One focus of this section was on the question of how dynamics and steady states are affected when intermediates and/or reactions are added or removed.  What we presently can say about dynamics, however, is mostly limited to the situation of global stability.  We therefore are interested in other dynamics, such as oscillations.  

For instance, is the presence of a Hopf bifurcation preserved when reactions or intermediates are added or removed?  To pose a concrete question, we revisit the mixed-mechanism network in Eq.~\ref{eq:mixed-network}. Suwanmajo and Krishnan showed that it admits a Hopf bifurcation; {\em is this still true if any of the irreversible reactions is made reversible?}  Can the answer be explained by some general, yet-to-be-proven theorem that ``lifts'' Hopf bifurcations  when reactions or intermediates are added?

Some partial results in this direction were given recently by Banaji.
Regarding the question posed above, Banaji proved that making any irreversible reactions reversible indeed preserves oscillations -- as do several other operations that add species or reactions~\cite{banaji-inheritance}.




 

\section{Dynamics of multisite phosphorylation of a single substrate} \label{sec:dynamics}
This section highlights what is known -- and what is open -- about the dynamics of the dual-site phosphorylation systems in Figure \ref{fig:multisite-aims}.  
We are not attempting to be thorough -- there are many more multisite systems, and we point the  reader to~\cite{scaffold,Enzyme-sharing,compartments,jolley,kapuy,SK}.

  \subsection{Distributive $n$-site systems, including Figure \ref{fig:multisite-aims}(A)} \label{sec:a}
Most studies on the mathematics of multisite phosphorylation
have focused on
the case of a sequential and fully {\em distributive}
mechanism
 \cite{a6maya, Enzyme-sharing,KathaMulti,ManraiGuna,Markevich,MAPK,ortega,thomson2009unlimited,WangSontag,TSS}:
\begin{equation}
  \label{eq:distributive-network}
  \begin{split}
  \xymatrix{
   
    S_0 + K \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_{0} K \ar
    [r] ^-{{}} 
    &S_1 + K
    \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{}
    S_{1} K 
    \ar
    [r]^-{} 
    & S_2 + K
    \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{}
    \cdots
    \ar
    [r]^-{} 
    & S_n + K
    \\
   
    S_n+F 
    \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{}
    \cdots
    \ar @<.4ex>
    [r]^-{} 
    & S_2 + F \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{} S_{2} F \ar
    [r] ^-{{}} 
    &S_1 + F
    \ar @<.4ex> @{-^>} [r] ^-{}
    &\ar @{-^>} [l] ^-{}
    S_{1} F \ar
    [r]^-{} 
    & S_0 + F
  }
  \end{split}
\end{equation}

This network is multistationary for $n \geq 2$; this result can be proven via a monomial
  parametrization of the steady states \cite{KathaMulti,TSS}.  
As for the maximum number of positive steady states, this number is at  most $2n-1$; this was proven by Wang and Sontag using an elimination similar to those in Section~\ref{sec:intermediates} \cite{WangSontag}.   A rationale to
  obtain rate constants for which the system has the
  maximum number of positive steady states was described in
  \cite{FHC}. 




As for multiple {\em stable} steady states, the distributive {\em dual-site} network (the $n=2$ case of 
the network in Eq.~\ref{eq:distributive-network}) admits bistability~\cite{bistable}.  

\subsection{Distributive systems, Michaelis-Menten approximation: Figure \ref{fig:multisite-aims}(B)} \label{sec:M-M} 
The Michaelis-Menten (MM) approximation of the distributive dual-site 
network
is a two-dimensional ODE system~\cite{WangSontag08}.  
The validity of this approximation for phosphorylation systems has been called into question~\cite{salazar}.  Nevertheless, analyzing this system is valuable,
 as dynamical properties of the MM approximation
sometimes can be ``lifted'' to the fully distributive system~\cite{WangSontag08}.  Indeed, the MM approximation is bistable, and Hell and Rendall ``lifted'' bistability to the full system~\cite{bistable}.

Oscillations, however, are precluded in the MM approximation, because every trajectory converges to some steady state.  This was proven by Wang and Sontag~\cite{WangSontag08} via monotone systems theory, and then re-proven by Bozeman and Morales \cite{no-osc} via the simpler Bendixson's Criterion.  

\subsection{Processive systems, including Figure \ref{fig:multisite-aims}(D)} \label{sec:pro-lit} 
As discussed in \S \ref{sec:intermediates},
 multisite phosphorylation systems that are fully {\em processive} 
 are globally stable~\cite{ConradiShiu,EithunShiu}, and thus,
in contrast to distributive systems, do {\em not} admit bistability or oscillations.  
This result revealed that processive systems lack the more interesting information-processing capabilities of distributive systems -- an idea that existed in the experimental literature, but, until recently, lacked a proof.




\subsection{The ERK-mechanism system: Figure \ref{fig:multisite-aims}(C)} \label{sec:c}
Figure \ref{fig:multisite-aims}(C) refers to the network of ERK regulation by dual-site phosphorylation by MEK and dephosphorylation by MKP3, here
\begin{equation}
  \label{eq:erk-1}
  \xymatrix@C=2ex{
    S_{00}+K 
    \ar @<.4ex> @{-^>} [r] ^-{}
    &
    \ar @{-^>} [l] ^-{}
    S_{00}K 
    \ar [r]^-{}
    & 
    S_{01}K
    \ar [r] ^-{k_{\rm cat}}
    &
    S_{11}+K 
    &  
    S_{11}+F   
    \ar @<.4ex> @{-^>} [r] ^-{}
    &
    \ar @{-^>} [l] ^-{}
    S_{11}F \ar [r] ^-{}
    &
    S_{10}F
    \ar [r] ^-{l_{\rm cat}}
    & 
    S_{00}+F
    \\
    S_{01} K
    \ar @<.4ex> @{-^>} [r] ^-{k_{\rm off}}
    &
    \ar @{-^>} [l] ^-{}
    S_{01}+K
    & 
    &
    &
    S_{10}F     
    \ar @<.4ex> @{-^>} [r] ^-{l_{\rm off}}
    &
    \ar @{-^>} [l] ^-{}
    S_{10}+F
    &
    &
    \\
    S_{10} + K     
    \ar @<.4ex> @{-^>} [r] ^-{}
    &
    \ar @{-^>} [l] ^-{}
    S_{10} K 
    \ar [r] ^-{}  
    &
    S_{11}+E   
    & 
    &
    S_{01} + F 
    \ar @<.4ex> @{-^>} [r] ^-{}
    &
    \ar @{-^>} [l] ^-{}
    S_{01}F  
    \ar [r] ^-{}
    &
    S_{00}+F~.
    &
  }
  \end{equation}
Here $S_{ij}$ represents the substrate with $i$ phosphate groups attached at the first site and $j$ at the second.  This notation reflects the fact that this network is {\em non-sequential}: the first site to be phosphorylated (via $ S_{00}K \to S_{01}K$) is also the first to be dephosphorylated ($S_{11}F \to S_{10}F$).

The ERK network
admits bistability and oscillations~\cite{long-term}.  
However, when 
 $k_{\rm cat} >> k_{\rm off}$ and
 $l_{\rm cat} >> l_{\rm off}$, then the network reduces to 
a network {\em without} bistability or oscillations, namely, the fully processive network in Figure~\ref{fig:multisite-aims}(D). 
(That network comprises only the reactions in Eq.~\ref{eq:erk-1} above.)

Accordingly, Rubinstein {\em et al.} posed the question, {\em How do bistability and oscillations in the ERK network emerge from the processive limit?}
More precisely, when the ``processivity levels''
$\frac{k_{\rm cat}}{k_{\rm cat} + k_{\rm off}}$ and $ \frac{l_{\rm cat}}{l_{\rm cat} + l_{\rm off}}$ are arbitrarily close to 1, is the ERK network still bistable and oscillatory?




\subsection{Mixed-mechanism systems: Figure \ref{fig:multisite-aims}(E)} \label{sec:e}
Only recently have there been studies of mixed mechanisms (partially
processive, partially distributive)~\cite{Aoki}. One such network, in
which phosphorylation is processive and dephosphorylation is
distributive, was depicted earlier in Eq.~\ref{eq:mixed-network}
(Example~\ref{ex:mixed}).
This network is {\em not} multistationary (Suwanmajo and Krishnan proved this 
via an elimination procedure like those described in
Section~\ref{sec:intermediates})~\cite{SK}. 
It follows, by a standard application of the Brouwer fixed-point
theorem, that there is always a unique steady state (for all rate
constants and conservation-law values).  This proves half of a
conjecture that we posed~\cite{ConradiShiu}.

We were surprised when the other half of our conjecture was disproven:
in contrast with processive systems (\S \ref{sec:pro-lit}), mixed
systems need {\em not} be globally stable.  In fact, the network in
Eq.~\ref{eq:mixed-network} exhibits oscillations~\cite{SK}!   
The significance of this result, as explained by its discoverers,
Suwanmajo and Krishnan, is that the network in
Eq.~\ref{eq:mixed-network} ``could be the simplest enzymatic
modification scheme that can intrinsically exhibit
oscillation''~\cite[\S 3.1]{SK}. Accordingly, this network's capacity
for oscillations is an interesting future direction. 

  As a starting point,
  we display in Figure~\ref{fig:osci} such oscillations arising from the
  mixed-mechanism network.  The rate constants we used,
  from~\cite{SK}, are displayed in Table~\ref{tab:rcs}.
  Suwanmajo and Krishnan displayed oscillations for 
  $(K_{\mbox{tot}},F_{\mbox{tot}})=(17.5,5)$, while those in Figure~\ref{fig:osci} 
  used   $(K_{\mbox{tot}},F_{\mbox{tot}})=(100,8)$.  These results suggest that 
  oscillations exist over a wide range of $K_{\mbox{tot}}$ -- and (possibly) $F_{\mbox{tot}}$ -- values,
  extending what was found earlier when only $K_{\mbox{tot}}$ was allowed to vary;~\cite{SK}.



\begin{table}[!h]
  \centering
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline
    $k_{1}$ &  $k_{2}$ &  $k_{3}$ & $k_{4}$ &  $k_{5}$ &  $k_{6}$ & $k_{7}$ & $k_{8}$ &  $k_{9}$ &  $k_{10}$ \\ \hline 
    1 & 1 & 1 & 1 & 100 & 1 & $0.9$ & 3 & 1 & 100 \\ \hline
  \end{tabular}
  \caption{\label{tab:rcs}
    Rate constants, from~\cite[Supplementary Information]{SK}, leading to oscillations 
   in the network in
Eq.~\ref{eq:mixed-network}. }

 
\end{table}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.9\linewidth]{OscialltionsMixedReview}
  \caption{\label{fig:osci}
    Oscillations in the mixed-mechanism network (Eq.~\ref{eq:mixed-network}),
    where rate constants are from Table~\ref{tab:rcs}~\cite{SK},
    and, in contrast with \cite{SK}, 
 
   
    $(K_{\mbox{tot}}, F_{\mbox{tot}},S_{\mbox{tot}})=(100,8,40)$ and the initial values are
    $x_1=0.16309$, $x_2=85.978$,
    $x_3=7.0112$, $x_4=7.0112$, $x_5=1.1$, $x_6=0.13971$,
    $x_7=7.7902$, $x_8=16.854$, $x_9= 0.070112$.
  } 
\end{figure}



\section{Future directions} \label{sec:future}
To add to the future directions already mentioned, here we highlight two more goals:
(1) characterizing bistability (so, going beyond multistationarity) and the corresponding parameter values, and (2) developing criteria for Hopf bifurcations that are well suited to PTM networks.  

Many researchers, ourselves included, have devoted much energy to
analyzing 
multistationarity in PTM networks
(and more generally in networks involving enzymes, substrates, and intermediates~\cite{messi})
and others~\cite{mss-review},
 but what we really are interested in
here, in terms of biology, is bistability.  Indeed, if an
experimentalist discovered a new PTM network and wanted to know what
dynamical properties it has, we have several techniques for assessing
whether this network is multistationary, but embarrassingly few for
determining whether it is bistable \cite{which-small}.  We would also
want to go further, to find witnesses for bistability (values of rate
constants, conservation-law values, and species concentrations at
which the system is bistable) and descriptions of the parameter-space
regions for which the system is bistable (or even just
multistationary).  Progress in this direction has been made for
    distributive dual-site and other PTM networks \cite{fein-024,fein-043}. Such
  developments have been 
 aided by techniques from areas of mathematics such as matroid theory~\cite{KathaMulti,messi}, degree theory~\cite{a6maya, CFMW} and numerical algebraic geometry~\cite{NAG-life-sci,case-study,decomposing}. We expect these and other areas to contribute in the future.

Turning our attention to oscillations, we need more tools for establishing that a network admits sustained oscillations or for precluding them.
Returning to a question posed at the beginning of this Perspective, which PTM networks admit oscillations or at least a Hopf bifurcation?  This area has seen some recent progress~\cite{Errami,jolley,long-term,SK}, but deserves more attention.






\section{Conclusions} \label{sec:end}
What we hoped to convey in this Perspective is that the structure of a PTM network constrains its dynamics (Figure~\ref{fig:multisite-aims}).  
This is an old idea.  The question of how a (general) network constrains its dynamics was the main focus of the mathematical study of reaction systems (including chemical reaction network theory~\cite{FeinOsc}) in the years prior to the scope of this Perspective.  
In other words, the goal was to obtain results without knowing the rate constants.  

An updated version of this goal, 
given our interest in networks coming from biology, 
is to obtain results when not only are the rate constants unknown, but also aspects of the network structure, such as the precise enzymatic mechanism (Figure~\ref{fig:proc-dist}).  
Progress toward this goal has been aided by the methods highlighted in this Perspective --
including using steady-state parametrizations and results that clarify the effect of intermediates or extra reactions.  These techniques form the beginning of a theory of reaction systems well suited to any PTM or other biological signaling networks -- those we already know and also those yet to be discovered.


\subsection*{Acknowledgements} AS was supported by the NSF
  (DMS-1312473/DMS-1513364).
The authors also thank three anonymous referees whose comments helped improve this work.


\bibliographystyle{plain}
