A parameter-free model discrimination criterion based on steady-state coplanarity
Heather A. Harrington, Kenneth L. Ho, Thomas Thorne, Michael P. H. Stumpf
aa r X i v : . [ q - b i o . Q M ] A ug A PARAMETER-FREE MODEL DISCRIMINATION CRITERION BASED ONSTEADY-STATE COPLANARITY
HEATHER A. HARRINGTON ∗ , KENNETH L. HO ∗ , THOMAS THORNE , AND MICHAEL P.H. STUMPF Abstract.
We describe a novel procedure for deciding when a mass-action model is incompatiblewith observed steady-state data that does not require any parameter estimation. Thus, we avoidthe difficulties of nonlinear optimization typically associated with methods based on parameterfitting. Instead, we borrow ideas from algebraic geometry to construct a transformation of the modelvariables such that any set of steady states of the model under that transformation lies on a commonplane, irrespective of the values of the model parameters. Model rejection can then be performedby assessing the degree to which the transformed data deviate from coplanarity. We demonstrateour method by applying it to models of multisite phosphorylation and cell death signaling. Ourframework offers parameter-free perspective on the statistical model selection problem, which cancomplement conventional statistical methods in certain classes of problems where inference has tobe based on steady-state data and the model structures allow for suitable algebraic relationshipsamong the steady state solutions. keywords: chemical reaction networks, Gr¨obner bases, mass-action kinetics, singular values,ordinary differential equations, algebraic statistics.many branches of science and engineering, one is often interested in the problem of model se-lection: given observed data and a set of candidate models for the process generating that data,which is the most appropriate model for that process? Such a situation commonly arises when theinner workings of a process are not completely understood, so that multiple models are consistentwith the current state of knowledge. For mechanistic models, e.g., ordinary differential equation(ODE) or stochastic dynamical models, most selection techniques involve parameter estimation,which typically requires some form of optimization, exploration of the parameter space, or formalinference procedure [1, 2]. For sufficiently complicated models, however, this task can become in-feasible, owing to the nonlinearity and multi-modality of the objective function (which penalizesany differences between the data and the model predictions), as well as the high dimensionality ofthe parameter space [3].Here, we present a framework for the discrimination of mass-action ODE models (and suitablegeneralizations thereof) that does not require or rely upon such estimated parameters. Our method(Fig. 1) operates on steady-state data and combines techniques from algebraic geometry, linearalgebra, and statistics to determine when a given model is incompatible with the data under all choices of the model parameters. The core idea is to use the model equations to construct atransformation of the original variables such that any set of steady states of the model under thattransformation possesses a simple geometric structure, irrespective of parameter values. In thiscase, we insist that the transformed steady states lie on a plane, which we detect numerically; if ∗ These authors contributed equally Division of Molecular Biosciences, Imperial College London, Wolfson Building, London, SW72AZ, UK Courant Institute of Mathematical Sciences and Program in Computational Biology, New YorkUniversity, 251 Mercer Street, New York, NY 10012, USA
E-mail address : [email protected], [email protected],[email protected], [email protected] . ata coplanar Data not coplanar
Data coplanar Data not coplanar
Models
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Observed Data (Steady state measurements)ˆ x . . .ˆ x . . .... . . .ˆ x m . . . Steady state invariantsModels
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariantsModels
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariantsModels
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariantsModels
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariants
Models
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Steady state invariantsCalculate elimination idealAssess coplanarity
Calculate elimination idealAssess coplanarityReduce number of variablesto include only observablesCharacterize steady states of modelsTransform model variables,parameters, and dataCalculate elimination idealAssess coplanarityReduce number of variablesto include only observablesCharacterize steady states of modelsTransform model variables,parameters, and dataCalculate elimination idealAssess coplanarityReduce number of variablesto include only observablesCharacterize steady states of modelsTransform model variables,parameters, and dataCalculate elimination idealAssess coplanarityReduce number of variablesto include only observablesCharacterize steady states of modelsTransform model variables,parameters, and data
Models
Model 1 Model 2 . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .
Steady state invariantsCalculate elimination idealAssess coplanarity
Calculate elimination idealAssess coplanarityReduce number of variablesto include only observablesCharacterize steady states of modelsTransform model variables, parameters, and dataparameters, and data
Data not coplanarModel compatibleModel incompatible Data not coplanarModel compatibleModel incompatible
Models
Model 1 . . . Model L. . . . . . . . .. . . . . . . . .. . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariantsModels
Model 1 . . . Model L. . . . . . . . .. . . . . . . . .. . . . . . . . .
Observed Data (Steady state measurements). . .. . .. . .. . .
Steady state invariants
Model
Model 1 . . . Model L. . . . . . . . .. . . . . . . . .. . . . . . . . .
Model (Polynomial di ff erential equations). . .. . .. . .. . . Model
Model 1 . . . Model L. . . . . . . . .. . . . . . . . .. . . . . . . . .
Model (Polynomial di erential equations)˙ x = . . .˙ x = . . .... . . .˙ x N = . . . Model
Model 1 . . . Model L. . . . . . . . .. . . . . . . . .. . . . . . . . .
Model (Polynomial di erential equations). . .. . .. . .. . .
Figure 1.
Parameter-free method for model discrimination (see text for details).the observed data are not coplanar under the transformation induced by a given model, then wecan confidently reject that model.The idea of transformation to coplanarity is not new, but previous efforts were limited, in part, byits systematic detection and quantification. For example, in [4], it was necessary to first manuallyreduce the dimension of the transformed space to three so that coplanarity could be assessedvisually. Other related research using similar methods include [5–7]. The current work extendsthis by devising a numerical scheme for quantifying the deviation from coplanarity that generalizesto higher dimensions and allows for statistical interpretation. Thus, we provide a richer and morepowerful framework for the application of this basic technique. Chemical reaction network theory(CRNT) [8,9] and stoichiometric network analysis [10] likewise embrace a parameter-free philosophyand can also be exploited for model selection [11–13].It is worth noting that our method provides a necessary but (generally) not sufficient conditionfor model compatibility: a model that is compatible with the data must provide a transformation tocoplanarity, but a model that achieves coplanarity is not necessarily compatible, due to additionaldegrees of freedom introduced in the transformation process. This is in contrast to traditional pproaches based on parameter fitting, which provide a sufficient but not necessary condition sincelocal extrema in the cost function surface may prevent a suitable fit. These two approaches aretherefore complementary and can be used together for improved model selection.The remainder of this paper is organized as follows. First, we introduce the concept of steady-state invariants [4, 5], polynomials that vanish at steady state and which depend only on exper-imentally accessible variables. Then we illustrate how to use steady-state invariants to deducecoplanarity requirements for model compatibility and how to detect such coplanarity numerically;we also discuss invariants in the context of standard parameter fitting techniques. Next, we applyour method to models of multisite phosphorylation and cell death signaling. Finally, we end withsome generalizations and concluding remarks.1. Steady-State Invariants
Consider a chemical reaction network model N X j =1 s ij X j k i −→ N X j =1 s ′ ij X j , i = 1 , . . . , R (1)in the species X , . . . , X N , where s ij and s ′ ij are the stoichiometric coefficients of X j in the reactantand product sets, respectively, of reaction i , with rate constant k i . Under mass-action kinetics, themodel has dynamics ˙ x i = R X j =1 k j (cid:0) s ′ ji − s ji (cid:1) N Y k =1 x s jk k , i = 1 , . . . , N, (2)where x i is the concentration of species X i (throughout, we follow the convention that lowercaseletters denote the concentrations of the corresponding species indicated in uppercase). Theseequations provide a quantitative description of the model and can, in principle, be used to testits validity by assessing the degree to which they are satisfied by observed data. Unfortunately, inpractice, the required variables are rarely all available. In particular, the velocities ˙ x = ( ˙ x , . . . , ˙ x N )can be difficult to measure, so we can often consider only the steady state ˙ x = , as we will do here.Furthermore, certain species may be experimentally inaccessible due to technological limitations;we eliminate these variables from the equations if possible.For simple models, this elimination can be done by hand, but a more systematic approach isrequired in general. One such approach is to use Gr¨obner bases [14], a central tool in computationalalgebraic geometry that provides a generalization of Gaussian elimination for multivariate polyno-mial systems. Here, we follow the general procedure of Manrai and Gunawardena [4]. Let Q [ a ] bethe polynomial ring consisting of all polynomials in the parameters a = ( k , . . . , k R ) with coeffi-cients from the rational numbers Q , and let K be its fraction field, comprising all elements of theform f /g , where f, g ∈ Q [ a ]. Clearly, each ˙ x i ∈ K [ x ], the ring of all polynomials in x = ( x , . . . , x N )with coefficients in K . Note that the parameters a have been absorbed into the coefficient field K ;thus, by performing all operations over K , we can treat a symbolically, i.e., without specifying anyparticular parameter values.To characterize the steady state ˙ x = , we construct the ideal J = h ˙ x i generated by ˙ x , consistingof all polynomials P Ni =1 f i ˙ x i , where each f i ∈ K [ x ]. Clearly, J contains all elements of K [ x ] thatvanish at steady state. To obtain only those elements of J that do not depend on the variables x , . . . , x i , we consider the i th elimination ideal J i = J ∩ K [ x obs ], where x obs = ( x i +1 , . . . , x N )denotes the “observable” variables. Here, it is useful to introduce Gr¨obner bases, which are specialsets of generators with the so-called elimination property that if g = ( g , . . . , g M ) is a Gr¨obner basisfor J under the lexicographic ordering x > · · · > x N , then J i = h g obs i , where g obs = g ∩ K [ x obs ]are precisely those elements of g containing only the variables x obs . The polynomials g obs generate ll elements of K [ x obs ] that vanish at steady state and so characterize the projection of the steadystate onto the variables x obs .Procedurally, we compute a reduced Gr¨obner basis g of J with respect to a suitable lexicograhicordering using standard algorithms, then obtain g obs by subselection. For numerical conveniencewe further rescale each polynomial in g obs so that all coefficients belong to Q [ a ] (i.e., we multiplythrough by their common denominator). Then the elements of g obs = ( I , . . . , I N inv ) have the form I i ( x obs ; a ) = n i X j =1 f ij ( a ) N obs Y k =1 x t ijk k , i = 1 , . . . , N inv , (3)where we have applied the relabeling x obs = ( x , . . . , x N obs ). Clearly, each I i is a polynomial in x obs that vanishes at steady state; we call such polynomials steady-state invariants (or sometimesjust invariants for short).Note in general that steady-state invariants may fail to exist since J i may be empty. Moreover,invariants and their properties (e.g., degrees) can depend delicately on the choice of monomialordering. Some manual intervention is therefore often required to obtain useful invariants. Wewill not treat this important (but subtle) issue here, instead focusing on the analysis of giveninvariants, however they are obtained. This also has the advantage of separating the computationof invariants from their interpretation, in principle allowing the use of invariants from varioustheories. Steady-state invariants, if they exist, describe relationships between observable variablesthat hold at steady state for any given realization of parameter values, regardless of other factorssuch as initial conditions.For full details on the computational procedure employed, see the accompanying Sage worksheet,which contains code for all computations performed ( Materials and Methods ). For further back-ground on algebraic geometry and Gr¨obner bases (including the potential problems of obtainingthem), see [14]; for other methods of variable elimination, see, e.g., [15]. Similar algebraic ideashave also appeared in the context of phylogenetics [16, 17].2.
Model Discrimination
We start with a set of steady-state measurements ˆ x obs ,i for i = 1 , . . . , m , and a given model withsteady-state invariants I = { I , . . . , I N inv } .2.1. Data Coplanarity.
An invariant, I ∈ I , can be written somewhat simplified as I ( x obs ; a ) = n X j =1 f j ( a ) N obs Y k =1 x t jk k . (4)We first describe a procedure for deciding whether it is possible that the invariant is compatiblewith the data, i.e., I (ˆ x obs ,i ; a ) = 0 , i = 1 , . . . , m, (5)for some choice of a . We therefore rewrite (4) as I ( y ; b ) = n X j =1 b j y j , (6)where y j = Q N obs k =1 x t jk k and b j = f j ( a ), with y = ( y , . . . , y n ) and b = ( b , . . . , b n ). Let ϕ bethe map taking x obs to y . Then compatibility implies that the transformed variable ˆ y = ϕ (ˆ x obs )corresponding to any observation ˆ x obs , considered as a point in R n with coordinates (ˆ y , . . . , ˆ y n ),lies on the hyperplane defined by the coefficients b . In other words, compatibility with the dataˆ x obs ,i implies that the corresponding transformed data ˆ y i = ϕ (ˆ x obs ,i ) are coplanar. n general, it is possible that the invariant vanishes trivially, ( b = ), under some choice ofparameters for which coplanarity need no longer hold. To discount this case, we can check, forinstance, that the denominator of the corresponding g obs is never zero. Then I always has at leastone nonzero coefficient; hereafter in this section, we assume that the invariant is non-vanishing inthis sense.Let Y ∈ R m × n be the matrix whose rows consist of the ˆ y i . Then the data are coplanar if andonly if Y b = for some nontrivial column vector b = . Such a vector, by definition, resides in thenull space of Y , which can be found using the singular value decomposition Y = U Σ V T , wherethe diagonal elements of Σ give the singular values σ i ≥ V . In particular, the smallest singular value σ min bounds the norm k Y b k for any b = via σ min = min k b k =1 k Y b k , (7)so if σ min >
0, then the data cannot be coplanar [18]. More generally, σ min gives the least squaresdeviation of the data from coplanarity under the scaling constraint k b k = 1. This quantity dependsonly on the data and is therefore parameter-free.Note that this applies for any choice of b , regardless of whether it can be realized by the originalparameters a . In this sense, the condition of small σ min provides a necessary but not sufficientcriterion for model compatibility. The additional degrees of freedom introduced by neglecting thefunctional forms f j effectively linearizes the compatibility condition (5), allowing for a simple, directsolution.To account for the presence of noise, suppose that we know each component ˆ x k of a measurementˆ x obs only up to an error ∆ˆ x k , with∆ˆ x k = ǫ ˆ x k Z, k = 1 , . . . , N obs , (8)where Z ∼ N (0 ,
1) is a standard normal random variable. We imagine that the noise parameter ǫ is given, for example, by instrument error. Then from the perturbation equation y + ∆ y = ϕ ( x obs + ∆ x obs ) , (9)we find, expanding to first order, that the error is propagated to the transformed variables as∆ˆ y = ∇ ϕ (ˆ x obs )∆ˆ x obs , where ∇ ϕ is the Jacobian of ϕ , with elements ( ∇ ϕ ) ij = ∂y i /∂x j . Therefore,∆ˆ y j = ǫ j Z, ǫ j = ǫ N obs X k =1 ( ∇ ϕ ) jk ˆ x k . (10)We now consider the effect of the ∆ˆ y i on σ min under the null hypothesis that the underlying ˆ y i are coplanar with coefficients b (of unit norm). For this, we study the vector Y b , whose entries areperturbed from zero to n X j =1 b j ∆ˆ y j = n X j =1 b j ǫ j Z (11)for each transformed datum ˆ y . Since k b k = 1 by assumption, if we rescale each row of Y by itscorresponding effective error ǫ eff = max j =1 ,...,n | ǫ j | ≥ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X j =1 b j ǫ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (12) hus obtaining Y ′ , then each entry of Y ′ b has the form µ i Z with | µ i | ≤
1, for i = 1 , . . . , m . Wehence define the coplanarity error ∆ = σ min (cid:0) Y ′ (cid:1) ≤ (cid:13)(cid:13) Y ′ b (cid:13)(cid:13) , (13)which, from the discussion above, is bounded by the length of a normal random vector with variances µ i ≤
1, whose distribution function clearly dominates that of the length of a normal random vectorwith variances µ i = 1. But this latter quantity simply follows the χ distribution with m degrees offreedom. In other words, Pr (∆ ≥ x ) ≤ Pr ( X ≥ x ) , X ∼ χ m ;(14)if p α is the upper α -percentile for χ m (e.g., α = 0 . ≥ p α ) ≤ Pr ( X ≥ p α ) = α, (15)which gives an approximate criterion for rejecting coplanarity. As the amount of data increases,the approximation improves since σ min ( Y ′ ) → k Y ′ b k as m → ∞ by the symmetry of (10).Depending on the exact situation at hand, it may be appropriate to choose a more conservativesignificance level α or to invoke additional criteria in order to decide whether a model is accept-able. In the examples below, however, we will see that whether a model can be rejected is oftenfairly obvious, and in such cases we will simply use the asymptotic arguments based on the χ m distribution.2.2. Invariant Minimization.
Steady-state invariants can also be used in conjunction with stan-dard parameter fitting techniques. The basic approach is to minimize the Frobenius norm of thematrix θ ∈ R N inv × m , with entries θ ij = I i (ˆ x obs ,j ; a ), over the parameters a . This readily provides asufficient condition for model compatibility since any a producing a small norm provides parametersthat fit the data by construction. However, the condition is not necessary since suitable parametersmay fail to be found even for compatible models due to the intricacies of the objective function.Clearly, prior knowledge of a can be used to guide the optimization away from such difficulties.Assuming that the model and its parameters are correct, each invariant I (ˆ x obs ; a ) = 0 in prin-ciple. However, due to noise, I (ˆ x obs ; a ) = ǫ (ˆ x obs ; a ) Z , where ǫ (ˆ x obs ; a ) = ǫ n X j =1 f j ( a ) N obs X k =1 ( ∇ ϕ ) jk ˆ x k (16)by (11). Therefore, if we use I (ˆ x obs ; a ) /ǫ (ˆ x obs ; a ) as the entry of θ corresponding to invariant I and datum ˆ x obs , then the invariant error θ ( a ) = k θ ( a ) k F ∼ χ N inv m . (17)This can be used to compute the likelihood L ( a ) = Pr( θ ( a )) and allows, e.g., various likelihood-based selection schemes [19, 20], assuming that the optimization can be performed. Here, we usethe Akaike information criterion (AIC), A = 2 R − L max , (18)where L max = max a L ( a ), which penalizes model complexity; the preferred model is the one withthe minimum AIC [21]. 3. Results
We apply our methods to two illustrative biological processes for which competing models exist:multisite phosphorylation and cell death signaling. .1. Multisite Phosphorylation.
We focus first on phosphorylation, a key cellular regulatorymechanism that has been the subject of extensive study, both experimentally [22–24] and theoret-ically [4, 5, 25–27]. Following [4], we consider a two-site system with reactions, K + S u a u −− ⇀↽ −− b u KS u c uv −−→ K + S v , (19a) F + S v α v −− ⇀↽ −− β v F S v γ vu −−→ F + S u , (19b)where u, v ∈ { , } are bit strings of length two, encoding the occupancies of each site (0 or 1 forthe absence or presence, respectively, of a phosphate), with u having less bits than v ; S u is thephosphoform with phosphorylation state u ; K is a kinase, an enzyme that adds phosphates; and F is a phosphatase, an enzyme that removes phosphates. Each enzyme can be either processive (P),where more than one phosphate modification may be achieved in a single step, or distributive (D),where only one modification is allowed before the enzyme dissociates from the substrate ( c = 0for K , γ for F ). This mechanistic diversity generates four competing models: PP, PD, DP,and DD; where the first letter designates the mechanism of the kinase, and the second, that of thephosphatase.As in [4], we consider only the concentrations x obs = ( s , s , s , s ) as observable and use theordering, ( ks , ks , ks , f s , f s , f s , k, f, s , s , s , s ) , (20)with which we are able to eliminate all other variables except f from the dynamics of each model.The remaining Gr¨obner basis polynomials are of the form p ( f, x obs ) = f · q ( x obs ), where f = 0unless there is no phosphatase in the system, which we assume not to be the case, so we take onlythe observable part q ( x obs ). It is easy to check that the resulting denominators are always of onesign.Each model has three steady-state invariants. Matched appropriately, the invariants for modelPP share the same transformed variables y = ϕ ( x obs ) as those for PD; the same is true for DP andDD. Thus, in terms of the transformed data, only the kinase mechanism is discriminative. BetweenPP/PD and DP/DD, two invariants ( I and I ) are discriminative in principle, though only one( I ) succeeds numerically: for simulated data from the PP/PD models, provided that the noiselevel is sufficiently low, lack of coplanarity on I is able to correctly reject the DP/DD models atsignificance level α = 0 .
05 (∆ ∼ versus ∆ ∼ ǫ = 10 − , against a thresholdof p α = 11 . I ,which has transformed variables, y PP / PD = (cid:0) s s , s s , s s , s s , s , s s (cid:1) , (21a) y DP / DD = (cid:0) s s , s s , s s , s , s s (cid:1) (21b)for PP/PD and DP/DD, respectively, i.e., y PP / PD has the additional variable s s over y DP / DD .Therefore, PP/PD models can be made to fit DP/DD data simply by setting the coefficient corre-sponding to s s to zero, which is in fact what we observe. No model is rejected on the basis ofdata generated from it.We emphasize that these results are specific to the particular ordering chosen. Indeed, one canmake the phosphatase mechanism discriminative instead by reversing the order of the variables x obs in (20). The exhaustive analysis of such orderings is beyond the scope here; rather, we aim toillustrate the potential uses (and usefulness) of this type of approach using concrete examples.Although the condition of coplanarity is technically valid only at steady state, there should never-theless be some convergence over time to coplanarity for any compatible model. We hence compute∆ for the PP/PD and DP/DD models along time course trajectories simulated from model PP at arious levels of ǫ (Fig. 2A). For low noise, the results confirm convergence for invariants previouslyidentified as compatible (all I i for PP/PD; I and I for DP/DD), with stagnation for incompatibleinvariants ( I for DP/DD); this does suggest wider applicability of this method, provided that thedata are approaching steady state reasonably fast. As the noise increases, however, ∆ decreasesinversely proportionally, until the stagnation point hits the basal error level of ∆ ∼ ǫ ∼ − (Fig. 2B).To further discriminate between all four models we next turn to invariant minimization. Therequired optimization involves highly nonlinear functions, so success should be expected only if wehave good initial estimates of the model parameters. This is a rather strong demand. In such acase, however, minimization is indeed capable of identifying the correct model from the data solong as ǫ . − (Fig. 2C). These results reinforce our belief that the algebraic approach proposedhere naturally complements conventional (i.e. parametric) reverse engineering schemes such asoptimization or inference procedures.3.2. Cell Death Signaling.
We next apply our methods to receptor-mediated cell death signaling,the so-called extrinsic apoptosis pathway, which plays a prominent role in cancers and other diseases[28–31]. Specifically, we consider the assembly of the death-inducing signaling complex (DISC), amulti-protein oligomer formed by the association of FasL, a death ligand, with its cognate receptorFas [32, 33].We investigate two models of DISC formation. The first [34], which we call the crosslinkingmodel , is based on the successive binding of Fas ( R ) to FasL ( L ), L + R k f −− ⇀↽ −− k r C , (22a) C + R k f −− ⇀↽ −− k r C , (22b) C + R k f −− ⇀↽ −− k r C , (22c)where C i is the complex FasL:Fas i . The second [6], which we call the cluster model , posits threeforms of Fas (inactive, X ; active and unstable, Y ; active and stable, Z ) and specifies receptorcluster-stabilization events driven by FasL, X k o −− ⇀↽ −− k c Y, (23a) Z k u −→ Y, (23b) jY + ( i − j ) Z k ( i ) s −−→ ( j − k ) Y + ( i − j + k ) Z, (23c) L + jY + ( i − j ) Z k ( i ) l −−→ L + ( j − k ) Y + ( i − j + k ) Z, (23d)where the last two reactions represent entire families generated by taking i = 2 or 3, with j = 1 , . . . , i and k = 1 , . . . , j . The cluster model is capable of bistability, whereas the crosslinking model exhibitsonly monostable behavior [6].The two models are structurally very different, and discriminating between them requires somecare. Hence, following [6], we establish a correspondence between the models by considering theapoptotic signal ζ transduced by the DISC, defined as ζ = c + 2 c + 3 c for the crosslinkingmodel and ζ = z for the cluster model. We assume that ζ is experimentally accessible; othervariables assumed accessible include λ , the total concentration of FasL ( λ = l + c + c + c and λ = l for the crosslinking and cluster models, respectively), and ρ , the total concentration of Fas( ρ = r + c +2 c +3 c and ρ = x + y + z , respectively). Eliminating all other variables via the orderings BC Figure 2.
Discrimination of multisite phosphorylation models. ( A ) Coplanarityerror ∆ of the steady-state invariants of the PP/PD (left) and DP/DD (right) modelsalong time course trajectories simulated from the PP model, corrupted by variouslevels of noise (lined, ǫ = 10 − ; dashed, ǫ = 10 − ; dotted, ǫ = 10 − ). At each noiselevel, the errors for three invariants are shown (blue, I ; green, I ; red, I ). ( B )Coplanarity error ∆ of DP/DD invariants on PP data at steady state as a functionof the noise level ǫ ; invariants colored as in (A). The shaded region indicates theregime over which the DP/DD models can be rejected at significance level α = 0 . C ) Invariant error AIC A for each model (blue, PP; green, PD; red, DP; cyan, DD)on data generated from the PP (top left), PD (top right), DP (bottom left), andDD (bottom right) models.( c , c , λ, ρ, ζ ) and ( y, λ, ρ, ζ ) for the crosslinking and cluster models (after appropriate variablesubstitutions) we obtain one non-vanishing steady-state invariant for each model. The dimensionsof the transformed spaces are 5 and 15 for the crosslinking and cluster models, respectively.As for phosphorylation, we compute the coplanarity error for each invariant on time course datasimulated from each model at various noise levels. Although results are inconclusive for data from igure 3. Discrimination of cell death signaling models. ( A ) Coplanarity error ∆ ofthe steady-state invariants of the crosslinking (left) and cluster (right) models alongtime course trajectories simulated from the cluster model, corrupted by various levelsof noise (blue, ǫ = 10 − ; green, ǫ = 10 − ; red, ǫ = 10 − ). ( B ) Coplanarity error ∆ ofmodel invariants (blue, crosslinking; green, cluster) on cluster data at steady stateas a function of the noise level ǫ . The shaded region indicates the regime over whichthe crosslinking model can be rejected at significance level α = 0 .
05. ( C ) Invarianterror AIC A for each model (blue, crosslinking; green, cluster) on data generatedfrom the crosslinking (left) and cluster (right) models.the crosslinking model, the coplanarity criterion can reject the crosslinking model on the basis ofcluster model data at α = 0 .
05, provided that ǫ . − (Fig. 3A and B). The minimization protocolalso correctly identifies the model from the data over the same range of noise levels (Fig. 3C).4. Discussion
In this paper, we have presented a novel model discrimination scheme based on steady-statecoplanarity that does not require known or estimated parameter values. Thus we are able tosidestep the parameter inference problem common to many fields including systems biology [3, 35].Such algebraic methods are not always effective, however: steady-state invariants may not exist,and even when they do, the additional degrees of freedom introduced by effective linearizationcan cause the method to fail. A promising solution to the problem when invariants cannot be alculated using Gr¨obner bases may be to employ invariants from CRNT [36]. Our results alsosuggest a somewhat low tolerance for noise, which can restrict its applicability. Significantly, ourmethod has the unique feature that it can be applied with complete ignorance of parameter values,and is therefore a useful additional tool in the analysis of inverse problems involving dynamicalsystems.Rather than competing directly with current model discrimination techniques, we expect thatcoplanarity will form one end of an entire spectrum of methods, to be used when no parameterinformation is available. At the other end lie methods based on parameter estimation (includinginvariant minimization), which, for dynamical systems, can depend delicately on qualitative andquantitive aspects of the systems under consideration [37, 38]. The intermediate regime comprisestechniques that can leverage partial knowledge, for instance, constraints on certain parameter valuesor qualitative features of the dynamics [39]. Along this spectrum, naturally, the discriminativepower increases with the amount of prior information available. In this broader context, coplanaritycan be used to efficiently reject candidate models before employing more demanding parameterestimation tools. Thus, it can serve as a preprocessor to thin out the model space. The realadvantages and limitations of any inferential procedure become apparent once their performancecan be evaluated in real-world applications. This is perhaps particularly true for this currentapproach. Certainly a range of theoretical and computational issues surround algebraic methodswhich will likely impact their applicability. Here we have found that a pragmatic approach yieldssome useful insights for small and intermediate-sized problems.Finally, we remark that the presented scheme is but the simplest of a potential new class ofparameter-free selection methods based on the detection of geometric structure. In this view,transformation to coplanarity is just one of many low-dimensional descriptions of such structure.The existence of low-dimensional representations has recently been predicted in neuronal signaling[40], and can ultimately be attributed to the inherent robustness of biological systems [41, 42].5. Materials
Gr¨obner Basis Calculation.
All reduced Gr¨obner bases are computed over the field K of ra-tional functions in the parameters a with rational coefficients, under a suitable lexicographic order-ing with the observables x obs located at the end of the variable list, using the computer algebra sys-tem Singular
Data Generation.
For each model parameters are drawn independently from a log-normaldistribution with median µ ∗ = e µ = 1 and multiplicative standard deviation σ ∗ = e σ = 2, where µ and σ are the mean and standard deviation, respectively, of the underlying normal distribution.Using these parameters m = 100 time course trajectories are computed for each model via integra-tion of the model ODEs over the time interval 0 ≤ t ≤ µ ∗ = 1 and σ ∗ ǫ = 10 − to 10 − , for each ǫ , multiplying the nominaldata by random log-normal samples with µ ∗ = 1 and σ ∗ = 1 + ǫ .5.3. Invariant Minimization.
Invariant error likelihood maximization is performed in two phases.First, an approximate optimal parameter set is obtained by minimizing the Frobenius norm of thematrix η ∈ R N inv × m , where each entry corresponds to an invariant-datum pair as in θ , but withvalues I (ˆ x obs ; a ) /M (ˆ x obs ; a ), where M (ˆ x obs ; a ) = n X j =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f j ( a ) N obs Y k =1 ˆ x t jk k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (24) his is then taken as an initial parameter estimate to compute L max . All optimizations are per-formed using L-BFGS-B [43] through SciPy, with lower and upper bounds of 0 .
01 and 100, respec-tively, for each variable. The minimization of k η k F is seeded with initial value 1 for all variables.5.4. Computational Platform.
Acknowledgments
HAH and MPHS gratefully acknowledge the Leverhulme Trust. KLH was funded in part byNational Science Foundation grant DGE-0333389. TT and MPHS also acknowledge funding fromthe Biotechnology and Biological Research Council; MPHS is a Royal Society Wolfson ResearchMerit award holder. We thank Carsten Wiuf, Elisenda Feliu, and Sarah Filippi for their commentson the manuscript.
References [1] Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MPH (2009) Approximate Bayesian computation scheme forparameter inference and model selection in dynamical systems.
J R Soc Interface
Bioinformatics
FEBS J
Biophys J
Biophys J
PLoS Comput Biol
PLoS ONE
Chem Eng Sci
Chem Eng Sci
Cell Biophys
IEE Proc Syst Biol
J Mol Catal A: Chemical
J Mol Catal A: Chemical
Ideals, Varieties, and Algorithms (Springer, New York).[15] Feliu E, Wiuf C (2012) Variable elimination in chemical reaction networks with mass-action kinetics.
SIAM JAppl Math
J Classif
SIAM Rev
Regression Diagnostics: Identifying Influential Data and Sources ofCollinearity (John Wiley & Sons, Inc., Hoboken, NJ).[19] Casella G, Berger RL (2001)
Statistical Inference (Duxbury Press).[20] Schwarz G (1978) Estimating the dimension of a model.
Ann Statist
IEEE Trans Automat Control
Nature
TrendsBiochem Sci
24] Seger R, Krebs EG (1995) The MAPK signaling cascade.
FASEB J
Proc NatlAcad Sci USA
Proc NatlAcad Sci USA
Nature
Science
Nature
Nature
Oncogene
Science
Cell Death Differ
Math Biosci Eng
Elements of Computational Systems Biology , eds Lodhi HM, Muggleton SH (JohnWiley & Sons, Inc., Hoboken, NJ), pp 21–48.[36] Karp RL, P´erez Mill´an M, Dasgupta T, Dickenstein A, Gunawardena J (2012) Complex-linear invariants ofbiochemical networks.
J Theor Biol , in press.[37] Erguler K, Stumpf MPH (2011) Practical limits for reverse engineering of dynamical systems: a statisticalanalysis of sensitivity and parameter inferability in systems biology models.
Mol BioSyst
PLoS Comput Biol
BMC Syst Biol
Proc Natl Acad Sci USA
Science
Nat Rev Genet
ACM Trans Math Softw23:550–560.