Numerical algebraic geometry for model selection and its application to the life sciences

Researchers working with mathematical models are often confronted by the related problems of parameter estimation, model validation and model selection. These are all optimization problems, well known to be challenging due to nonlinearity, non-convexity and multiple local optima. Furthermore, the challenges are compounded when only partial data are available. Here, we consider polynomial models (e.g. mass-action chemical reaction networks at steady state) and describe a framework for their analysis based on optimization using numerical algebraic geometry. Specifically, we use probability-one polynomial homotopy continuation methods to compute all critical points of the objective function, then filter to recover the global optima. Our approach exploits the geometrical structures relating models and data, and we demonstrate its utility on examples from cell signalling, synthetic biology and epidemiology.


Introduction
Across the physical, biological, and social sciences, mathematical models are formulated and studied to better understand real-world phenomena.Often, multiple models are developed to explore alternate hypotheses.It then becomes necessary to choose between different models, for example, based on their fit with noisy experimental data.This is the problem of model selection, a fundamental scientific problem with practical implications (1)(2)(3).
The standard approach to model selection is to first estimate all model parameters and hidden variables from the data, then select the model with the smallest best-fit error, up to some penalty on model complexity (4,5).Thus, at its core, model selection is intimately tied to parameter estimation, which, for a given model f (a, x) = 0 in the parameters a and variables x with measurable "outputs" y = g(x), say, can be written as the following optimization problem: where y denotes the observed data, i.e. measured outputs.Unless f and g are convex, solving (1) is a nonconvex problem, which can be challenging as standard local solvers run the risk of getting trapped in local minima (especially in high dimensions).This can be mitigated somewhat with techniques such as simulated annealing (6,7) or convex relaxation which has been successful for model invalidation (8)(9)(10), but there is generally no guarantee that a global minimum will be found.When f and g are polynomial, however, problem (1) can be solved globally by finding all roots of an associated polynomial system.In this case, ideas from computational algebra and algebraic geometry can be effective; see, e.g., (11)(12)(13)(14) for applications of Gröbner bases in systems biology and (15) for applications of algebraic geometry to statistical inference.Such symbolic methods tend to be computationally expensive, which limits their use in practice.Thus, although they provide a solution in principle, new algorithms and techniques are yet desired.
In this paper, we aim to fill this gap by proposing a framework for global parameter estimation for polynomial deterministic models using numerical algebraic geometry (NAG), a suite of tools for numerically approximating the solution sets of multivariate polynomial systems via adaptive multi-precision, probabilityone polynomial homotopy continuation (16,17).Our approach scales well in dimension relative to classical symbolic methods (18), and, while it comes with a higher computational cost than standard local optimization, it has a probability-one guarantee to recover the global optima, solving problem (1) in the strong sense.This allows us to reason rigorously about model fit and to address the related problems of model selection and parameter estimation from a maximum-likelihood perspective.
We demonstrate our techniques on examples from biology, where polynomial models often arise as the steady-state descriptions of mass-action chemical reaction networks.Although some limitations remain, we believe that this work achieves its primary purpose of introducing NAG as a valuable complement to existing tools for model evaluation and analysis.Additionally, this paper highlights specific challenges that arise when using polynomial methods for model inference, such as dealing with positivity constraints and non-isolated solutions, and provides guidance for tackling these challenges.
The remainder of the paper is organized as follows.In the next section, we state precisely the problems with which we are concerned: model validation, model selection, and parameter and hidden variable estimation.We then present NAG algorithms for solving each problem.Finally, we showcase our approach on a few examples including cell death activation, synthetic bio-circuits, HIV progression, and protein modification.

Problem Statement
Consider a model whose dynamics are described by the system of polynomial differential equations where a = (a 1 , . . ., a k ) are parameters (e.g.rate constants in a deterministic mechanistic model, such as a chemical reaction network with mass-action kinetics), x = (x 1 , . . ., x n ) are variables, and f = (f 1 , . . ., f r ) are polynomials in x and a with measurable outputs y = g(x) where y = (y 1 , . . ., y m ), m ≤ n, and g = (g 1 , . . ., g m ) are polynomials in x.While the parameters a 1 , . . ., a k are treated as fixed variables in our exposition, we separate them from x 1 , . . ., x n to respect how such variables are treated differently in experimental and computational settings.
In algebraic geometry a variety is a solution set of a system of polynomial equations; we use this terminology for our next two definitions.The real model variety is the solution set of the system f = 0 (3) corresponding to the steady states of the model.Now, consider, for simplicity, the case of a single data point ŷ = (ŷ 1 , . . ., ŷm ).(See SI Appendix for multiple data points.)The real data variety is then the affine linear space (V D ) R := {(a, x, y) ∈ R k+n+m : y i = ŷi , ∀i = 1, . . ., m}, with dim(V D ) R = k + n.We consider the case when the data includes some extrinsic (measurement) noise; we assume the errors { 1 , . . ., m } on the observed data variables are uncorrelated random variables and each error i is normally distributed with known variance σ i (which can be obtained by instrument calibration).Using this geometric framework, the problems of (1) model validation, (2) model selection, and (3) parameter estimation, can be described precisely in terms of the real algebraic varieties (V M ) R and (V D ) R .

Problem 1: Model Validation
For model validation, we want to determine whether a deterministic polynomial model M is compatible with the data according to a given significance level α.This is akin to asking whether the model is a "good fit" for the data.A natural goodness-of-fit statistic is: When the variances differ, d 2 is the minimum squared weighted Euclidean distance between (V M ) R and (V D ) R .When all variances are equal to one, the value of ( 5) is just the minimum squared distance.For the remainder of the text, we will assume the latter, knowing that we can simply rescale.
Optimization problem 5 can be derived directly from the log-likelihood function log L(a, x, y) as demonstrated in Section 2.1 of the SI.This provides a bridge to standard statistical model selection tools such as AIC (19) and BIC (20).
If the data are generated from the model M with normally distributed extrinsic noise, the distribution function of d 2 is dominated by that of the chi-squared distribution with m degrees of freedom, χ 2 m , where m is the number of measurable outputs.Thus, if p α is the upper α-percentile for χ 2 m , then Pr(d 2 ≥ p α ) ≤ Pr(U ≥ p α ) = α where U ∼ χ 2 m .We reject the model M as incompatible if d 2 > p α ; otherwise we say that the model M is compatible.
If the real model and data varieties intersect, that is, (V M ) R ∩ (V D ) R = ∅, then d 2 = 0, and we also say that the model is compatible with the data.If there are restrictions on (a, x, y), for example if all parameters and variables are required to be non-negative, then finding d 2 becomes a constrained optimization problem (see SI Appendix).

Problem 2: Model Selection
For model selection, we are given a set of models {M 1 , . . ., M s } and want to determine the model of best fit.In this setting, we use the statistic d 2 to make a selection, either by choosing the model with the smallest goodness-of-fit statistic or by using d 2 in conjunction with a complexity penalty, similar to the Bayesian or Akaike information criteria (2).
If the statistic d 2 evaluates to zero for all (or even multiple) models, then we are unable to make a selection between the models.This can be remedied by designing experiments that yield more informative data.For example, measuring more variables can reduce the dimension of (V M ) R ∩ (V D ) R ; the most informative situation is when (V Mi ) R ∩ (V D ) R = ∅ for all models.This indeterminacy can also be resolved by taking multiple measurements and minimizing the joint squared distance (see SI Appendix).

Problem 3: Parameter Estimation
Parameter estimation can be achieved by finding the point (a * , x * , y * ) ∈ (V M ) R that minimizes the value of (5).The point (a * , x * , y * ) is the maximum likelihood estimate under the given noise assumptions.The parameter estimate is then a * ; the estimate of the hidden variables, x * ; and the estimate of the "de-noised" outputs, y * .Of course, if the data and model varieties intersect, there will be one or more (possibly infinite) choices for (a * , x * , y * ).Otherwise, it is a matter of solving a polynomial system that yields the point(s) on (V M ) R nearest (V D ) R .This is described in more detail in the next section.

Geometry
In each of the problems stated above, we seek to minimize the distance between (V M ) R and (V D ) R or to find the intersection of these two sets.Standard methods for solving non-linear optimization problems are local in nature, i.e., only guaranteed to converge to a local minimum which may or may not be the global minimum.However, using NAG, we find all local extrema over C, necessarily including the global minimum.Due to this global benefit, NAG has been used before in statistical inference in the field of algebraic statistics (21,22).
Let V M ⊆ C k+n+m be the (complex) Zariski closure of (V M ) R and V D ⊆ C k+n+m be the (complex) Zariski closure of (V D ) R .We will refer to V M and V D as the model variety and data variety, respectively.
The problem of determining the intersection of V M ∩ V D is simply a matter of solving the polynomial system obtained by taking the union of the polynomials defining V M and the polynomials defining V D .This is handled directly by NAG.If the intersection is nonempty and positive-dimensional (complex curves, surfaces, etc), real points can be found using the polynomial homotopy method described in (23), a method based on symbolic algorithms in (24), and, more classically, on the decision method in (25).
The problem of finding the points on the two varieties nearest one another can also be stated in terms of a polynomial system, on which we then call NAG solvers to find solutions.A well-known necessary condition for local extrema is given by Lagrange multipliers.In the main text we assume that r + m = codim V M ; however, when this is not the case, the number of equations can be reduced (see SI Appendix).
then there exists λ := (λ 1 , . . ., λ c ), such that (a, x, y, λ) is a solution to the system Solving this system with NAG will provide us with all local extrema of (6) over V M , from which we may easily select the pair of nearest points.The geometry of zero sets of systems such as ( 7)-( 9) are described in (26).

Numerical Algebraic Geometry
Given a polynomial system F consisting of r polynomials and N variables, NAG packages, such as Bertini (17), PHCpack (27), HOM4PS-2.0(28), use polynomial homotopy continuation to provide probability-one numerical methods for finding approximations of all isolated complex solutions of F = 0 (points) as well as witness points on all positive-dimensional irreducible components of the solution set of F = 0.These methods are probability-one in that the computations include some randomization, and this randomization will yield theoretically correct results so long as the random numbers are not chosen from some measure zero set in the parameter spaces of potential choices (16,29).
If x ∈ R N is a real solution of F = 0, it is either isolated among the complex solutions or it lies on a positive-dimensional complex irreducible component.In the former case, the methods of NAG will find x and recognize it as real.In the latter case, x can be difficult to uncover.However, for our purposes, we usually only need to verify the existence of a real solution.In this case, we can find witness points on all positive dimensional components and then use the procedure in Section 2.1 of (23) to verify the existence of real points.

Algorithms
We present three related algorithms to address model validation, model selection, and parameter estimation (see Fig. 1).
The aim of the first algorithm, Algorithm 1, is to find the pair of points that minimize the distance between (V M ) R and ( then this is obtained by solving (7)-( 9); otherwise, additional techniques are required.
We note also simply that these techniques are indeed probability-one algorithms.The computations will necessarily be carried out in finite time and the steps proceed linearly (no loops), so the methods will necessarily finish.

Input
Model M, data D = {ŷ}, tolerance α Output yes or no The computation of the intersection V M ∩ V D in Steps 1 and 2 can be determined in several ways.The simplest way is by considering dimensions: if dim(V M ) + dim(V D ) exceeds the ambient dimension, they will almost surely intersect.If the ambient dimension is larger than the sum of the variety dimensions, they typically do not intersect.To compute the intersection (or check to see if it is nonempty), one could substitute ŷ − g(x) for (8) in the system (7)- (8).
In Step 2, if dim(V M ∩V D ) = 0, then the intersection of the two varieties consists of finitely many points; the condition ( one needs more sophisticated methods, such as those in (23).
To find the pair (z 1 , z 2 ) in Step 3, one may solve the polynomial system (7)- (9).If there is a positivedimensional set of (complex) extremal points, then the procedures in (23) could be used to determine if the set contains a real point.
If there are constraints on the variables or parameters, for example, if we seek to minimize (6) over the non-negative part of (V M ) R , then the algorithm is updated as follows: if the dim(V M ∩ V D ) = 0 or V M ∩ V D = ∅ then the Karush-Kuhn-Tucker equations as described in Section 3.1 of the SI Appendix should be used, while if dim(V M ∩ V D ) > 0, or if there is a positive-dimensional set of extremal points at Step 3, then the algorithm should return possibly.There are methods (30)(31)(32) for finding real points, curves, and surfaces within complex components of dimension 2 or less, but little is known about higher dimensions.

Algorithm 2: Model selection
In this case, there are several competing models, each with its own polynomial system.The algorithm proceeds much as in Algorithm 1, but iterated for the several models under consideration.If a threshold α is provided, one should first reject models that do not adequately support the data (d 2 ≥ p α ), then choose the model yielding the minimum value of d 2 (up to some complexity penalty).Various conclusions may be drawn, e.g., no model adequately fits the data or three models adequately fit the data and model M 1 provides the best fit.

Algorithm 3: Parameter estimation
Again, this algorithm is similar to the first.The input consists of only one model M and data D. It is assumed that there are unknown parameters and the goal is to find values of these parameters producing the best fit between M and D. The output of Steps 2 and 3 need to be adjusted appropriately.The output of Step 4 is simply z 2 since there is no α to be used for rejection.The method also simultaneously estimates hidden/unknown variables and "de-noised" outputs.

Simple Example
To illustrate Algorithm 1, consider a simple model with three variables x, y, z and three parameters a, b, c satisfying Let α = 0.1 and assume, for this example, that the variances on the errors are σ 2 i = 0.1 for i = 1, 2, 3.The model variety V M is simply an ellipsoid (Fig. 2A) where a, b, c describe the principal axes.Suppose we know that a, b, c = 1.For the case when the outputs are x and y and x = 0, ŷ = 0 and α = 0.1, Step 2 indicates that the data does fit the model, i.e. the model is compatible with the data.In this case, there are two real points in the 0-dimensional intersection V M ∩ V D (see Fig. 2B).For the same set-up, but with data x = 0, Step 3 indicates that this data possibly fits this model.Since there is a positive-dimensional intersection (Fig. 2C), it is possibly compatible (depending on constraints imposed by the user).For the same model and α, but different data x = 1.7, ŷ = 0, Step 4 yields points (1, 0, 0) and (1.7, 0, 0), so ||z 1 − z 2 || 2 > 0.4605 = p α , and the model is rejected (Fig. 2D).However when the data are x = 1.01, ŷ = 0, Algorithm 1 indicates model compatibility (Fig. 2E).Previous algebraic methods that required Gröbner basis calculations would result in a zero set and thus those approaches are not useful here.

Results
We demonstrate our methods on problems in biomedicine: cell death activation, synthetic biology, epidemiology and multisite phosphorylation.Each of the forthcoming applications can be written as a mass-action chemical reaction network, which has the form ẋ = f (a, x), studied at steady state: f (a, x) = 0 as in the problem statement.Throughout these real-world examples, we emphasize the pivotal computations in these methods, such as determining the dimension of the intersection V M ∩ V D and finding points in the intersection (see Fig. 1B).Our first two examples, cell death signaling and genetic toggle bio-circuits, demonstrate how to handle positive dimensional intersections using two different approaches, while the remaining two examples, HIV and MAP kinase signaling, highlight analysis of zero-dimensional and empty intersections.In the following examples, we are interested in results that can be interpreted biologically, therefore we restrict our analysis to non-negative real solutions.

Cell death activation
We demonstrate model compatibility (Algorithm 1) on an example from receptor-mediated programmed cell death, which is initiated by the activation of death receptors upon the detection of extracellular death ligands (33)(34)(35)(36).We consider, in particular, the "cluster" model of (12), which was inspired by crystallographic data (37) and describes the recruitment of receptors by ligands into local self-activating clusters capable of bistability.
The cluster model is a system of 3 degree-4 polynomials in the form of (2) in 3 variables (representing various receptor states) and 6 rate parameters, supplemented by ligand and receptor conservations (see SI Appendix).We assume that we can measure the total ligand and receptor concentrations, which may be considered experimental inputs, as well as the concentration of active receptors.We do not assume access to the concentrations of other individual receptor states nor to any of the rate parameters.
A steady-state data point was simulated from the model with all parameters and initial concentrations drawn independently and identically from the log-normal distribution ln N (0, 4), then combined and corrupted with i.i.d.noise from N (0, 0.1) to obtain ŷ.The real model and real data varieties intersect in the positive orthant with a distance zero and hence the model is indeed compatible with the data.

Synthetic biology and experimental design
We demonstrate an example from synthetic biology with excess intersection (dim(V M ∩ V D ) > 0).A goal in synthetic biology is to design or modify existing biological systems with new features according to specific design criteria.Reverse engineering of biological systems often includes modules (such as feedback loops) and how these are interconnected are described by different circuits (models).Understanding differences between bio-circuit implementations is crucial; therefore we compare three bistable bio-circuits models analyzed in (38): monomer-dimer toggle circuit (M 1 ), dimer-dimer toggle circuit (M 2 ), and single operator gene circuit (M 3 ), which were initially presented in (39)(40)(41).The model variables include genes (X i ) and proteins (P i ) where i = 1 for M 3 or i = 1, 2 for M 1 , M 2 , as well as these species complexes (e.g., P i P i , X j P i P i ).These variables interact following mass-action kinetics and form systems of polynomial differential equations where M 1 , M 2 , and M 3 have 7, 8 and 6 model variables, respectively, and 10, 12, and 9 kinetic parameters, respectively.The models can be reduced (given in SI Appendix) by assuming that the total amount of gene 1 (X 1tot ) and gene 2 (X 2tot ) is conserved.
Suppose that the total amounts X 1tot and X 2tot and specific protein synthesis and degradation parameters k bas1 , k bas2 , k deg1 , and k deg2 are known.Since protein concentrations are often measurable, we assume that our data are P 1 , P 2 , and their complexes P 1 P 1 , and P 2 P 2 .The aim is to select the best model M 1 , M 2 , and M 3 given the data.We simulate steady-state data from the dimer-dimer toggle model (M 2 ) and add Gaussian noise from N (0, 0.1).We find that all three have positive dimensional intersections, where the dimension of the intersections are: Clearly, all three models are compatible with the data, thus one can only select a "best fit" model using data-independent measures, e.g.number of parameters, dimension, etc.
In fact, the dimensions of the model and data varieties can help us design more informative experiments for model selection.For example, the variety associated to the monomer-dimer toggle model lives in a 15 dimensional ambient space, V M1 ⊆ C 15 , and has dimension 6, while the data variety V D ⊆ C 15 has dimension 12; thus by the Dimension Theorem, we can determine that dim(V M1 ∩V D ) ≥ 6+12−15 = 3.This dimension calculation provides guidance towards the number of minimal additional variable and parameter values that must be measured to ensure V M1 ∩ V D = ∅, i.e., at least four more variables and parameter values must be known.
Suppose we can experimentally measure four forward biochemical reaction rate constants (e.g., k cF , k kF , k nF , and k kR ), then V M1 is cut down by four dimensions and does not intersect V D .We get similar results for the model varieties associated to M 2 and M 3 provided we measure rate constants specific to these models (see SI Appendix).Now that all the intersections are empty, we run Algorithm 2 and find that the sum of squares (Eq.( 5)) for each model are as follows: d 2 1 = 0.2262, d 2 2 = 7.34 × 10 −7 , and d 2 3 = 0.3040.Therefore, we select the M 2 model, which is indeed the true model.

HIV Progression
We demonstrate parameter estimation (Algorithm 3) on an example coming from epidemiology.We use a model that includes long-term HIV dynamics from initial viremia, latency, and virus increase (42), based on (43).In the model (see SI Appendix), the HIV virus inhibits the CD4+T cell population while promoting macrophage proliferation, which in turn houses the replicating virus.As macrophages proliferate, the virus reservoir increases, so the model offers a description of HIV patient progression to AIDS.Model variables x are uninfected CD4+T cells (T ), infected CD4 + T cells (T i ), uninfected macrophages (M ), infected macrophages (M i ), and HIV virus population (V ).
Hernandez-Vargas et al. show that the model can have two real equilibria, one of which is stable, representing patients that are "long-term non-progressors" (42).The parameters a are (s 1 , s 2 , k 1 , . . ., k 6 , δ 1 , . . ., δ 5 ), where s i represents synthesis of T cells and macrophages, k are rate constants describing interactions between variables x, and δ i represents natural death.For this example, y = x.We estimated the natural death of the virus, parameter δ 5 , using the long-term non-progressors steady-state value (Table 3 of ( 42)) and adding noise to each variable ∼ N (0, 1).By Algorithm 3, the data variety and model variety do not intersect.We find the closest point and estimate δ * 5 = 2.99876 (true value of δ 5 = 3).

Multisite phosphorylation with experimental data
We examine phosphorylation mechanisms of cellular signaling with experimental data, and demonstrate model selection (Algorithm 2) and parameter estimation (Algorithm 3).We focus on phosphorylation, a key cellular regulatory mechanism that has been the subject of extensive study, both experimentally and theoretically ((44) and references therein).An area of interest is the mechanism by which a kinase phosphorylates a two-site substrate, either distributively, where the kinase can add at most one phosphate before dissociating, or processively, where it can add both phosphates in sequence.The MAPK/ERK pathway is a well-known system for studying phosphorylation, whereby MEK (kinase) phosphorylates ERK (its substrate).Aoki et al. (45) showed experimental evidence, while working with polynomial models, that the mammalian MAPK/ERK pathway acts distributively in vitro but processively in vivo.
We compare these distributive and processive models against the in vivo data reported in the same study.The distributive model consists of 12 molecular species and 17 mass-action reactions, while the processive model has 14 species and 18 reactions (species correspond to variables, each reaction corresponds to a parameter).The data take the form of 36 concentration measurements of three aggregate phosphoforms over a range of 12 EGF stimulation levels.All model parameters are calibrated using in vitro estimates by (45), except the parameter k 1 representing EGF loading, which we estimate for both models (see SI Appendix).
Next we perform model selection by running Algorithm 2 on each data point individually and select independently for each run the preferred model (more details in SI Appendix).Under low EGF stimulation, the best model estimates are nearly identical with a slight preference for distributive.At high EGF stimulation, the models are identical with no preference for one model over the other.These results can be justified by noting that the main distinction between distributivity and processivity is nonlinear switching behavior (i.e., a sigmoidal response curve), and this only occurs at intermediate stimulations.However, at medium EGF stimulation (see Fig 3), there is a preference for the processive model, which supports the findings in (45).

Conclusion
The problem of determining whether given real-world data fits one or more given mathematical models is challenging.When a model is defined by algebraic (polynomial) functions, the methods of NAG may be employed to study the geometry underlying the model and data.In particular, these methods are useful for model variables observed at steady-state.We demonstrated this numerical and geometric framework for comparing models with experimental dose-response data in MAPK/ERK pathway and highlighted that the intermediate EGF doses were the most informative for model selection, complementing another finding that model selection results can be very sensitive to experimental parameters (46).
Despite the difficulties associated with positive dimensional components, and limitations in analysis, we can reproduce compatible models, and furthermore can predict additional information, such as measurements of parameters, that are necessary for selecting models.Our geometric investigations of positive dimensional components may perhaps relate to algebraic analyses for biochemical models, such as model identifiability or matroids for experimental design (14,47,48).
There are further directions to be considered in this vein, aside from making the existing computational methods more efficient.First, there would be great value in developing strictly real geometric methods for solving polynomial systems such as those that appear in this article.Some such techniques exist, but only in very special circumstances.Second, there would be much value in developing effective numerical methods for treating inequalities.It should be noted that the methods described in (49) and the references therein will incorporate such constraints, though the cost of such computations restricts their use to relatively low dimensions.Finally, there is likely much to be gained from considering the geometry underlying models not defined by algebraic functions.Algebraic geometry provides very clean, well-understood structures, paving the way for numerical methods.Differential geometry or topology could lead to similarly useful techniques for model selection and parameter estimation.
6 Materials and Methods

Numerical algebraic geometry
General references for NAG include (16,29), with the latter doubling as a user manual for the software package Bertini.For computations, we used Bertini 1.4 and Macaulay2 version 1.6.

Data generation
Data simulated from cell death activation, synthetic biology and HIV models were performed in MATLAB R2014b using ode15s.
Find point in ?
Supporting Information for: Numerical algebraic geometry for model selection and its application to the life sciences Elizabeth Gross

Geometry
In this section we briefly review the some of the fundamental geometric concepts needed for the methods in the main text.

Numerical algebraic geometry for isolated solutions
Numerical algebraic geometry refers to the use of numerical methods, particularly homotopy continuationbased methods, to compute approximations to solutions of polynomial systems.In other words, given a polynomial system F : C N → C n with n equations in N variables, numerical algebraic geometry seeks to 1 arXiv:1507.04331v3[q-bio.QM] 1 Apr 2016 find numerical approximations of all z ∈ C N such that F (z) = 0.It may be the case that the solution set of F has infinitely many such points (curves, surfaces, etc.), in which case the data structure that encodes the solutions is called a witness set.See the description of the numerical irreducible decomposition below for more on this.For simplicity, we assume in this section that N = n.The books [1,2] and the references therein provide much more detailed explanations than those included in this section, from more mathematical and computational perspectives.
The core technique for most numerical algebraic geometry algorithms is homotopy continuation.The idea of homotopy continuation is to cast the polynomial system to be solved, say F : C N → C N , as a member of a parameterized family of polynomial systems, H : C N × C → C N , called a homotopy H(z, t) with parameter t ∈ C, that includes one polynomial system G : C N → C N that is easily solved and has the special property of having "enough" isolated solutions.In this document, we use the Bertini standard assumption that H(z, 1) = G(z) and H(z, 0) = F (z), i.e., t marches from 1 to 0. There are several canonical options for the construction of such a homotopy, and the reader is encouraged to consult the general references above and the references therein for further details.
As t varies from 1 to 0, some results from algebraic geometry tell us that the solutions of the polynomial system H(z, t) = 0 vary continuously and generally stay distinct until t = 0, where they may converge to solutions of F (z) = 0 or diverge.More specifically, there is a measure zero subset of t ∈ C, meaning a finite set of points in this particular parameter space, over which two or more solutions coalesce.Such occurrences are thus probability zero events and, furthermore, can be detected on the fly and avoided.Said more technically, there is a Zariski open, dense set of the parameter space above which the solution set is finite and consists of a fixed number of solutions.Here, "Zariski" refers to the Zariski topology, for which basic open sets are the complements of solution sets of polynomial systems.
In practice, t is moved in discrete increments, not continuously.For each solution at t = 1, a path of solutions is tracked using numerical predictor/corrector methods as t advances to 0. Implementations typically utilize adaptive step lengths and adapative precision.There are far too many details about this procedure to give a thorough explanation here.Instead, refer to the references above (and those therein) for further details.
Ultimately, the output of this procedure is a superset of numerical approximations of the isolated solutions of F (z) = 0, possibly including approximations to points lying on positive-dimensional components, if any.It is important to note that this procedure necessarily works over C and finds all complex solutions.Real solutions could be buried somewhere within the complex solutions, and it is particularly difficult to extract these outside of the zero-dimensional case.However, methods do exist to extract such a real point [3][4][5]; here we use the method of [5], which is guaranteed (with probability one) to minimize the distance of a prescribed real point to each real connected component of the variety defined by F (z) = 0.

Numerical algebraic geometry for positive-dimensional solution sets
For solution sets of positive dimension (curves, surfaces, etc.), there is an extension of homotopy continuation referred to as the numerical irreducible decomposition (NID).As opposed to the case of systems of linear equations (at most one solution component of one dimension), there may be many components of many different dimensions.For example, one solution set might consist of seven components of dimension four, five surfaces, three curves and 15 isolated points.Furthermore, components may be singular, meaning that the Jacobian matrix is rank-deficient throughout the component.
Technical definitions of dimension, irreducible component, and the like go a bit beyond the scope of this paper.It is enough to know that each "piece" of a solution set has a fixed dimension (e.g., a curve has dimension one, a surface two, etc.) and by the dimension of a solution set of a polynomial system of equations, we mean the maximum of the dimensions of the irreducible components.
The numerical irreducible decomposition of a solution set of a polynomial system consists of a catalog of the dimensions and degrees of each irreducible component, along with a set of witness points on each component.By degree, we mean the number of points in the intersection of a component with a randomlychosen affine linear space of complementary dimension.The witness points on a component are then exactly these points (and thus depend on the choice of linear space).One fundamental result from algebraic geometry is that an irreducible component will almost always intersect a complementary-dimensional linear space exactly in a set of points and that the number of points is the same for almost any choice of linear space.Again, this can be stated as a probability one guarantee or with Zariski open, dense sets, but we choose not to be that technical here.

Software for numerical algebraic geometry
Various numerical algebraic geometry software packages have been produced over the years.Currently, there are three main options: PHCpack [6], HOM4PS-3 [7], and Bertini [8], each with their own benefits and drawbacks.In this article, we used Bertini exclusively.

Geometry specific to presented algorithms
Theorem 1 of the main text should be familiar to those trained in multivariate calculus; this is essentially the method of Lagrange multipliers.Geometrically, this system of equations forces the normal directions of the objective function and the constraints to line up in the same direction (up to some scalar(s), the Lagrange multiplier(s)).
When working with an irreducible component of the solution set of a system of polynomial equations, it is often useful to deal with a complete intersection.Said simply, the idea is that computations can be more difficult if there are more equations than necessary.To see why this might be true, let us consider an example from linear algebra.Suppose we have a single linear equation in three variables, defining a plane.Now suppose we consider a system of two equations consisting of that equation and twice that equation (having the same solution set).Then the matrix of coefficients of this linear system is not full rank (not a desirable situation) and we have two equations defining a geometric object that could be described by a single equation.
In the nonlinear setting, the situation is quite similar.Having "too many" equations leads to an undesirable rank-deficient Jacobian matrix.Suppose polynomial system F : C N → C n has an irreducible component X of dimension N − m (codimension m).Then, again with probability one, the polynomial system A • F has X as an irreducible component but has the "correct" number of equations, where A is a random constant matrix with n columns and m rows.Here, "correct" means the number of equations matches to codimension m.We will refer to this method as squaring up.
Finally for this section specific to the algorithms developed in the main text, we require the user to check two geometric facts, the meaning of which may not be entirely clear.
1. V M ∩ V D refers to the intersection of the model and data varieties, as defined in the main text.To find the intersection of two solution sets, it is sufficient to simply solve the system consisting of all equations appearing in the systems for V M and V D , i.e., the union of those two polynomial systems.There are more sophisticated methods, but this is sufficient.
2. The user must determine whether the intersection of (V M ) R and (V M ) R is empty.By this, we simply mean that one should search for real points in the intersection just described, e.g., using the method of [5].

Choice of test statistics and parameter estimates
In this section, we justify our procedures for model validation and parameter estimation.

Maximum likelihood
Here we justify the assertion that the test statistic given in the main text is related to likelihood maximization.
Consider first, for simplicity, the case of a single data point ŷ = (ŷ 1 , . . ., ŷm ), which we assume is a perturbation ŷ = ξ + of some unknown true value ξ = (ξ 1 , . . ., ξ m ), where each component i of the error = ( 1 , . . ., m ) is an independent zero-mean Gaussian random variable with variance σ 2 i .We are interested in computing the probability that ŷ comes from a given model as defined by a model variety V M .A point on V M has the form (a, x, y), where a = (a 1 , . . ., a k ) are the model parameter values, x = (x 1 , . . ., x n ) are the variable values, and y = (y 1 , . . ., y n ) are the outputs.The probability that ŷ comes from a given point (a, x, y) ∈ V M , i.e., that ŷ is a perturbation of y where (a, x, y) ∈ V M for some a and x, is then This is also called the likelihood L(a, x, y) of (a, x, y) and we wish to find its maximizer over all (a, x, y) ∈ V M .This can equivalently be done by considering the log-likelihood, which gives log L(a, x, y) i by normality.The maximizer (a * , x * , y * ) can therefore be found by solving the optimization problem where the optimum is precisely the test statistic.The values a * , x * , and y * are the maximum likelihood estimates for, respectively, the parameters (estimation), the unobservable variable values (inference/recovery), and the true output values (filtering/denoising), The test statistic d 2 itself also has a useful interpretation as follows.Suppose that ŷ comes from a point (a, x, y) ∈ V M .Then σ 2 i by definition.But regarding each ŷi as a random variable, each term (ŷ i − y i )/σ i in the summation above is standard normal.Hence the right-hand side has a chi-squared distribution with m degrees of freedom (χ 2 m ).The inequality should be interpreted by regarding d 2 as a random variable subject to the same source of randomness.This can be written somewhat clearer as where we have made explicit the underlying dependence of both sides on the same random realization ω.
The inequality then holds for each value of ω.Consequently, we conclude that where p α is the upper α-percentile for χ 2 m .This can be used to test the hypothesis that ŷ comes from V M .The test statistic is also related to the log-maximum-likelihood as log L(a * , x * , y which is a useful quantity for model selection via, e.g., the Akaike or Bayesian information criteria.Now consider the case of multiple data points {ŷ (j) } p j=1 .As before, we assume that each ŷ(j) = (ŷ j) , where each (j) i is an independent zero-mean Gaussian random variable with variance σ 2 j,i .Instead of searching for one point on V M , we now have to search for p points (a, x (j) , y (j) ) for j = 1, . . ., p all with the same parameter values (since they come from the same fixed model realization).The probability that ŷ(j) comes from (a, x (j) , y (j) ) for j = 1, . . ., p is then Pr(ŷ (1) , . . ., ŷ(p) | a, x (1) , y (1) , . . ., x (p) , y ≡ L(a, x (1) , y (1) , . . ., x (p) , y (p) ) by independence.This is essentially the same as before except that we now loop over each coordinate of each data point.Therefore, the maximum likelihood estimates (a * , x (1) * , y (1) * , . . ., x (p) * , y (p) * ) can be obtained by solving 1) ,y (1) ,...,x (p) ,y (p) )∈VM The same arguments go through and we find that Pr(d 2 ≥ p α ) ≤ α for p α the upper α-percentile for χ 2 mp .The log-maximum-likelihood is related to d 2 as log L(a * , x (1) * , y (1) * , . . ., x (p) * , y (p) * ) = 1 2

Algorithm modifications
The algorithms presented in the main text are in their simplest form.Some applications require modifications, particularly if there are constraints on the variables or parameters.

Solving the constrained optimization problem
In many common settings, there exist constraints on the variable and parameter spaces.For example, in chemical reaction networks, the rate parameters are all assumed to be positive.Thus, when positivity or other constraints are present, instead of finding the weighted squared distance between two varieties, we are finding the weighted squared distance between two semi-algebraic sets, i.e. sets defined by polynomial equalities and inequalities as opposed to just polynomial equalities.Indeed, if we let S M ⊂ (V M ) R denote the semi-algebraic set associated to the model, e.g.
, then the appropriate statistic is: In the case when a bound on the statistic d 2 is sufficient, then no additional work is needed.One can solve the system from Proposition 1, keeping in mind that the weighted squared distance between the closest pairs of points returned would be a lower bound on d 2 .If the closest point to (V D ) R in (V M ) R is also an element of S M , then the squared distance would be exactly the statistic d 2 .
When the exact value of d 2 is needed, then one should solve the Karush-Kuhn-Tucker (KKT) system of equations.Let F 1 , . . ., F r , h 1 , . . ., h s be polynomials in the ring R[a 1 , . . ., a k , x 1 , . . ., x n , y 1 , . . ., y m ].Let S M be the semi-algebraic set of all (a, x, y) ∈ R k+n+m that satisfies F i (a, x, y) = 0 for i = 1, . . ., r h i (a, x, y) ≤ 0 for i = 1, . . ., s Let λ 1 , . . ., λ r , µ 1 , . . ., µ s be indeterminates (these are the KKT multipliers).The KKT system is f = 0 (1) In order for (a, x, y) to be a critical point, the point (a, x, y) must be a solution to this system and satisfy the inequalities h i ≤ 0 for i = 1, . . ., s and µ i ≥ 0 for i = 1, . . ., s.Thus, we can find the global minimum by using numerical algebraic geometry to solve the system defined by equations ( 1) -( 5), then filtering the solutions appropriately.This combination of numerical algebraic geometry and the KKT equations was first employed in [9].
Alternatively, in some situations, it may be more efficient to minimize the objective function over (V M ) R and then check the boundary of S M .We describe this method for when there are non-negative constraints on all indeterminates.After solving the constrained optimization problem on (V M ) R , we assume, in an orderly fashion, that one of the indeterminates, say x i , is zero.In this case, we are required to minimize the distance of V M ∩ {x i = 0} to V D .This involves solving a smaller and related constrained optimization problem.In total, if there are N indeterminates, there are 2 N − 1 combinations to set to zero.This amount to solving 2 N − 1 Lagrange multiplier systems.
One observation is that the number of systems that need to be solved explodes when N is large.Note however that the dimension of V M ∩ {x i = 0} is less than or equal to the dimension of V M , with the inequality being strict when V M {x i = 0}.Geometrically, there are no longer degrees of freedom in the variable x i so the dimension is reduced.Furthermore, we can expect the dimension to be reduced by one.As we impose additional constrains, x j = 0 for i = j, the dimension may drop further.
In practice, there are diminishing returns as you begin setting x i to zero.That is, there exists some positive integer k such that V M ∩ {x i1 = . . .= x ik = 0} is zero dimensional for some indexing set i 1 < • • • < i k .In this case, measuring the distance of V M ∩ {x i1 = . . .= x ik = 0} to V D using a Lagrange multiplier method is unnecessary, and, as we continue imposing additional constraints on V M , either the intersection is non-empty and every observable variable is set to zero, or the model variety becomes empty.Thus, it becomes redundant or unnecessary to check additional cases.
For example, in Section 4.4.3, the MAPK model variety V M is one-dimensional and V M ∩ {x i = 0} becomes zero-dimensional for each i.For the cases where x i = x j = 0 for i = j, the intersection of the corresponding linear spaces with V M becomes empty.Even though there are 2 16 − 1 = 65, 535 boundary cases to check, there are really only 16 relevant cases.See Section 4.4.3 for more explicit details on how this calculation carried out.

Removing extinction components
Given a model, it is quite common that the model variety is not irreducible but instead is the union of several irreducible components.In applications, it may be preferred to remove from consideration components that lie entirely in a coordinate hyperplane, since, in such components, one or several of the parameters and/or variables are equal to zero.For example, in a chemical reaction network, such a component is called an extinction component [10] since it captures the situation where one or more of the reactants have "run out."It is common to want to avoid extinction components when estimating parameters.
Removing components where a parameter or variable is equal to zero throughout the set from consideration can be done algebraically with saturation.In particular, if I M is the defining ideal of the model V M one should compute This procedure can be performed using the saturate command in Macaulay2.
To estimate parameters such that the best estimate corresponds to a point not on an extinction component of V M , one should modify Algorithm 3, replacing V M with V(I M ain ).

Results
We provide details for the calculations of examples in the main text.All code is available at: http://www.math.sjsu.edu/∼egross/NAGModelSelection/AuxillaryFiles.

Cell death activation
We provide the details of the calculations for model compatibility of the cell death cluster model.This subsection includes detailed information regarding the solving schemes available in NAG software that were utilized.A summary for the practioner can be found at the end of the subsection.
The model describes activation of apoptosis by death receptor Fas mechanisms [11].The model includes constitutive receptor opening and closing, pairwise open Fas stabilization, higher-order open Fas stabilization enabled by FasL, and ligand-induced receptor opening.According to its conformational states, Fas is assumed to be one of three species: closed (X 1 ); open, unstable (X 2 ); and open, stable (X 3 ), i.e., active and signaling.Furthermore, let the ligand FasL be denoted by L. Then the model has the reactions for i ∈ {2, 3}, j = 1, . . ., i, and k = 1, . . ., j.The first reaction describes spontaneous receptor opening and closing.The second reaction describes constitutive destabilization of open Fas.The third reaction describes cluster-stabilization by open Fas, independent of the presence of FasL.The fourth reaction describes clusterstabilization events enabled by FasL.
Assuming mass-action kinetics, the reactions can be translated as follows: where l x 2 x 3 l, where v i are the reaction velocities for the variables x i , and lowercase letters denote the concentrations of their uppercase counterparts.
The model parameters for the cell death cluster model are l , k l ) and the variables are The outputs are y = (λ, ρ, ζ) representing, respectively, the total ligand concentration, the total receptor concentration, and the total downstream "death signal", as given by the equations We set the model equations ( ẋ1 , ẋ2 , ẋ3 ) to zero, and, together with equations ( 6)-( 8), we obtain the defining equations for the model variety V M .The ambient space that V M is contained in has dimension 14.This space has coordinates defined by both the model parameters a, the variables x and the outputs y.
Given an observable data point ŷ = ( λ, ρ, ζ), we define the data variety as: for the clustering model.Note that V D has dimension 11 since there are no degrees of freedom in the variables λ, ρ, or ζ.We first compute a numerical irreducible decomposition (NID) of V M using Bertini; this will aid in understanding V M ∩ V D .One can verify, after computing the NID for V M , that V M is a 9-dimensional complex algebraic set of degree 10 (file name: Cluster Model NID ).Now suppose we are given the following data point (taken from the model without noise): One can then verify using the NID that V M ∩ V D = ∅ (file name: Cluster Model Data NID).That is, the intersection of the model variety and data variety is nonempty.Specifically, V M ∩ V D is a 6-dimensional complex algebraic set of degree 5. Adding noise to the coordinates of ŷ taken from N (0, 0.1) did not affect the dimension or degree.Again, we got that V M ∩ V D has dimension 6 and degree 5.
Since we are interested in model compatibility, our goal is to find at least one nonnegative point in The above computation at least provides evidence that this is the case, but it may be possible that V M ∩ V D does not contain any nonnegative real points or any real points at all for that matter.We will approach this problem using the methods described in [5].If V M ∩ V D contains a real nonnegative point then we cannot reject the model and may conclude that the model is compatible with the data.If (V M ) R ∩ (V D ) R does not contain a real point then we can try and use a more general Lagrange multiplier method similar to the one employed in Section 4.4.3 in dealing with model selection.
We first randomly select a real, positive point: = 6.491154749564521 x 1 = 7.317223856586703 x 2 = 6.477459631363067 x 3 = 4.509237064309449 where each coordinate is chosen uniformly on the interval [0, 10].This point will determine the observable, i.e. output, variables λ, ρ, and ζ using equations ( 6)- (8).Call this point (a , x , y ) ∈ R 14 .For the time being, the left and right endpoints of each subinterval [0, 10] have been chosen arbitrarily.It is unclear, at this time, how the endpoints or length of the interval affects the performance in finding nonnegative real points contained in (V M ) R ∩ (V D ) R using the methods described below.Our aim then is to solve the constrained optimization problem: Geometrically, we are minimizing the distance between the chosen point (a , x , y ) and We will refer to the system defining (V M ) ∩ (V D ) as f * (a, x, y).Squaring up the polynomial system will be a necessary step when utilizing the perturbed regenerative solving scheme.Squaring up was briefly described in Section 1.4 but we will give additional details here.First notice from previous computations that V M ∩ V D has codimension 14 − 6 = 8.Thus, there exists a nonempty Zariski open set A ⊆ C 8×9 such that for every matrix A ∈ A, we have V M ∩ V D ⊆ V(Af * (a, x, y)).This means we may take 8 random C-linear combinations of the equations defining V M ∩ V D and, with probability one, still cut out at least V M ∩ V D .It is sufficient to sample the elements of A uniformly along the complex unit circle.
As a side note, we may take additional steps to reduce complexity by replacing the matrices A with matrices of the form [I 8 | b] ⊆ A where I 8 denotes the 8 × 8 identity matrix and b denotes a 8 × 1 column vector who elements are sampled uniformly along the complex unit circle.This has the effect of adding a random multiple of the last function to each of the other functions.One may verify a posteriori if a point x ∈ V(Af * (x)) is also in V M ∩ V D by function evaluation of f * (a, x, y).
The polynomial system for the optimization problem ( 10) is: We call λ = (λ 1 , . . ., λ 8 ) the Lagrange multipliers for V M ∩ V D .This is a system of 22 variables and 22 equations.Written more compactly we may write equations ( 11)-( where J h (a, x, y) is the Jacobian matrix of h(a, x, y).The Lagrange multipliers in equations (16) need not be real since we are taking C-linear combinations of f * (a, x, y).
Regeneration is a numerical algebraic geometry method we found to be most relevant to solve equations ( 15)-( 16) and is implemented in Bertini.Regeneration uses a linear product homotopy scheme in which each equation is built up one at a time.Depending on the ordering and structure of each equation, that can lead to huge computational savings.However, regeneration is restrictive in that it may not find all singular isolated solutions to equations ( 15)- (16).For the finer details of regeneration see reference [2].
One small adjustment we can do to solve this problem is to first solve a slightly perturbed problem (this is the strategy employed in [5]).Indeed, there is a nonempty Zariski open set Γ such that for every γ ∈ Γ the solutions to: will be nonsingular.The benefit to first solving this system is that regeneration will now find all solutions.
After the solutions to equations ( 17)- (18) have been computed, one can use a parameter homotopy to compute all isolated solutions (15)-( 16), which may contain singular solutions.We briefly describe the parameter homotopy employed following regeneration.The solutions of equations ( 17)-( 18) lead to the solutions of equations ( 15)-( 16) through a collection of homotopy paths where each homotopy path as functions of t are solutions to the straight-line homotopy function: for t ∈ (0, 1] ⊂ R. As t → 0, we obtain numerical approximations to the solutions of equations ( 15)-( 16).Additional details on parameter homotopies can be found in [1,2].A basic implementation of parameter  Timing summaries for the clustering model can be found in Table 1.In all cases, we have employed the use of intrinsically defined variables (see Appendix F.1.2 of [2]).These timings include computing the numerical irreducible decomposition of V M , the numerical irreducible decompsition of V M ∩ V D , computing the solutions to equations ( 17)-( 18) using regeneration, and computing the solutions to ( 15)-( 16) using the parameter homotopy.We found it most appropriate to utilize Bertini in serial for each computation except for regeneration which was done in parallel.Serial runs were done using a Apple MacBook Pro with 2.4 GHz Intel "Core i5" processor.Parallel runs were done using 24 (2.67 GHz Xeon-5650) compute nodes on the CentOS 5.11 operating system.In total, after reviewing the solutions, there are three solutions that correspond to real points contained in (V M ) R ∩ (V D ) R .Among the three real solutions, two solutions are nonnegative.Solutions are listed in Table 2.One can verify that these are indeed solutions to (V M ) R ∩(V D ) R by function evaluation of f * (a, x, y).
We may conclude from these computations that the clustering model V M is compatible with the observable data ŷ.It may be the case that there is a very large number of data points, ŷ, for which we would like to determine model compatibility.Fortunately, by employing a parameter homotopy scheme we may solve this problem rapidly where given each data point, the computations will be on the same order as the parameter homotopy solve in Table 1.
In summary the steps for model compatibility for the clustering model are as follows: 1. Determine the dimension of V M ∩ V D using the NID.

Using the information gathered in
Step 1, select a random point whose coordinates are sampled The model dimer-dimer toggle circuit (M 2 ) system is: The model single-operator positive feedback circuit (M 3 ) system is: In this example, we suppose that the total amounts X 1tot and X 2tot and specific protein synthesis and degradation parameters k bas1 , k bas2 , k deg1 , and k deg2 are known and we assume that our data are measurements of P 1 , P 2 , and their complexes P 11 , and P 22 , i.e. y = (P 1 , P We add Gaussian noise from N (0, 0.1) and then find dim(V Mi ∩ V D ) for i = 1, 2, 3. We can compute the dimension of each intersection using the dim command in Macaulay2 or by computing a numerical irreducible decomposition in Bertini; we find: Some further computations are required to find dim((V Mi ) R ∩(V D ) R .Specifically, we need to find real points in each intersection and determine whether or not those points are smooth.Computing the dimension of the real part of the intersections is more work than necessary for Algorithm 2, however, it provides an illustrative example on how to work with real varieties and the algorithm in [5].
Let f (i) = 0 be the polynomial system defining V Mi ∩ V D for i = 1, 2, 3 and let w (1) ∈ R 17 , w (2) ∈ R 20 , w (3) ∈ R 15 be random points.Let x (i) be the vector of indeterminates (unknown parameters and variables) for ith model, and let c i be the codimension of V Mi ∩ V D .We can find a real point on every component of each interesection, by solving the system: Table 4: Reactions for HIV model.The published parameter value is used from [17], see references therein.

Description
We estimated the natural death of the virus, parameter δ 5 , using Algorithm 3 with main component V 1 in place of V M (see Section 3.2).For V D , we used the long-term non-progressors steady-state value (Table 3 of [17]) and added noise to each variable ∼ N (0, 1).In particular, the data variety V D is defined by the equations: For s 1 , s 2 , k 1 , . . ., k 6 , δ 1 , . . ., δ 4 , we treated these parameters as known using the values from Table 1 of [17].
The varieties V 1 and V D do not intersect, which we can confirm with Macaulay2.Using Bertini, we solve the following system from Theorem 1:

Multisite phosphorylation
Here we describe the details for the multisite phosphorylation model selection and parameter estimation computations.First we describe the relevant biology, next we present the mathematical models of the distributive and processive mechanisms, then we apply our model selection method using data from [19].
We also estimate the relationship between the EGF concentration and activation of the pathway described by the parameter k 1 (see Table 11).

Biology of MAP Kinase system
Many cellular decisions are governed by molecular post-translational modifications.One type of modification, phosphorylation, is the addition of a phosphate group to a site of a substrate by an enzyme called a kinase.Some proteins (substrates) require multiple phosphate groups to be added by the kinase before the protein in activated/de-activated by these modifications.One well-studied signaling system is the MAP Kinase pathway, with kinase MEK and its substrate ERK; however the mechanism by which the phosphate group is added has been debated.Either MEK could phosphorylate ERK, disassociate and then phosphorylate again, called distributive; or MEK could bind and phosphorylate in sequence, called processive.Aoki et al [19] showed experimentally (with mathematical models) that this mechanism is different in vitro than in vivo.This experiment included 12 different levels of EGF stimulus ranging from 0.0244140625 ng/mL to 50 ng/ML.EGF actives cRAF which then phosphorylates MEK and finally doubly phosphorylates ERK.The data are measurements of three replicates of nonphosphorylated ERK (np-ERK), tyrosine monophosphorylated ERK (pY-ERK), and doubly phosphorylated ERK (pTpY-ERK), at each stimulus level.These data are given as percentage of total ERK (ERK), so we use the concentration measurement for each of these ERK states.

Mathematical models
The model variables and parameters are given in  The variables for the processive model are the same as for the distributive model except for two additional variables x 13 , x 14 .The reaction velocities are given in Table 6 and the corresponding equations are given in Table 7.Note in Table 7 that there are various conserved species concentrations in addition to the ODEs.We use the in vitro parameters estimates from Table S2 in reference [19] for k 2 , . . ., k 27 , c 1 , c 2 and the conserved quantities MEK tot , cRAF tot , ERK tot , as given in Table 8.The remaining parameter, k 1 , describes the rate of MEK phosphorylation and depends on the level of EGF stimulation, which varies throughout the data.Thus, we left it as a free variable and estimated it as a byproduct of distance minimization (11).The output variables are np-ERK, pY-ERK, and pYpT-ERK, which are sums of species concentrations.For the distributive model, the output equations are:   algebraic geometry approach, we should solve the system ( 27)- (30), and then solve the system again 16 more times, setting one of x 1 , . . ., x 12 , a 1 , y 1 , y 2 , y 3 to zero each time.We know that this will be sufficient to check the boundary, since the intersection of V Md with any coordinate hyperplane is zero dimensional.
The data set we used can be found in the file (aokidata.txt).There are 36 sets of data points where each set consists of a triple ŷ = ( np-ERK, pY-ERK, pYpT-ERK).Each data point defines a data variety, which we will denote (V Di ) R for the ith data point.
The polynomial system ( 27)-(30) is a classic example of a parameterized system of polynomial equations with parameter ŷ (notice that this is different than the rate parameters of the ODE).Let J T f (a, x, y) be the Jacobian matrix of f (a, x, y).The theory for parameter homotopies states that there is a nonempty Zariski open set U ⊂ C 3 such that for every y ∈ U the nonsingular solutions of: f (a, x, y) = 0 (31) 0 y T − ( y ) T = J T f (a, x, y)λ T lead to solutions of equations ( 27)-( 30) by homotopy paths.Each homotopy path is the set of zeros of the straight-line homotopy function of t ∈ (0, 1] as t varies from 1 to 0: f (a, x, y) 0 y T − t( y ) T + (1 − t) y T − J T f (a, x, y)λ T .
The benefit of employing a parameter homotopy is that after solving equations ( 31)-( 32) with a more general method a priori, we significantly reduce the computation required in solving equations ( 27)- (30) for each i.In practice, the entries of y are sampled uniformly along the complex unit circle.More details on how parameter homotopies fit into numerical algebraic geometry can be found in [1] and [2].
One additional reformulation that we can do to reduce computation is to define some of the variables intrinsically.This is common if one or more variables can be written as a linear combination of some of the other variables.Specifically, we know from Table 7 that: where cRAF tot is defined as a constant in Table 8.Thus, we can "remove" x 2 from our computations.Partial derivatives are no longer necessary with respect to the variable x 2 , and x 2 is no longer defined explicitly when tracking homotopy paths.
Table 9 summarizes the sequence of reductions made in the number of paths by imposing a {(a, x, y), λ}homogeneous structure followed by intrinsically defining the variable x 2 along with the number of paths required using the standard total degree homotopy [1], [2].13 record the distances between the data and model varieties for all 36 data points.A missing "interior" distance in Table 12 and Table 13 indicate there were no positive real critical points found for the given EGF level and replicate.However, we may still compute a distance to the boundary of the semi-algebraic sets corresponding to each model.Bertini input files, shell scripts, and MATLAB scripts are available within the supplementary files to analyze model selection and parameter estimation.The distances are summarized graphically in Figure S1.
Timing summaries for both the processive and distributive model can be found in Table 10.These timings include the numerical irreducible decomposition required to compute the dimension of each component in the model variety, solving equations ( 31)-( 32) required to employ a parameter homotopy scheme, and the parameter homotopy to solve equations ( 27)-( 30) for all i.Timings to compute the dimension of the model variety and the data solve were done in serial using a Apple MacBook Pro with 2.4 GHz Intel "Core i5" processor.The initial solve for the parameter homotopies were done in parallel using 96 (2.67 GHz Xeon-5650) compute nodes on the CentOS 5.11 operating system.
variety Intersection?Compute dimension of intersectionSelection requires more informationdim(V Mi ∩ V D )Find closest points between and(V Mi ∩ V D )Positive EmptyCompute distance d between andEx: Syn Bio (Could recover parameters) Ex: MAP kinase, HIV Is d less than tolerance?No, reject model Yes, model compatible Compare d for all models and select Ex: Syn Bio, MAP kinase

Figure S1 :
Figure S1: Schematic of numerical algebraic geometry framework corresponding to Algorithms I-III.(A) Input to algorithms include model (system of polynomials) translated into a model variety (red), and steadystate data translated into a data variety (blue).(B) Flow chart of model compatibility, parameter estimation and model selection methods.Examples (green) are described in Results section.

Figure
Figure S3: MAP Kinase model selection using 36 data points from 3 model variables (aggregate phosphoforms of ERK) taken in triplicate at each of 12 EGF stimulation levels.

Figure
Figure S1: Distance Plot three replicates per level Distance to model varieties Distance to processive model variety Distance to boundary of model varieties Distance to distributive model variety 1, Brent R. Davis 2 , Kenneth L. Ho 3 , Daniel J. Bates 2 , Heather A. Harrington 4 1 Department of Mathematics, San José State University, 2 Department of Mathematics, Colorado State University, 3 Department of Mathematics, Stanford University, 4 Mathematical Institute, University of Oxford.

Table 1 :
Expected timings collected over 20 runs.The table includes the average time and standard deviations associated to the four computations described in this section.

Table 5 :
Description of variables and parameters for distributive and processive MAP Kinase models

Table 6 :
Reaction velocities for the MAP Kinase distributive and processive model.The processive model uses the additional reaction velocities v 18 , v 19 , v 20 .

Table 7 :
Equations for distributive and processive MAP Kinase models

Table 8 :
Parameter values for MAP Kinase models

Table 10 :
Expected timings for the MAPK model collected over 20 'random' runs.

Table 12 and
Table

Table 12 :
Distance to (smaller) distributive model variety

Table 13 :
Distance to (larger) processive model variety