Decomposing the parameter space of biological networks\\ via a numerical discriminant approach
<author name> %% \affil{<number>
<author name> %% \affil{<number>
Roberta Graff\affil{1
Heather A. Harrington
Dhagash Mehta
Helen M. Byrne
Jonathan D. Hauenstein
Many systems in biology, physics and engineering can be described by systems of ordinary differential equation containing many parameters. When studying the dynamic behavior of these large, nonlinear systems, it is useful to identify and characterize the steady-state solutions as the model parameters vary, a technically challenging problem in a high-dimensional parameter landscape. Rather than simply determining the number and stability of steady-states at distinct points in parameter space, we decompose the parameter space into finitely many regions, the steady-state solutions being consistent within each distinct region. From a computational algebraic viewpoint, the boundary of these regions is contained in the discriminant locus. We develop global and local numerical algorithms for constructing the discriminant locus and classifying the parameter landscape. %For low-dimensional parameter spaces, we provide a global picture of the steady state behavior within the parameter landscape. %For high-dimensional parameter spaces, we develop a local picture. We showcase our numerical approaches by applying them to molecular and cell-network models.\\ keywords: parameter landscape, numerical algebraic geometry, discriminant locus, dynamical systems, cellular networks
April9_arXiv_HAH.tex
\documentclass[11pt]{amsart}

\usepackage{amsmath,amsfonts}
\usepackage[table]{xcolor}
\usepackage{graphicx,amssymb,verbatim,epsf,amsmath,setspace,float,dsfont, soul}
\usepackage[numbers, comma, round, sort&compress]{natbib}
%\numberwithin{equation}{section}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\usepackage{tabularx}

\newtheorem{theorem}{Theorem}[section]

\usepackage[font=footnotesize,labelfont=bf]{caption}
\usepackage{anysize}
\marginsize{2.3cm}{2.3cm}{2cm}{3.3cm}

%\usepackage{subcaption}
%\usepackage{mhchem}
%\usepackage[version=3]{mhchem}
\usepackage{multicol}
\makeatletter
%\renewcommand{\thefigure}{S\@arabic\c@figure}

%\usepackage{algorithm}
%\usepackage{algpseudocode}
%\usepackage{algorithmicx}
%
%\renewcommand{\algorithmicrequire}{\textbf{Input:}}
%\renewcommand{\algorithmicensure}{\textbf{Output:}}

%\newcommand\qed{{\hspace*{\fill}$\Box$\vskip12pt plus 1pt}}
%\newenvironment{proof}{{\noindent\bf Proof.\ }}{\qed}

% Turn changes on/off
%\def\tcr{}

%% MACROS
\newcommand{\C}{\mathbb{C}}
\newcommand{\bC}{\mathbb{C}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\bR}{\mathbb{R}}
\newcommand{\PP}{\mathbb{P}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\bfa}{\mathbf{a}}
\newcommand{\cL}{\mathcal{L}}
\newcommand{\cX}{\mathcal{X}}
\newcommand{\cY}{\mathcal{Y}}
\newcommand{\cA}{\mathcal{A}}
\newcommand{\cM}{\mathcal{M}}
\newcommand{\cD}{\mathcal{D}}
\newcommand{\M}{\mathbf{M}}
\newcommand{\V}{\mathcal{V}}
\newcommand{\sP}{\mathcal{P}}


%\begin{document}

%\title{Decomposing parameter space of biological networks using a numerical discriminant locus method}

\title[Decomposing the parameter space of biological networks]{Decomposing the parameter space of biological networks\\ via a numerical discriminant approach}

%% Enter authors via the \author command.  
%% Use \affil to define affiliations.
%% (Leave no spaces between author name and \affil command)

%% Note that the \thanks{} command has been disabled in favor of
%% a generic, reserved space for PNAS publication footnotes.

%% \author{<author name>
%% \affil{<number>}{<Institution>}} One number for each institution.
%% The same number should be used for authors that
%% are affiliated with the same institution, after the first time
%% only the number is needed, ie, \affil{number}{text}, \affil{number}{}
%% Then, before last author ...
%% \and
%% \author{<author name>
%% \affil{<number>}{}}

%% For example, assuming Garcia and Sonnery are both affiliated with
%% Universidad de Murcia:
%% \author{Roberta Graff\affil{1}{University of Cambridge, Cambridge,
%% United Kingdom},
%% Javier de Ruiz Garcia\affil{2}{Universidad de Murcia, Bioquimica y Biologia
%% Molecular, Murcia, Spain}, \and Franklin Sonnery\affil{2}{}}

\author{
Heather A. Harrington, Dhagash Mehta,\\
Helen M. Byrne, Jonathan D. Hauenstein
}
\address{Mathematical Institute, The University of Oxford, Oxford, UK \newline\indent
Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN, USA.}
%%%%%%%%%%%%%%%%%%%
\begin{document}
%% The \maketitle command is necessary to build the title page.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}
Many systems in biology, physics and engineering can be described by systems of ordinary differential equation containing many parameters. When studying the dynamic behavior of these large, nonlinear systems, it is useful to identify and characterize the steady-state solutions as the model parameters vary, a technically challenging problem in a high-dimensional parameter landscape. Rather than simply determining the number and stability of steady-states at distinct points in parameter space, we decompose the parameter space into finitely many regions, the steady-state solutions being consistent within each distinct region.  From a computational algebraic viewpoint, the boundary of these regions 
is contained in the discriminant locus.  
We develop global and local numerical algorithms for 
constructing the discriminant locus and  classifying the parameter landscape. 
%For low-dimensional parameter spaces, we provide a global picture of the steady state behavior within the parameter landscape.
%For high-dimensional parameter spaces, we  develop a local picture.  
We showcase our numerical approaches by applying them to  molecular and cell-network models.\\
keywords: parameter landscape, numerical algebraic geometry, discriminant locus,  dynamical systems, cellular networks
 \end{abstract}

\maketitle
%
%{\color{red}
%\begin{itemize}
%\item should we say more about what we do for high-dimensional parameter spaces? text seems to end prematurely in abstract.[Helen]
%\end{itemize}}

%% When adding keywords, separate each term with a straight line: |


%\section{Significance Statement}
%A key problem in analyzing mathematical models is understanding how the global dynamics change as the system parameters vary. The difficulty of this task is compounded by high dimensional parameter spaces. Many methods rely upon sampling a finite number of points in parameter space and performing local bifurcation analyzes in a neighbourhood of the points; there is no guarantee that such techniques will fully characterize the system dynamics. We present methods that exploit {\em numerical algebraic geometry} to decompose the parameter space into distinct regions that together provide a global view of the steady states and their stability. We showcase our framework on two biological network models.

%{\color{red}
%\begin{itemize}
%\item what do we mean by 'model selection' and 'global best fits?' 
%unsure how this fits into story that we are trying to tell 
%\item what is the optimization problem that we are trying to solve?
%\item where does statistical inference and maximum likelihood appear in this work?
%\item do we identify global optima for the two examples that we study? unclear (to me, anyway) what we are trying to optimise and/or what model selection is being done \item need to make clear (i) that we are restricting attention to systems of ODEs for which the kinetic terms are polynomial in the dependent variables and parameters, and (ii) that we are interested in establishing how the existence and stability of the steady state solutions change within a high-dimensional parameter space [Helen]
%\end{itemize}}

%% Optional for entering abbreviations, separate the abbreviation from
%% its definition with a comma, separate each pair with a semicolon:
%% for example:
%% \abbreviations{SAM, self-assembled monolayer; OTS,
%% octadecyltrichlorosilane}

% \abbreviations{}

%% The first letter of the article should be drop cap: \dropcap{}
%\dropcap{I}n this article we study the evolution of ''almost-sharp'' fronts

%% Enter the text of your article beginning here and ending before
%% \begin{acknowledgements}
%% Section head commands for your reference:
%% \section{}
%% \subsection{}
%% \subsubsection{}


%\textcolor{blue}{para 1. global behavior of dynamical systems, large number of
%parameters and variables, want to understand how parameters affect
%number of steady-states. - Heather}

%{\color{red} 
%\begin{itemize}
%\item [Helen] I think introduction needs to be more bio-physically motivated - why
%do we care about identifying and characterising steady state solutions of a system of ODEs? how will this help us to understand behaviour of a system?
%\item we can probably also trim text
%\item how much detail do we need to include re existing methods for identifying and characterising steady state solutions? make clear that there are lots of methods/approaches but they all have weaknesses - what is so special about discriminant locus method? 
%\end{itemize}
%}
%
\section{Introduction}
The dynamic behavior of many biophysical systems can be mathematically modeled 
with systems of differential equations that describe 
how the state variables interact and evolve over time.  The differential equations 
typically include parameters that represent physical processes such as kinetic rate constants, the strength of cell-cell interactions, and external stimuli, and the behavior of the state variables
may change as the parameters vary.  Typically, determining and classifying all steady state
solutions of such nonlinear systems, as a function of the parameters, is a difficult problem. However, when the equations
are polynomial, or can be translated into polynomials, which is the case for many biological, physics and engineering systems,
computing the steady state solutions 
becomes a problem in computational algebraic geometry.  With this setup, it is possible to compute the regions of the parameter space that give rise to different numbers of steady-state solutions. Here we present an approach 
based on the numerical algebraic geometric paradigm 
of computing points rather than defining equations and inequalities, surpassing current symbolic approaches, with a theoretical guarantee of finding all steady-states.


Due to the ubiquity of such problems, many methods have been proposed
for identifying and characterizing steady state solutions 
over a parameter space.  Some examples include using 
a cylindrical algebraic decomposition \cite{Collins} with
related variants \cite{LR07,xia2007discoverer,brown2003qepcad} and
computing the ideal of the discriminant locus using resultants or Gr\"obner basis methods, 
(see \cite{CLO05,GKZ,Sturmfels}).
Unfortunately, each of these methods
has drawbacks that prevent them from being applicable to models 
for real-world applications due to their algorithmic 
complexity, symbolic expression swell, and inherent sequential
nature.  On the other hand, numerical continuation methods as implemented in, for example,
AUTO \cite{doedel1981auto} and MATCONT \cite{dhooge2003matcont}, 
compute bifurcations for a parametric model and
are local in that one ``continues'' from a given initial point.
Recently, it has become possible to compute all solutions 
over the complex numbers $\bC$ to a system of polynomial equations 
using homotopy continuation and, 
more generally, numerical algebraic geometry,
(see \cite{bates2013numerically,SWbook,wampler1990numerical}).
%%%
%
%{\color{red}
%\begin{itemize}
%\item say what $\bC$ means? similarly $\bR$? [Helen]
%\item are we ok with e.g.? sometimes we say 'for example' - need to decide and be consistent [Helen]
%\end{itemize}}
%
%%%
Such methods have been implemented in software packages
including Bertini \cite{bates2008software}, 
HOM4PS-3 \cite{HOM4PS3}, and PHCpack \cite{verschelde1999algorithm}
with Paramtopy~\cite{Paramotopy} extending Bertini to 
study the solutions at many points in parameter space.
Typically, these methods work over $\bC$ while the solutions
of interest in biological models are in a subset of the real numbers $\bR$,
e.g., one is interested in steady-states in the positive orthant where the variables are~biologically~meaningful.

We present a {\em numerical discriminant locus} method for decomposing parameter space into distinct solution regions. Recall that when solving the roots of the equation $ax^2+bx+c=0$, where $a$, $b$, and $c$ are parameters, and $x$ is the variable, 
the discriminant locus defined by $\Delta := b^2-4ac = 0$ is the 
boundary separating regions in which the two distinct solutions for $x$ are 
real ($\Delta>0$) and nonreal ($\Delta<0$).  We propose three methods
for decomposing the parameter space.  The first, aimed at one-dimensional parameter
spaces, is a sweeping approach that will locate the finitely-many
real points in the discriminant locus building upon~\cite{3Dtumor,PV10}
thereby yielding finitely many regions, which are open intervals in this case, 
where the number of steady-state solutions
is consistent.  The second, for low-dimensional parameter spaces, 
provides a complete decomposition of the parameter space into finitely many regions
after decomposing the discriminant locus.  
Since computing and decomposing the discriminant locus may be impractical
for high-dimensional parameter spaces, our third method 
uses the sweeping approach to compute a local decomposition of the 
parameter space near a given point in the parameter space.  
When decomposing a high-dimensional parameter space is desirable, one could bootstrap together the local analyses to generate a more complete, or global, view of the parameter space.


%%Our sweeping approach evades possible numerical issues from \cite{PV10} 
%%which arise as the number of state variables increase 
%%based on monitoring the determinant of the Jacobian matrix 
%%with respect to the state variables.  That is, the determinant of
%%the Jacobian matrix can be ``close'' to zero for matrices which are 
%%``far'' from being rank deficient.  For example, 
%%consider the matrix $A_n = 2^{-1}\cdot I_n$,
%%where $I_n$ is the $n\times n$ identity matrix.
%%For any $n\geq1$, $A_n$ is a half a unit away from the nearest singular 
%%matrix while $\det(A_n) = 2^{-n}$.  To avoid this situation,
%%one monitors the condition number as in \cite{3Dtumor}.  
%%
%%Our sweeping approach is also different than \cite{3Dtumor} 
%%in that we utilize a perturbation to regularize the sweeping path.
%%That is, rather than track directly along the real line  
%%and possibly pass through a parameter point lying exactly on the 
%%discriminant locus where the Jacobian matrix is rank deficient,
%%we track along a nearby path where the Jacobian matrix has full rank.  
%%This avoids possible issues associated with tracking 
%%where the Jacobian matrix is rank deficient, but still permits
%%the location of the singularities by monitoring the condition number
%%which then could be further analyzed if desired.

%Since our approach is based on the numerical algebraic geometric paradigm 
%of computing points rather than defining equations and inequalities, it surpasses current symbolic approaches, with a theoretical guarantee
%of finding all solutions.


% numerical algebraic
%geometric approach, we demonstrate below 
%that our approach is practical on several examples.  
%
%Several numerical methods have also been proposed.  
%For example, 
%A homotopy sweeping approach was presented in \cite{PV10} 
%for tracking along curves aiming to locate singular solutions.
%This approach tracks along a curve and monitors the determinant of the Jacobian matrix
%to locate singular solutions.  Apart from possibly tracking
%through singularities, this approach can lead to
%unnecessary numerical issues as the size of the problems 
%increase via the determinant.  

%PUT THIS IN ALGORITHM? For example,  $\det(2^{-1}\cdot I_n) = 2^{-n}$ which can be arbitrarily small even though $2^{-1}\cdot I_n$ is far from being singular where $I_n$ is the $n\times n$ identity matrix.


%
%A classical example of a discriminant locus arises from studying the 
%solutions of the parameterized quadratic equation $f(x,a,b,c) = a x ^2 + b x + c = 0$
%where the discriminant is $\Delta = b^2 - 4 a c$.
%Restricting to $a,b,c\in\bR$, the discriminant is the boundary
%of the two regions, namely $\Delta > 0$ where $f = 0$ has two distinct real solutions
%and $\Delta < 0$ where $f = 0$ has two distinct complex solutions.
%Rather than computing explicit equations which define the regions,
%e.g., as in \cite{LR07} which can be computationally infeasible
%for real-world problems, 

In the next section, we introduce the discriminant locus and present
the proposed algorithms. Then we apply our algorithms to two biological models.


\section{Discriminant Locus and the Algorithms}
\label{sec:disc_and_alg}
%%That is, we decompose the space into finitely many regions (nonempty, open, and connected subsets).  In this context, "cells" are just compact regions, but the regions we will compute need not be compact.  Thus, regions are a better word choice.  

We consider autonomous systems of differential equations of the form
\begin{equation} 
\frac{d}{dt} \boldsymbol{x} = \boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) \label{eq:ode} 
\end{equation}
where $ \boldsymbol{x} = (x_1,\dots,x_N)$ 
is the collection of state variables,
$\boldsymbol{p} = (p_1,\dots,p_s)$ is the collection of system parameters,
and $ \boldsymbol{f}(\boldsymbol{x},\boldsymbol{p})$
is a system of $N$ polynomial or rational functions.
For $\boldsymbol{p}\in\bR^s$, the steady state solutions to~Eq.~\ref{eq:ode} 
are defined as $ \boldsymbol{x}\in\bR^N$ such that $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$.  
We are particularly interested in the case when 
$\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$ 
has finitely many solutions in $\bC^N$, all of which are nonsingular, 
for generic~$\boldsymbol{p}$ as this is typically the case for biological networks.

The parameter space $\sP\subset\bR^s$ for Eq.~\ref{eq:ode} 
consists of the values of the parameters $\boldsymbol{p}$ which one
aims to consider, here those that are biologically meaningful.
For simplicity, one should consider $\sP = \bR^s$ 
or the positive orthant in $\bR^s$.  
The quantitative behavior of the steady state solutions 
is constant on regions inside the parameter space $\sP$,
e.g., the number of physically realistic steady state solutions is the same for all parameter values in the region.
One can also refine the quantitative behavior, by considering, for example, the number of positive steady state solutions that are locally stable.
The boundaries of these regions are contained in the {\em discriminant~locus}.

Classically, the discriminant locus consists of the parameter
values for which nongeneric behavior occurs.  
Suppose that $\boldsymbol{p}\in\sP$ is such that $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$ has the generic behavior,
that is, $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$ has the expected number of nonsingular isolated solutions.
By the implicit function theorem, this generic behavior extends to an
open neighborhood containing~$\boldsymbol{p}$.  One can keep increasing the size
of this neighborhood in the parameter space until it touches the discriminant locus.
When focused on counting the number of steady state solutions, 
the discriminant locus is defined by
the closure of all values of $\boldsymbol{p}\in\sP$ such that there 
exists~$\boldsymbol{x}$ with $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$
and $J_{\boldsymbol{x}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p})$ is not invertible, where $J_{\boldsymbol{x}}\boldsymbol{f}$ is the Jacobian matrix of $\boldsymbol{f}$
with respect to the state variables.

One can add to the classical discriminant locus other values of the parameters
based on the quantitative property of interest,
which we collectively call the {\em discriminant locus} for the problem.
For example, the number of positive steady state solutions
occurs at either the classical discriminant locus or the closure
of the set of $\boldsymbol{p}\in\sP$ such that there exists $\boldsymbol{x}$ with 
$\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p}) = 0$ and $x_i = 0$
for some $i\in\{1,\dots,N\}$.  If a function $f_i$ is rational, 
one also has to consider the closure of the set of $\boldsymbol{p}\in\sP$ where
the numerators of $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p})$ are zero and the denominator of $f_i$ is zero.  
Since the discriminant locus is contained in a hypersurface in $\sP$,
the parameter set $\sP$ with the discriminant locus removed consists of regions 
for which the quantitative behavior of the solution set is constant.

\subsection{Algorithm and Implementation}

Our approach builds on methods in numerical algebraic geometry, e.g.,
see \cite{bates2013numerically,SWbook,wampler1990numerical} for a general overview.
In particular, rather than aiming to find equations that vanish
on the discriminant locus and inequalities that describe 
the regions in parameter space, we search for points lying on the discriminant
locus and describe regions based on giving a point lying in the region.  
That is, since each region is connected, one can start with one point 
in a region and trace out the boundaries, which are contained in the discriminant locus, for that region.  In this way, rather than finding explicit
equations and inequalities that define the boundaries of a region,
one starts with a point in the region and has a method that traces
out the boundaries as needed.  This paradigm is similar to the numerical ``cell decomposition'' approach for real 
curves \cite{RealCurves} and surfaces \cite{RealSurfaces1,RealSurfaces2},
rather than a symbolic approach where the boundaries
are described by polynomials \cite{Collins,LR07} (see \cite{niu2008algebraic,hanan2010stability,hernandez2011towards,GHRS15} for applications to biology).

The first approach, which we call {\em perturbed sweeping},
is for one-dimensional parameter spaces, i.e., $s = 1$.  
In this case, the discriminant locus for the problem of interest
consists of at most finitely many points.  Sweeping refers to tracking
the solution curves $\boldsymbol{x}(p)$ where $\boldsymbol{f}(\boldsymbol{x}(p),p) = 0$ as $p\in\bR$ varies via
continuation.  To avoid numerical issues with attempting to track through the
discriminant locus, we propose to sweep along a perturbed path.  
For $i = \sqrt{-1}$ and $\epsilon\in\bR$, we consider the perturbed solution curves
$\boldsymbol{x}_\epsilon(p)$ defined by $\boldsymbol{f}(\boldsymbol{x}_\epsilon(p),p+\epsilon i) = 0$.

\begin{theorem}\label{thm:perturbed}
With the setup described above, for all but finitely many $\epsilon\in\bR$, 
all perturbed solution curves are smooth.
\end{theorem}
\begin{proof}
Since there are only finitely many points in the discriminant locus, 
there can be only finitely many values of $\epsilon\in\bR$ such that there 
exists $\delta\in\bR$ with $\delta + \epsilon i$ in the discriminant locus.
\end{proof}

Since $\boldsymbol{x}_\epsilon(p)\rightarrow \boldsymbol{x}(p)$ as $\epsilon\rightarrow0$, we are able
to recover information about the actual solution curves with the distinct
numerical advantage of tracking {\em smooth} solution curves. 
%By monitoring the condition number, as in \cite{3Dtumor},  as oppose to the determinant in \cite{PV10}, one is able to 
%both robustly locate the points on the discriminant locus using this perturbed sweeping approach.  
%%%%%MOVED THE INTRODUCTION BIT HERE %%%%
Importantly, our sweeping approach evades possible numerical issues from \cite{PV10} 
which arise as the number of state variables increase 
and are associated with monitoring the determinant of the Jacobian matrix with respect to the state variables.  That is, the determinant of
the Jacobian matrix can be ``close'' to zero for matrices which are 
``far'' from being rank deficient.  For example, 
consider the matrix $A_n = 2^{-1}\cdot I_n$,
where $I_n$ is the $n\times n$ identity matrix.
For any $n\geq1$, $A_n$ is half a unit away from the nearest singular 
matrix while $\det(A_n) = 2^{-n}$.  To avoid this situation,
one monitors the condition number as in \cite{3Dtumor}; moreover 
we utilize a perturbation to regularize the sweeping path.
Thus, rather than track directly along the real line  
and possibly pass through a parameter point lying exactly on the 
discriminant locus where the Jacobian matrix is rank deficient,
we track along a nearby path where the Jacobian matrix has full rank.  
This avoids possible issues associated with tracking 
where the Jacobian matrix is rank deficient, but still permits
the location of the singularities by monitoring the condition number.
If further refinement is needed, additional efficient
local computations can be employed, e.g., \cite{GS85,GH15}.

We now build on this perturbed sweeping approach to generate methods
for parameter spaces which are not one-dimensional.  
Our next approach, applicable for low-dimensional parameter spaces, 
is a global approach called a {\em global region decomposition}
that mixes projections, critical sets,
and perturbed sweeping.  We start in the two-dimensional case, $s = 2$,
for which the discriminant locus is contained in a curve.  
For a (sufficiently) random projection $\pi(\boldsymbol{p}) = \alpha_1 p_1 + \alpha_2 p_2$, 
consider intersecting the discriminant locus with the family of lines
defined by $\pi(\boldsymbol{p}) = t$.  Using the perturbed sweeping approach 
on the added parameter~$t$ along the discriminant locus
yields the critical points of the discriminant locus with respect to $\pi$.
The discriminant locus has the same quantitative behavior between
its critical points and thus the regions can be constructed by slicing
between the critical points and at each of the critical points, following
a modification of \cite{RealCurves}.

For a global region decomposition for higher dimensional parameter spaces
one reduces the number of dimension by considering the critical point set. 
For example, in the three-dimensional case, $s = 3$, 
and $\pi(\boldsymbol{p})$ consisting to two (sufficiently) random projections,
one first considers the so-called critical curve of the discriminant locus
with respect to $\pi$, namely the set of points on the discriminant locus
where the Jacobian is singular with respect to $\pi$.  
One then decomposes the curve case as above yielding a region decomposition
in $\pi(\sP)$ which one then lifts to a region decomposition for $\sP$.

Since a global region decomposition is not practical for high-dimensional 
parameter spaces, we propose a local region decomposition method 
by combining perturbed sweeping with the classical approach 
of ray tracing.  Given a point $\boldsymbol{p}\in\sP$ not contained in the discriminant locus,
one can view the codimension one components of the discriminant locus by
using the perturbed sweeping approach along emanating from $\boldsymbol{p}$.  
Once points on the discriminant locus are found, one can use homotopy
continuation methods to track along the boundaries and locate other regions.
This method is local in the sense that small regions could be missed
but avoids the expense of computing critical points of projections
in the global approach.  By bootstrapping local decompositions from 
various $\boldsymbol{p}$, one can aim towards generating a 
global view of the parameter space if desired.

%%FIG1
\begin{figure}[htb!]
\includegraphics{FlowChart_fig1.eps}
\caption{
Flow chart of methods. {\bf(A)} Linear perturbation. {\bf(B)} Global method for decomposing $(a,b)$ parameter space into regions of number of real steady-states. {\bf(C)} Local method for high dimensional parameter space analysis.
}
\label{fig-1}
\end{figure}

\subsection{Quadratic Example}

To illustrate the perturbed sweeping and global region decomposition approach, 
we consider two examples of a parameterized quadratic equation.
The first has one parameter, namely $f(x,a) = x^2 + ax + a/4$.  
The classical discriminant for quadratic polynomials
yields $\Delta = a^2 - a$ with discriminant locus
$a = 0,1$.  In particular, $f = 0$ has 
a singular real solution when $a = 0$ or $a = 1$, 
two distinct real solutions when $a < 0$ or $a > 1$, 
and two distinct nonreal solutions when $0 < a < 1$.  
To avoid tracking through the singularities, 
we use the perturbed sweeping method 
using $f(x,a+\epsilon i) = 0$.  
In this case, for any fixed nonzero $\epsilon\in\bR$ and
$a\in\bR$, $f(x,a+\epsilon i) = 0$ always has 2 distinct solutions.
By taking $\epsilon$ near zero, say $\epsilon = 10^{-6}$, 
we sweep along the smooth curve parameterized by $a$ and 
observe the expected solution behavior as shown in Figure~\ref{fig-1}.  


The second example has two parameters, namely $f(x,a,b) = x^2 + ax + b/4$
and we aim to decompose the parameter space based on the number of real
and positive solutions, which is typical in biological problems.  
The classical discriminant for quadratic polynomials
yields $\Delta = a^2 - b$ with the sign condition adding to this
via the equation $f(0,a,b) = b/4 = 0$, i.e., $b = 0$.  
Thus, the discriminant locus in this case consists of two irreducible curves
which cut the parameter space $(a,b)\in\bR^2$ into four regions 
where the number of real, positive, and negative
solutions are constant on these regions as shown in Figure~\ref{fig-1}.  

By using the projection $\pi(a,b) = a$, the perturbed sweeping method
finds the critical point $a = 0$ of the discriminant locus.
For any $a < 0$ or $a > 0$, there are three regions in $b$, namely
$b < 0$, $0 < b < a^2$, and $b > a^2$.  For $a = 0$, there are two regions
in $b$, namely $b < 0$ and $b > 0$ which must extend away from $a = 0$ by the implicit
function theorem.  Starting from one point in each of these regions, continuation 
is used to merge these $7$ regions into $4$ distinct regions of $\bR^2$.

\section{Results}\label{sec:results}

We showcase our methods by applying them to two illustrative biological models: a detailed ODE model of the gene and protein signaling network that induces long-term memory proposed by Pettigrew {\em et al.}~\cite{Pettigrew:2005ba}, and a network model of cell fate specification in a population of interacting stem cells. Since cellular decision making often depends on the number of accessible (stable) steady-states that a system exhibits, we seek to identify distinct regions of parameter space that can elicit different system behavior.

\subsection{Molecular network model}

A gene and protein network for long term memory was proposed by Pettigrew et al \cite{Pettigrew:2005ba} and investigated using bifurcation and singularity analysis by Song et al. \cite{Song:2006ku}. 
%The long-term facilitation of sensorimotor neuron synapses can be sensitive to actions of serotonin. 
The model developed is of the form of Eq.~\ref{eq:ode} where $\boldsymbol{f}$ 
consists of a total of $15$ polynomial and rational functions,
$\boldsymbol{x} \in \bR^{15}$ is the vector of 15 model variables, 
and $\boldsymbol{p} \in \bR^{40}$ is a vector of 40 parameters. 

\begin{figure}[h!]
\includegraphics{Song_fig2.eps}
\caption{Perturbation sweeping method applied to molecular network of long term memory. $\lambda$ is the bifurcation parameter. Solid lines denote a stable steady-state, dashed lines denote an unstable steady-state. {\bf(A)} Parameter $k_{34} = 0.0022.$ {\bf(B)} Parameter $k_{34} = 0.015.$ {\bf(C)} Parameter $k_{34} = 0.03.$ {\bf(D)} Parameter $k_{34} = 0.1.$
}
\label{fig-bifs}
\end{figure}

The full system $\boldsymbol{f}$ is given in Appendix~\ref{appendix:Song}
for which the equations $\boldsymbol{f}(\boldsymbol{x},\boldsymbol{p})=0$
has $432$ isolated nonsingular solutions for general $\boldsymbol{p}$.
The variable of interest required for long term facilitation is the steady-state of protein kinase A (PKA) in response to the extracellular stimulus parameter $\lambda$. Our aim is to demonstrate the sweeping perturbation method on a large model and reproduce the results in Fig.~5 of \cite{Song:2006ku}. We verify the two stable states and one unstable state, reproducing this figure. However we also find another solution not reported, demonstrating the power of this method (see top branch of Fig.~\ref{fig-bifs}). On inspection, this additional steady state is not on the same branch and is not biologically feasible so we can reject it as nonphysical. From this exhaustive first step of identifying and characterizing all the steady-states, one must exercise caution and systematically check each solution. 
%
%, and changes to the steady-state dynamics as another parameter is changed (for example, irreversible vs reversible . Since the perturbed sweeping method is global in the sense that our method guarantees finding all steady-states, our method found an additional stable steady-state (see solid upper branch in Fig.~\ref{fig:bifs}). Upon further analysis, this stable state is on a separate branch so that a standard continuation methods could not find it and it was not reported by Song et al. \cite{Song:2006ku}. While this method is global in the sense it finds all solutions, one must investigate each steady-state; in this model, the high value of the steady-state is not biologically meaningful in context of the molecular network of long term memory. 

\subsection{Cellular network model}
Most multicellular organisms emerge from a small number of stem-like cells which become increasingly specialized as they proliferate until they transition to one of a finite number of differentiated states \cite{clevers}.  
We propose a caricature model of cell fate specification for a ring of cells, and investigate how cell-cell interactions, mediated by diffusive exchange of a key growth factor, may affect the number of (stable) configurations or patterns that the differentiated cells may adopt.  The model serves as a good test case for discriminant locus methods since, by construction, there is an upper bound on the number of feasible steady states ($2^N$ solutions for a ring of $N$ cells) and some of these patterns are equivalent due to symmetries inherent in the governing equations. 
%
%\begin{itemize}
%\item meaning of above sentences not clear to me
%\item when we talk about bistability, are we referring to stem-like and differentiated states as the two stable steady states? or are we envisaging all cells becoming differentiated, and each cell having two possible fates to choose from?
%\item apologies for the questions - at present it is not clear what we are modelling - we need to be clear in order to interpret results. 
%\item am conscious that we do not account for cell proliferation in model so we need to tread slightly carefully. in view of this, may be more realistic to view steady states as terminal, differentiated states and to assume that each cell could  differentiate into one of two states.
%\item need also to say what $x_i$ represents - concentration of generic growth factor/protein associated with cell $i$, cell-cell communication effected by diffusion of growth factor between neighbouring cells, with cell $i$ communicating with cells $i \pm 1$. 
%\item for example, $x_i$ could represent subcellular levels of notch in intestinal epithelial cells [refs to intestinal crypt - Hans Clevers et c], high and low levels of notch being associated with secretory and absorptive phenotypes respectively 
%\item references:


In more detail, we consider a ring of interacting cells ($i=1, \dots, N$) and denote by $x_i(t) \geq 0$ the concentration within cell $i$ of a growth factor or protein (e.g., notch), whose value
determines that cell's differentiation status \cite{fre,sprinzak,yeung,Visvader:2016fu,Grun:2015ev}. The subcellular dynamics of $x_i$ are represented by a phenomenological function $q(x_i)= - (x-\varepsilon) (a + x) (1 - x)$ with $0 < \varepsilon < a < 1$, this function guaranteeing bistability of each cell in the absence of cell-cell communication. The bistability represents two distinct cell fates; e.g., high and low levels of notch may be associated with differentiation of intestinal epithelial cells into secretory and absorptive phenotypes differentiation \cite{fre,sprinzak,clevers}.
We assume further that cell $i$ communicates with its nearest neighbors (cells $i\pm1$) via diffusive exchange of $x_i$ and denote by parameter $g \geq 0$ the coupling strength. Thus, our cell network model can be written 
%

\begin{equation}
\frac{d x_i}{d t} = q(x_i) + g \cdot \sum_{j= i-1}^{i+1} (x_j - x_i), 
\text{ for all $i=1, \dots, N$.}
\label{eq:cell-network}
\end{equation}
%
\begin{figure}[h!]
\includegraphics[width=\textwidth]{fig3.eps}
\caption{
Global region decomposition applied to coupled cell network model. {\bf(A)} Ring of cell network. Each cell $x_i$ has bistable dynamics at stable states at $0$ and $1$, and unstable state given by $a$. The cells are coupled to neighbor cells by a coupling strength parameter~$g$. {\bf(B)} Region decomposition with parameters $a,g$, similar to a classification diagram, with the number of real stable steady-stages denoted by different colors. The number of stable real steady-states is given for each network where $N = 3,4,5$. {\bf(C)} Steady-state values in state space for $N = 3$ cell network are plotted as parameter $a$ is varied between $[\varepsilon,1]$ or $g$ is varied between 0 and 0.08. The two black dots are always a stable steady-states, the red is always an unstable steady-state, these steady-states are independent of parameter values $a$ and $g$. As the parameter value varies, the color changes according to the color bar.  {\bf(D)} Bifurcation diagram showing the 2-norm of $x$ as parameters are varied  
}
\label{fig-cells}
\end{figure}
We impose periodic boundary conditions so that $x_{N+1} \equiv x_1$ and $x_0 \equiv x_N$, as shown in Fig.~3A.
We analyze the model by applying the global decomposition method for $N=3,4,5$ cells, and construct classification diagrams in $(a,g)$ parameter space (Fig.~3B). We notice that intermediate values of $a$ generate the largest number of real stable steady-state; small and high values of $a$ yield fewer real stable steady-states. Interestingly, all cells synchronize for intermediate to strong values of the coupling parameter $g$ (two stable states in blue region in Fig~\ref{fig-cells}B). We conclude that strong cell-cell communication reduces the number of stable steady state configurations that a population of cells can adopt and, thus, cell-cell communication could be used robustly to drive the cells to a small number of specific states. When coupling is weak ($0 <g<0.1$), the interacting cells have more flexibility in terms of their final states, with the 5-cell network admitting up to 32 stable steady-states (Fig~\ref{fig-cells}B(iii)).  Weaker cell-cell communication allows more patterns to emerge and may be appropriate when it is less important that neighboring cells share the same phenotype. We find that the regions of $(a,g)$ parameter space that give rise to more than two (synchronized) steady states also increase in size as the number of cells increases. 

%\begin{itemize}
%\item {\color{red} not sure about next paragraph - unclear what we are trying to say.
%agree that we need to talk to subplots (C) and (D) but ... think we need to discuss this tomorrow, Heather - sorry!}
%\item {\color{red} need also to include additional references in bibliography}
%\end{itemize}
In addition to decomposing parameter space into regions based on the multiplicity of steady solutions, the method also provides valuable information about how solution structure and stability change as system parameters vary. For example, in Fig.~3C for the $N=3$ cell network, we show how the values and stability of the steady-states for $(x_1,x_2,x_3)$ change as $a$ varies with $g = 0.025$ and as $g$ varies with 
\mbox{$a = 0.4$}.  In Fig.~3D, we plot bifurcation diagrams as $a$ and $g$ vary 
as before for the $N=4$ cell network; instead of presenting particular components $x_i$ ($i=1,2,3,4$), we plot the 2-norm ($\|\boldsymbol{x}\| = (x_1^2 + x_2^2 + x_3^2 + x_4^2)^{1/2}$) to capture the multiplicity of solutions. We note that for $N=3$ and $N=4$ there are always two stable and one unstable steady-states, independent of $(a,g)$ parameter values (as shown by the black and red points in Fig.~3C, and by the solid blue lines in Fig.~3D). We demonstrate the local region decomposition method to analyze a generalization of Eq.~\ref{eq:cell-network} to be cell specific by considering different cell-cell couplings (see Appendix~\ref{appendix:cell}). Through this caricature stem-cell model, we have classified the steady-states, determined which are parameter-independent, explored changes in steady-state behavior as parameters are varied through the distinct regions of the parameter landscape. 

%When n is odd, uniform spacing: 1, 2+n, 2+2n. When n is even, spacing of both n and n/2 increases.

%\begin{itemize}
%\item not sure what more to say about results and/or results for $n$ odd/even
%\item from what I understand, we can conclude that strong cell-cell communication reduces number of stable steady state configurations that cells can adopt and so cell-cell communication possible mechanism for generating robust distributions of differentiated cells; weaker cell-cell communication allows more patterns to emerge and may be appropriate when less important that neighboring cells share same phenotype.  
%\end{itemize}

%Helen write things about intuition



\section{Conclusion}

We have presented a suite of new methods, based on computational algebraic geometry, for decomposing the parameter space associated with a dynamical system into distinct regions based on the multiplicity and stability of its steady-state solutions. The methods enable us to understand the parameter landscape of high-dimensional, ordinary differential models with large numbers of parameters. 
These methods have considerable potential: they could be used to analyze differential equation models associated with a wide range of real-world problems in biology, science and engineering which cannot be tackled with~existing~approaches.

We have demonstrated that groups of cells, especially cells that can have bistable internal dynamics, when coupled, increase the number of real stable steady-states. First we have considered all the cells to be homogeneous, and even with this simple model, we can gain new insight into how stem cells may transit into fully differentiated cells, with much richer dynamics, just by considering the interaction of the cells as a network. We considered different couplings between cells using a local region decomposition in the Appendix.  We remark that it would be interesting to understand how a mutation in a single cell (for example, set $a$ very small) compares to real biological systems. In the future, it would be 
beneficial to consider different networks of 
cells (e.g., two growth factors, generalizations to 2D lattices) 
and the dynamical regimes these systems can admit. 
%\medskip
%HEATHER/HELEN: I DO NOT UNDERSTAND "increase the number of real stable steady-states".
%As the coupling $g$ increases, the number of stable steady-states can only only decrease ($g = 0$: all $2^N$ are real stable steady-states)
%\medskip

%In the future, it would be fruitful to explore singularly perturbed problems and consider how this approach relates to the rich literature of symmetries in bifurcations \cite{golubitsky}.  
%
%%%%
%{\color{red} not sure what this means. need to relate to examples and/or provide more context.
%think conclusions are a bit repetitive - need to streamline them - what are the 3-4 key points?
%new method(s), wide applicability, ... [Helen]}
%%%%

By embracing and exploiting the complexity naturally present in 
biology, which often translates to nonlinearity in biological models, there is an incredible opportunity to 
develop and apply numerical algebraic geometric methods to study 
such models.  With nonlinearity arising in all areas of science and engineering,
the increasing complexity of various biological models provides a fertile 
testing ground for nonlinear approaches with a long-term goal
of making numerical algebraic geometric methods  
as ubiquitous as methods in numerical linear algebra.  





\section*{Acknowledgement}
We thank J. Byrne and G. Moroz for for helpful discussions.
JDH was supported in part by NSF ACI 1460032, Sloan Research Fellowship, and Army Young Investigator Program (YIP). HAH acknowledges funding from AMS Simons Travel Grant and EPSRC Fellowship EP/K041096/1.

\bibliographystyle{pnas}
\bibliography{bibliofile}
%\newpage
\appendix

\section{The Song et al. model}\label{appendix:Song}

For completeness, we reproduce below the system of 15 ordinary differential equations proposed by Song et al. \cite{Song:2006ku}. The parameters we vary are $\lambda=[5-HT]$ and $k_{34}$ = $k_{\mbox{ApSyn}}$. All other parameters are fixed to values given in \cite{Song:2006ku}.

{\scriptsize
\begin{eqnarray}
\frac{d[cAMP]}{dt} &=& 3.6 \left(\frac{[5-HT]}{[5-HT]+14}\right) - ([cAMP] - 0.06)\\
\frac{d[R]}{dt} &=& V_{\mbox{syn}} + k_{\mbox{fpka}}[RC] [cAPM]^2 - k_{\mbox{bpka}}[R][C]\nonumber \\
& & - k_{\mbox{Ap-Uch}}[R]([\mbox{Ap-Uch}]_{\tau} - \mbox{Ap-Uch}_{\mbox{basal}}) - k_{\mbox{dpka}}[R] \\
\frac{d [C]}{dt} &=& V_{\mbox{syn}} + k_{\mbox{fpka}}[RC] [cAPM]^2 - k_{\mbox{bpka}}[R][C]\nonumber \\
& &+ k_{\mbox{Ap-Uch}}[RC]([\mbox{Ap-Uch}]_{\tau} - \mbox{Ap-Uch}_{\mbox{basal}}) - k_{\mbox{dpka}}[C]\\
\frac{d [RC]}{d t} &=& k_{\mbox{bpka}} [R][C] - k_{fpka}[RC] [cAMP]^2 \nonumber \\
& & - k_{\mbox{Ap-Uch}}[RC]([\mbox{Ap-Uch}]_{\tau} - \mbox{Ap-Uch}_{\mbox{basal}}) - k_{\mbox{dpka}}[RC]\\
\frac{d [REG]}{d t} &=& k_{\mbox{translation}}[\mbox{mRNA}_{\mbox{REG}}] ERK_{\mbox{act}} \nonumber \\ 
& & - v_{\mbox{dreg}}\left(\frac{[REG]}{[REG] + K_{\mbox{dreg}}} \right) - k_{\mbox{dsm}}[REG] \\ 
\frac{d [pREG]}{dt} &=& v_{\mbox{rphos}}ERK_{\mbox{act}} \left(\frac{[REG] - [pREG]}{[REG]-[pREG] + K_{\mbox{rphos}}} \right) - v_{\mbox{dreg}} \left(\frac{[pREG]}{[REG] + K_{\mbox{dreg}}} \right) \nonumber \\
& & - k_{\mbox{dsm}}[\mbox{pREG}] \\
\frac{d [\mbox{mRNA}_{\mbox{REG}}]}{d t} &=& v_{\mbox{mREG}} - v_{\mbox{dmreg}}[5-HT] \left(\frac{[\mbox{mRNA}_{\mbox{REG}}]}{[\mbox{mRNA}_{\mbox{REG}}] + K_{\mbox{dmreg}}}\right) - k_{\mbox{dmreg}} [\mbox{mRNA}_{\mbox{REG}}] \\
\frac{d [\mbox{Raf}]}{d t} &=& - k_{f, \mbox{Raf}} [5-HT] [Raf] + k_{b,\mbox{Raf}} [\mbox{Raf}^{p}] \\
\frac{d [\mbox{MAPKK}]}{d t} &=& -k_{f, \mbox{MAPKK}}[\mbox{Raf}^{p}] \left(\frac{[\mbox{MAPKK}]}{[\mbox{MAPKK}] + K_{\mbox{MK}}} \right) + k_{b,\mbox{MAPKK}}\left(\frac{[\mbox{MAPKK}^{p}]}{[\mbox{MAPKK}^{p}] + K_{\mbox{MK}}} \right)\\
\frac{d[MAPKK]^{pp}}{d t} &=& k_{f,\mbox{MAPKK}}[\mbox{Raf}^{p}]\left(\frac{[\mbox{MAPKK}^{p}]}{[\mbox{MAPKK}^{p}] + K_{\mbox{MK}}} \right)
- k_{b,\mbox{MAPKK}}\left(\frac{[\mbox{MAPKK}^{pp}]}{[\mbox{MAPKK}^{pp}] + K_{\mbox{MK}}} \right) \\
\frac{d [\mbox{ERK}]}{d t} &=& - k_{f,\mbox{ERK}} [\mbox{MAPKK}^{pp}] \left( \frac{[\mbox{ERK}]}{[\mbox{ERK}] + K_{\mbox{MK}}} \right)
+ k_{b, \mbox{ERK}} \left( \frac{[\mbox{ERK}^{p}]}{[\mbox{ERK}^{p}] + K_{\mbox{MK}}}\right)\\
\frac{d [\mbox{ERK}^{pp}]}{d t} &=&  k_{f,\mbox{ERK}} [\mbox{MAPKK}^{pp}] \left( \frac{[\mbox{ERK}^{p}]}{[\mbox{ERK}^{p}] + K_{\mbox{MK}}} \right)
- k_{b, \mbox{ERK}} \left( \frac{[\mbox{ERK}^{pp}]}{[\mbox{ERK}^{pp}] + K_{\mbox{MK}}}\right) \\
\frac{d P_{\mbox{pka}}}{d t} &=& k_{\mbox{phos1}} \mbox{PKA}_{\mbox{act}} (1 - P_{\mbox{pka}}) - k_{\mbox{dephos1}} \mbox{PPhos} P_{\mbox{pka}}\\
\frac{d P_{\mbox{erk}}}{d t} &=& k_{\mbox{phos2}} \mbox{ERK}_{\mbox{act}} (1 - P_{\mbox{erk}}) - k_{\mbox{dephos2}} \mbox{PPhos} P_{\mbox{erk}}\\
\frac{d [\mbox{AP-Uch}]}{d t} &=& k_{\mbox{ApSyn}}\left(\frac{P^{2}_{\mbox{pka}}}{P^{2}_{\mbox{pka}} + K^{2}_{\mbox{pka}}} \right)
\left(\frac{P^{2}_{\mbox{erk}}}{P^{2}_{\mbox{erk}} + K^{2}_{\mbox{erk}}} \right) + k_{\mbox{ApSynBasal}} - k_{\mbox{deg}}[\mbox{Ap-Uch}]
\end{eqnarray}
}
where
{\scriptsize
\begin{eqnarray*}
 PKA_{\mbox{act}} &=& [C], \\ 
 \left[\mbox{Raf}^{p} \right] &=& [ \mbox{Raf}_{\mbox{tot}} ] - [ \mbox{Raf} ], \\
 k_{\mbox{bpka}}  &=& \left(\frac{k_{\mbox{baspka}}}{1+[pREG]/K_{\mbox{reg}}} \right), \\ 
 \left[\mbox{ERK}^{p}\right] &=& [\mbox{ERK}_{\mbox{tot}}] - [\mbox{ERK}] - [\mbox{ERK}^{pp}],  \\
 \mbox{ERK}_{\mbox{act}} &=& [\mbox{ERK}^{pp}] + \mbox{ERK}_{\mbox{basal}}, \\ 
\left[\mbox{MAPKK}^{p}\right] &=& [\mbox{MAPKK}_{\mbox{tot}}] - [\mbox{MAPKK}] - [\mbox{MAPKK}^{pp}]. 
\end{eqnarray*}
}
Following \cite{Song:2006ku}, we fix the model parameters at the following values:
{\scriptsize
\begin{eqnarray*}
 V_{\mbox{syn}} &=& 0.002 \mbox{min}^{-1},
 k_{\mbox{fpka}} = 105.0 \mbox{min}^{-1},\nonumber \\
 k_{\mbox{Ap-Uch}} &=& 0.007 \mbox{min}^{-1}, 
 \mbox{Ap-Uch}_{\mbox{basal}} = 0.10, \nonumber \\
 k_{\mbox{dpka}} &=& 0.00048 \mbox{min}^{-1}, 
 \tau = 250 \mbox{min}, \nonumber \\
 k_{\mbox{baspka}} &=& 12 \mbox{min}^{-1}, 
 k_{\mbox{reg}} = 0.00064,\nonumber \\
 k_{\mbox{translation}} &=& 4 \mbox{min}^{-1},
 v_{\mbox{rphos}} = 1 \mbox{min}^{-1}, K_{\mbox{rphos}} = 1.5, \nonumber \\
 v_{\mbox{dreg}} &=& 0.16 \mbox{min}^{-1}, 
 K_{\mbox{dreg}} = 0.0015, \nonumber \\
 k_{\mbox{dsm}} &=& 0.02 \mbox{min}^{-1}, 
 v_{\mbox{mREG}} = 2\time 10^{-5} \mbox{min}^{-1}, \nonumber \\
 v_{\mbox{dmreg}} &=& 0.00225 \mbox{min}^{-1}, 
 K_{\mbox{dmreg}} = 0.01,\nonumber \\
 k_{\mbox{dmreg}} &=& 3 \times 10^{-5} \mbox{min}^{-1}, 
 k_{f, \mbox{Raf}} = 0.00345 \mbox{min}^{-1}, \nonumber \\
 k_{f, \mbox{MAPKK}} &=& 0.7 \mbox{min}^{-1}, 
 k_{f, \mbox{ERK}} = 0.44 \mbox{min}^{-1}, \nonumber \\
 k_{b, \mbox{Raf}} &=& 0.001 \mbox{min}^{-1}, 
 k_{b, \mbox{MAPKK}} = 0.12 \mbox{min}^{-1}, \nonumber \\
 k_{b, \mbox{ERK}} &=& 0.12 \mbox{min}^{-1},
  [\mbox{Raf}_{\mbox{tot}}] = 0.5, \nonumber \\
  \mbox{MAPKK}_{\mbox{tot}} &=& 0.5, 
[ERK_{\mbox{tot}}] = 0.5 \nonumber \\
K_{\mbox{MK}} &=& 0.08, 
\mbox{ERK}_{\mbox{basal}} = 0.015, \nonumber \\
k_{\mbox{phos1}} &=& 0.01 \mbox{min}^{-1}, 
k_{\mbox{dephos1}} = 1.5 \mbox{min}^{-1}, \nonumber \\
\mbox{PPhos} &=& 0.1, 
k_{\mbox{phos2}} = 0.005 \mbox{min}^{-1},\nonumber \\
k_{\mbox{ApSyn}} &=& 0.02 \mbox{min}^{-1}, 
k_{\mbox{ApSynBasal}} = 0.0009 \mbox{min}^{-1}, \nonumber \\
K_{\mbox{pka}} &=& 0.2,
K_{\mbox{erk}} = 0.004, \nonumber \\
k_{\mbox{deg}} &=& 0.01 \mbox{min}^{-1}
\end{eqnarray*}
}
%\baselineskip24pt
We remark that the discriminant method can handle rational functions.
For this example, the denominators do not vanish near the regions of 
interest so they do not have any impact on the behavior of the solutions.
If the denominator also vanished when finding a solution to the system of equations from the numerators, then the parameter values for which this occurs would be added into the discriminant~locus.

\section{Cell network model}\label{appendix:cell}

The model in Eq.~\ref{eq:cell-network} used uniform coupling $g$ for
all nearest neighbors.  Using \cite{HS10}, 
as a set, the classical discriminant locus for this model 
with $N = 4$ has degree $72$.
That is, there is a degree $72$ polynomial $\Delta(a,g)$ such
that the classical discriminant locus is defined by $\Delta = 0$.
With $a = 0.4$ as in Fig~3D, the univariate polynomial equation
$\Delta(0.4,g) = 0$ has $45$ complex solutions, $25$ of which are real with $15$ positive with only~$4$ corresponding to a change in the 
number of stable steady-state solutions yielding the 
regions (intervals) for $g$ approximately:
$${\small
\begin{array}{r|ccccc} 
\mbox{regions} & [0,0.0197) & (0.0197,0.0206) & (0.0206,0.0411) & (0.0411,0.0533) & (0.0533,\infty) \\
\hline 
\mbox{\# stable steady-state solns} & 16 & 12 & 10 & 6 & 2 
\end{array}
}$$

We now consider a generalization of this model
which uses coupling strengths $g_{i,i+1}=g_{i+1,i}\geq0$ 
between cell $i$ and $i+1$ with cyclic 
ordering ($N+1\equiv 1$), namely
\begin{equation}
\frac{d x_i}{d t} = q(x_i) + \sum_{j= i-1}^{i+1} g_{i,j} \cdot (x_j - x_i), 
\text{ for all $i=1, \dots, N$.}
\label{eq:cell-network2}
\end{equation}
With $N = 4$, as a set, the classical discriminant locus for this generalized 
model has degree $486$.  We consider the case with $a = 0.4$ and 
$g_{4,1} = g_{1,4} = 0$ leaving three free 
parameters $g_{1,2}$, $g_{2,3}$, and $g_{3,4}$.  

The perturbed sweeping approach along the ray defined
by $g_{i,i+1} = i\cdot t$ for $i = 1,2,3$ and $t\geq0$
decomposes the space $t\geq0$ into $15$ regions, approximately
$$
\begin{array}{c|c} 
\mbox{regions} & \mbox{\# stable steady-state solns} \\
\hline
[0,0.0079) & 16 \\
(0.0079,0.0080)& 15 \\
(0.0080,0.0132)& 14 \\
(0.0132,0.0136)& 13 \\
(0.0136,0.0137) & 12 \\
(0.0137,0.0142)& 11 \\
(0.0142,0.0153)& 10 \\
(0.0153,0.0161)& 9 \\
(0.0161,0.0171)& 8 \\
(0.0171,0.0237)& 7 \\
(0.0237,0.0352& 6 \\
(0.0352,0.0360)& 5 \\
(0.0360,0.0407)& 4 \\
(0.0407,0.1264)& 3 \\
(0.1264,\infty)& 2 
\end{array}
$$
As a comparison, the classical discriminant with respect to $t\in\bC$ 
consists of $312$ distinct points of which $84$ are real with $42$ positive
with $14$ corresponding to a change in the number stable steady-state solutions.
Figure~\ref{fig-sweep} plots the region with $16$ 
stable steady-state solutions along various rays emanating from the origin.

\begin{figure}[h!]
\includegraphics[scale=0.4]{Sweeping.eps}
\caption{Plot of the regions with $16$ stable steady-state solutions
along various rays emanating from $g_{1,2} = g_{2,3} = g_{3,4} = 0$.}
\label{fig-sweep}
\end{figure}

\end{document}


April9_arXiv_HAH.bbl
\def\cprime{$'$}
\begin{thebibliography}{10}

\bibitem{Collins}
Collins GE
\newblock (1975) in \emph{Automata theory and formal languages ({S}econd {GI}
  {C}onf., {K}aiserslautern, 1975)}
\newblock (Springer, Berlin), pp 134--183. Lecture Notes in Comput. Sci., Vol.
  33.

\bibitem{LR07}
Lazard D, Rouillier F
\newblock (2007) Solving parametric polynomial systems.
\newblock \emph{Journal of Symbolic Computation} 42:636 -- 667.

\bibitem{xia2007discoverer}
Xia B
\newblock (2007) Discoverer: a tool for solving semi-algebraic systems.
\newblock \emph{ACM Communications in Computer Algebra} 41:102--103.

\bibitem{brown2003qepcad}
Brown CW
\newblock (2003) Qepcad b: a program for computing with semi-algebraic sets
  using cads.
\newblock \emph{ACM SIGSAM Bulletin} 37:97--108.

\bibitem{CLO05}
Cox DA, Little J, O'Shea D
\newblock (2005) \emph{Using algebraic geometry}, Graduate Texts in Mathematics
\newblock (Springer, New York) Vol.{} 185, Second edition.

\bibitem{GKZ}
Gel{\cprime}fand IM, Kapranov MM, Zelevinsky AV
\newblock (1994) \emph{Discriminants, resultants, and multidimensional
  determinants}, Mathematics: Theory \& Applications
\newblock (Birkh\"auser Boston, Inc., Boston, MA).

\bibitem{Sturmfels}
Sturmfels B
\newblock (2002) \emph{Solving systems of polynomial equations}, CBMS Regional
  Conference Series in Mathematics
\newblock (Published for the Conference Board of the Mathematical Sciences,
  Washington, DC; by the American Mathematical Society, Providence, RI)
  Vol.{}~97.

\bibitem{doedel1981auto}
Doedel EJ
\newblock (1981) Auto: A program for the automatic bifurcation analysis of
  autonomous systems.
\newblock \emph{Congr. Numer} 30:265--284.

\bibitem{dhooge2003matcont}
Dhooge A, Govaerts W, Kuznetsov YA
\newblock (2003) Matcont: a matlab package for numerical bifurcation analysis
  of odes.
\newblock \emph{ACM Transactions on Mathematical Software (TOMS)} 29:141--164.

\bibitem{bates2013numerically}
Bates DJ, Hauenstein JD, Sommese AJ, Wampler CW
\newblock (2013) \emph{Numerically solving polynomial systems with Bertini}
\newblock (SIAM) Vol.{}~25.

\bibitem{SWbook}
Sommese AJ, Wampler, II CW
\newblock (2005) \emph{The numerical solution of systems of polynomials}
\newblock (World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ)
\newblock Arising in engineering and science.

\bibitem{wampler1990numerical}
Wampler C, Sommese A, Morgan A
\newblock (1990) Numerical continuation methods for solving polynomial systems
  arising in kinematics.
\newblock \emph{Journal of mechanical design} 112:59--68.

\bibitem{bates2008software}
Bates DJ, Hauenstein JD, Sommese AJ, Wampler~II CW
\newblock (2008) in \emph{Software for algebraic geometry}
\newblock (Springer), pp 1--14.

\bibitem{HOM4PS3}
Chen T, Lee TL, Li TY
\newblock (2014) in \emph{Mathematical software---{ICMS} 2014}, Lecture Notes
  in Comput. Sci.
\newblock (Springer, Heidelberg) Vol.{} 8592, pp 183--190.

\bibitem{verschelde1999algorithm}
Verschelde J
\newblock (1999) Algorithm 795: Phcpack: A general-purpose solver for
  polynomial systems by homotopy continuation.
\newblock \emph{ACM Transactions on Mathematical Software (TOMS)} 25:251--276.

\bibitem{Paramotopy}
Bates DJ, Brake DA, Niemerg ME
\newblock (2015) Paramotopy: Parameter homotopies in parallel. (Available at
  {\tt paramotopy.com}).

\bibitem{3Dtumor}
Hao W, Hauenstein JD, Hu B, Sommese AJ
\newblock (2011) A three-dimensional steady-state tumor system.
\newblock \emph{Appl. Math. Comput.} 218:2661--2669.

\bibitem{PV10}
Piret K, Verschelde J
\newblock (2010) Sweeping algebraic curves for singular solutions.
\newblock \emph{J. Comput. Appl. Math.} 234:1228--1237.

\bibitem{RealCurves}
Lu Y, Bates DJ, Sommese AJ, Wampler CW
\newblock (2007) in \emph{Algebra, geometry and their interactions}, Contemp.
  Math.
\newblock (Amer. Math. Soc., Providence, RI) Vol.{} 448, pp 183--205.

\bibitem{RealSurfaces1}
Bates DJ, Brake DA, Hauenstein JD, Sommese AJ, Wampler CW
\newblock (2014) in \emph{Mathematical software---{ICMS} 2014}, Lecture Notes
  in Comput. Sci.
\newblock (Springer, Heidelberg) Vol.{} 8592, pp 246--252.

\bibitem{RealSurfaces2}
Besana GM, Di~Rocco S, Hauenstein JD, Sommese AJ, Wampler CW
\newblock (2013) Cell decomposition of almost smooth real algebraic surfaces.
\newblock \emph{Numer. Algorithms} 63:645--678.

\bibitem{niu2008algebraic}
Niu W, Wang D
\newblock (2008) Algebraic approaches to stability analysis of biological
  systems.
\newblock \emph{Mathematics in Computer Science} 1:507--539.

\bibitem{hanan2010stability}
Hanan W, Mehta D, Moroz G, Pouryahya S
\newblock (2010) Stability and bifurcation analysis of coupled fitzhugh-nagumo
  oscillators.
\newblock \emph{"Extended abstract" published in the Joint Conference of
  ASCM2009 and MACIS2009, Japan, 2009. arXiv:1001.5420}.

\bibitem{hernandez2011towards}
Hernandez-Vargas EA, Mehta D, Middleton RH
\newblock (2011) Towards modeling hiv long term behavior.
\newblock \emph{Proceedings of 18th IFAC World Congress, Milan, 2011.
  arXiv:1105.2823}.

\bibitem{GHRS15}
Gross E, Harrington HA, Rosen Z, Sturmfels B
\newblock (2016) Algebraic systems biology: A case study for the wnt pathway.
\newblock \emph{Bulletin of Mathematical Biology} 78:21--51.

\bibitem{GS85}
Golubitsky M, Schaeffer DG
\newblock (1985) \emph{Singularities and groups in bifurcation theory. {V}ol.
  {I}}, Applied Mathematical Sciences
\newblock (Springer-Verlag, New York) Vol.{}~51, pp xvii+463.

\bibitem{GH15}
Griffin ZA, Hauenstein JD
\newblock (2015) Real solutions to systems of polynomial equations and
  parameter continuation.
\newblock \emph{Adv. Geom.} 15:173--187.

\bibitem{Pettigrew:2005ba}
Pettigrew DB, Smolen P, Baxter DA, Byrne JH
\newblock (2005) Dynamic properties of regulatory motifs associated with
  induction of three temporal domains of memory in aplysia.
\newblock \emph{Journal of computational neuroscience} 18:163--181.

\bibitem{Song:2006ku}
Song H, Smolen P, Av-Ron E, Baxter DA, Byrne JH
\newblock (2006) Bifurcation and singularity analysis of a molecular network
  for the induction of long-term memory.
\newblock \emph{Biophysical journal} 90:2309--2325.

\bibitem{clevers}
Clevers H
\newblock (2015) {STEM CELLS. What is an adult stem cell?}
\newblock \emph{Science} 350:1319--1320.

\bibitem{fre}
Fre S, {et~al.}
\newblock (2005) {Notch signals control the fate of immature progenitor cells
  in the intestine.}
\newblock \emph{Nature} 435:964--968.

\bibitem{sprinzak}
Sprinzak D, {et~al.}
\newblock (2010) {Cis-interactions between Notch and Delta generate mutually
  exclusive signalling states.}
\newblock \emph{Nature} 465:86--90.

\bibitem{yeung}
Yeung TM, Chia LA, Kosinski CM, Kuo CJ
\newblock (2011) {Regulation of self-renewal and differentiation by the
  intestinal stem cell niche.}
\newblock \emph{Cellular and molecular life sciences : CMLS} 68:2513--2523.

\bibitem{Visvader:2016fu}
Visvader JE, Clevers H
\newblock (2016) {Tissue-specific designs of stem cell hierarchies.}
\newblock \emph{Nat Cell Biol} 18:349--355.

\bibitem{Grun:2015ev}
Gr{\"u}n D, {et~al.}
\newblock (2015) {Single-cell messenger RNA sequencing reveals rare intestinal
  cell types.}
\newblock \emph{Nature} 525:251--255.

\bibitem{HS10}
Hauenstein JD, Sommese AJ
\newblock (2010) Witness sets of projections.
\newblock \emph{Appl. Math. Comput.} 217:3349--3354.

\end{thebibliography}
