A Case Study on the Parametric Occurrence of Multiple Steady States
Russell Bradford, James H. Davenport, Matthew England, Hassan Errami, Vladimir Gerdt, Dima Grigoriev, Charles Hoyt, Marek Kosta, Ovidiu Radulescu, Thomas Sturm, Andreas Weber
AA Case Study on the Parametric Occurrence ofMultiple Steady States
Russell Bradford
University of Bath, U.K.
[email protected] James H.Davenport
University of Bath, U.K.
[email protected] England
Coventry University, U.K.
[email protected] Hassan Errami
University of Bonn, Germany [email protected] Gerdt
JINR, Dubna, Russia [email protected] Dima Grigoriev
CNRS & University of Lille, France [email protected] Hoyt b-it, Bonn, Germany [email protected] Marek Košta
Slovak Academy of Sciences [email protected] Radulescu
DIMNP UMR CNRS/UM 5235,University of Montpellier, France [email protected] Thomas Sturm
U Lorraine, CNRS, Inria & LORIA, Nancy, FranceMPI Informatics & Saarland University, Germany [email protected] Weber
University of Bonn, Germany [email protected]
ABSTRACT
We consider the problem of determining multiple steadystates for positive real values in models of biological net-works. Investigating the potential for these in models ofthe mitogen-activated protein kinases (MAPK) network hasconsumed considerable effort using special insights into thestructure of corresponding models. Here we apply combi-nations of symbolic computation methods for mixed equal-ity/inequality systems, specifically virtual substitution, lazyreal triangularization and cylindrical algebraic decomposi-tion. We determine multistationarity of an 11-dimensionalMAPK network when numeric values are known for all butpotentially one parameter. More precisely, our consideredmodel has 11 equations in 11 variables and 19 parameters, 3of which are of interest for symbolic treatment, and further-more positivity conditions on all variables and parameters.
1. INTRODUCTION
The occurrence of multiple steady states (multistationar-ity) has important consequences on the capacity of signaling . pathways to process biological signals, even in its elemen-tary form of two stable steady states (bistability). Bistableswitches can act as memory circuits storing the informa-tion needed for later stages of processing [26]. The responseof bistable signaling pathways shows hysteresis, namely dy-namic and static lags between input and output. Becauseof hysteresis one can have, at the same time, a sharp binaryresponse and protection against chatter noise.Bistability of signaling usually occurs as a result ofactivation of upstream signaling proteins by downstreamcomponents [3]. A different mechanism for producingbistability in signaling pathways was proposed by Marke-vich et al. [19]. In this mechanism bistability can becaused by multiple phosphorylation/dephosphorylation cy-cles that share enzymes. A simple, two-step phosphory-lation/dephosphorylation cycle is capable of ultrasensitiv-ity, a form of all or nothing response with no hysteresis(Goldbeter–Koshland mechanism). In multiple phosphory-lation/dephosphorylation cycles, enzyme sharing providescompetitive interactions and positive feedback that ulti-mately leads to bistability.Algorithmically the task is to find the positive real so-lutions of a parameterized system of polynomial or ratio-nal systems, since the dynamics of the network is given bypolynomial systems (arising from mass action kinetics) orrational functions (arising in signaling networks when someintermediates of the reaction mechanisms are reduced). Dueto the high computational complexity of this task [13] con-siderable work has been done to use specific properties ofnetworks and to investigate the potential of multistationar- a r X i v : . [ c s . S C ] A p r ty of a biological network out of the network structure. Thisonly determines whether or not there exist rate constants al-lowing multiple steady states, instead of coming up with asemi-algebraic description of the range of parameters yield-ing this property. These approaches can be traced back tothe origins of Feinberg’s Chemical Reaction Network Theory (CRNT) whose main result is that networks of deficiency 0have a unique positive steady state for all rate constants [12,8]. We refer to [6, 20, 16] for the use of CRNT and othergraph theoretic methods to determine potential existence ofmultiple positive steady states, with [17] giving a survey.Given a bistable mechanism it is also important to com-pute the bistability domains in parameter space, namelythe parameter values for which there is more than one sta-ble steady state. The size of bistability domains gives thespread of the hysteresis and quantifies the robustness of theswitches. The work of Wang and Xia [24] is relevant here:they used symbolic computation tools, including cylindricalalgebraic decomposition as we do below, to determine thenumber of steady states and their stability for several sys-tems. They reported results up to a 5-dimensional systemusing specified parameter values, but their method is exten-sible to parametric questions. Higher-dimensional systemswere studied using sign conditions on the coefficients of thecharacteristic polynomial of the Jacobian. In some casesthese guarantee uniqueness of the steady state [7].In this paper we use an 11-dimensional model of amitogen-activated protein kinases (MAPK) cascade [19] asa case study to investigate properties of the system usingalgorithmic methods towards the goal of semi-algebraic de-scriptions of parameter regions for which multiple positivesteady states exist.
2. THE MAPK NETWORK
The model of the MAPK cascade we are investigating canbe found in the Biomodels database [18]. We have renamedthe species names to x , . . . , x and the rate constants to k , . . . , k to facilitate reading:˙ x = k x + k x − k x x − k x x ˙ x = k x + k x + k x + k x − x x ( k + k ) − k x x ˙ x = k x + k x − k x x ˙ x = x ( k + k ) + x ( k + k ) − k x x − k x x ˙ x = k x + k x + k x + k x − x x ( k + k ) − k x x − k x x ˙ x = k x x − x ( k + k )˙ x = k x x − x ( k + k )˙ x = k x x − x ( k + k )˙ x = k x − k x + k x x ˙ x = k x x − x ( k + k )˙ x = k x − k x + k x x . (1)The Biomodels database also gives us meaningful values for k = 0 . , k = 1 , k = 0 . , k = 0 . ,k = 1 , k = 15 , k = 0 . , k = 1 ,k = 0 . , k = 1 , k = 0 . , k = 0 . ,k = 1 , k = 0 . , k = 0 . , k = 0 . . (2)Some of these values are measured and some are well-educated guesses. For the purpose of our study we assumethey are suitable.We add three linear conservation constraints introducingthree further constant parameters k , k , k : x + x + x + x + x = k x + x + x = k x + x + x + x + x + x + x + x + x = k . (3)Computations to produce these in MathWorks SimBiologyuse the left-null space of the stoichiometric matrix underpositivity conditions, see for example [21].Meaningful values for k , k , k are harder to obtainthan the constants in (2). The following are some realisticvalues estimated by ourselves on the basis of our understand-ing of the biological model: k = 100 , k = 50 , k ∈ { , } . (4)The long-term goal of our research is to treat all three ofthese together parametrically, although in the present workwe focus on situations with one free-parameter.The steady state problem for the MAPK cascade can nowbe formulated as a real algebraic problem. That is, we re-place the left hand sides of all equations in (1) with 0. Thistogether with the equations in (3) yields an algebraic systemwith polynomials in F ⊂ Z [ k , . . . , k ][ x , . . . , x ] . All entities in our model are strictly positive, which yieldsin addition a system P = { k , . . . , k , x , . . . , x } ⊂ Z [ k , . . . , k ][ x , . . . , x ]establishing a side condition on the solutions of F that allvariables x i and parameters k i of P be positive. In termsof first-order logic our specification of F and P yields aquantifier-free Tarski formula ϕ = (cid:94) f ∈ F f = 0 ∧ (cid:94) v ∈ P v > . (5)The estimations for the rate constants in (2) formally estab-lish a substitution rule σ = [0 . /k , . . . , . /k ], whichcan be applied to F , P , or ϕ in postfix notation. In this section we are going to analyze the system formultiple positive steady states. As we will not include apriori information about the stability of the fixed points,we do not only have to consider (at least) two stable fixedpoints but also unstable fixed points, i.e., we investigate theexistence of at least three different roots x ∈ ]0 , ∞ [ of F for given choices k ∈ ]0 , ∞ [ of parameters.We present two investigations: one using the Redlog pack-age in Reduce and the other using the Regular Chains Li-brary in Maple. Both will make use of Cylindrical Algebraicecomposition (CAD) [1] to solve the problem. The worst-case time complexity of CAD is doubly exponential. Ourapproaches admit, in principle, arbitrary numbers of inde-terminates. However, for the sake of realistic computationtimes we must restrict ourselves to one free parameter. Eventhen, the number of variables present is too large for con-temporary CAD implementations. We make progress bycombining CAD with additional symbolic methods. Ourfirst approach uses virtual substitution techniques and thesecond real triangularization. In both cases we have com-bined the corresponding methods by hand, but automationis clearly possible.
Real Quantifier Elimination (QE) can directly handlethe parametric existence of steady states, taking as input ∃ x . . . ∃ x ϕ , possibly with substitutions for some param-eters. However, we are not only interested in the existencebut also in the number of solutions. We are going to com-bine Virtual Substitution (VS) [25] with CAD. The formersmoothly eliminates the majority of the quantifiers whilethe latter allows us to count numbers of solutions via de-composition of the remaining low-dimensional spaces. Thatcombination of methods requires the solution of several QEruns with each problem and some combinatorial arguments.Throughout this subsection we are using Redlog [9]. Parameter-Free Computations.
We consider ϕ = ϕσ [100 /k , /k , /k ] whereall parameters have been substituted with rational numbers.The closed formula ¯ ϕ = ∃ x . . . ∃ x ϕ states the exis-tence of a suitable real solution. In a first step, we solve for i ∈ { , . . . , } the following eleven QE problems using VS: ϕ ( i )500 = VS( ∃ x . . . ∃ x i − ∃ x i +1 . . . ∃ x ϕ ) . Each ϕ ( i )500 is a univariate quantifier-free formula describingall possible real choices for x i for which there exist realchoices for all other variables such that ϕ holds. CAD caneasily decompose the corresponding one-dimensional spaces.It turns out that for each x i there are exactly three zero-dimensional cells a i , b i , c i ∈ R for which ϕ ( i )500 holds. Weextract all a i , b i , and c i as real algebraic numbers , i.e., uni-variate defining polynomials with integer coefficients plusisolating intervals. By combinatorial arguments it is nothard to see that the following holds for the set S of realsolutions of ϕ :3 ≤ | S | and S ⊆ (cid:81) i =1 { a i , b i , c i } . Notice that at this point we have proven multistationar-ity for k = 500. We can furthermore compute S byplugging the 3 candidates from the Cartesian product into ϕ . A straightforward approach requires arithmetic withreal algebraic numbers followed by the determination of thesigns of the results, which is quite inefficient in practice. Weuse instead a heuristic approach combining refinements ofthe isolating intervals of the real algebraic numbers with in-terval arithmetic. This excludes 3 − Traditionally, doubly exponential in the number of vari-ables. However recent progress on CAD in the presence ofequational constraints [10], such as (1) with 0 for left-handside, allows us to conclude it is actually doubly-exponentialthe number of variables minus the number of equational con-straints at different levels of the projection [11]. solutions. The three remaining candidates require no fur-ther checking since we already know that | S | ≥
3. Theoverall CPU time is 71.3 seconds for 11 runs of VS plus 11runs of CAD, followed by 16 hours for checking candidates. Our checking procedure is a file-based prototype starting aReduce process for every single of the 3 candidates; thereis considerable room for optimization.For k = 200 instead of 500 all eleven univariate CADcomputations yield unique solutions which can be straight-forwardly combined to one unique solution for the corre-sponding ϕ . The overall CPU time here is 66.4 secondsfor 11 runs of VS plus 11 runs of CAD. Machine float ap-proximations of all our solutions are given in Table 1. Parametric Analysis for k . We now consider ϕ k = ϕσ [100 /k , /k ] leaving k as a parameter. Again, we solve for i ∈ { , . . . , } elevenQE problems using VS: ϕ ( i ) k = VS( ∃ x . . . ∃ x i − ∃ x i +1 . . . ∃ x ϕ k ) . This time each ϕ ( i ) k is a bivariate quantifier-free formulain k and the corresponding x i . This time we construct atwo-dimensional CAD for each ϕ ( i ) k . The projection orderis important: we first project x i , then the CAD base phasedecomposes the k -axis, followed by an extension phase thatdecomposes the x i -space over the k -cells obtained in thebase phase. This is feasible with one limitation: we do notextend over zero-dimensional k -cells. In other words, weaccept finitely many blind spots in parameter space, whichwe can explicitly read off from the CAD so that in the endwe know exactly what we are missing.Figure 1 shows our CAD tree for ϕ (2) k . The first layer nextto the root shows the decomposition of the k -axis. Thefive zero-dimensional (rectangular) cells are the previouslymentioned blind spots, among which the smallest one withnegative value of k is not relevant. Those zero-dimensionalcells also establish the limits of the full dimensional (oval)cells in between. The cylinders over those one-dimensional k -cells each contain either one or three zero-dimensional x -cells where ϕ (2) k holds. We have deleted from the tree all x -cells where ϕ (2) k does not hold. We make two observa-tions, important for a qualitative analysis of our system:(i) For all positive choices of k —extending to infinity—there is at least one positive solution for x .(ii) There is a break point around k = 409 .
253 where thesystem changes from unique solutions to exactly threesolutions.Recall that for all floating point numbers given here as ap-proximations we in fact know exact real algebraic numbers.For instance, the exact break point is the only real zero inthe interval (409 , (cid:80) i =0 c i k i with integer coefficients c i as in Table 2 . (6)Figure 2 depicts all eleven CAD trees for ψ (1) k , . . . , ψ (11) k .They are quite similar to the one just discussed. Even thebreak point from one to three solutions for x i is identical forall i ∈ { , . . . , } so that we can generalize our observations: All QE-related computations have been carried out on a 2.4GHz Intel Core i7 with 3 GB RAM or cores on a computeserver with similar speed and memory limitations. able 1: The unique solution x (200) for k = 200 and the three solutions x (500)1 , . . . , x (500)3 for k = 500 . We haveactually computed real algebraic numbers, which are pairs of univariate polynomials and isolated intervals.For convenience we are giving machine float approximations here, which can be made arbitrarily precise. x (200) = (90 . , . , . , . , . , . , . , . , . , . , . x (500)1 = (17 . , . , . , . , . , . , . , . , . , . , . x (500)2 = (122 . , . , . , . , . , . , . , . , . , . , . x (500)3 = (323 . , . , . , . , . , . , . , . , . , . , . Figure 1: The pruned CAD tree for x . Ellipses and rectangles are full-dimensional and zero-dimensionalcells, respectively. We have removed cells where k is negative or where the input formula is false. (i) For all positive choices of k —extending to infinity—there is at least one positive solution for ( x , . . . , x ).(ii) There is a break point β around k = 409 .
253 wherethe system changes its qualitative behavior. We haveexactly given β as a real algebraic number in Equa-tion (6). For k < β there is exactly one positive so-lution for ( x , . . . , x ). For k > β there are at least3 and at most 3 positive solutions for ( x , . . . , x ).The overall computation time for our parametric analysis is4.3 minutes. It is strongly dominated by 2.8 minutes for thecomputation of one particular CAD tree, for ϕ (11) k . It turnsout that the suitable projection order with x i eliminated firstis computationally considerably harder than projecting theother way round. As a preprocessing step we apply CAD-based simplification of the ϕ ( i ) k with the opposite, faster,projection order. Here we use Qepcad B, which performsbetter than Redlog at simple solution formula construction. We now describe an alternative approach to the solutionusing regular chains methods. Regular chains are the tri-angular decompositions of systems of polynomial equations(triangular in terms of the variables in each polynomial).Highly efficient methods for working in complex space havebeen developed based on these; see [23] for a survey.Recent work by Chen et al. [4] proposes adaptations ofthese tools to the real analogue: semi-algebraic systems. They describe two algorithms to decompose any real poly-nomial system into finitely many regular semi-algebraic sys-tems. The first does so directly while the second, Lazy RealTriangularize (LRT) produces the highest dimension solu-tion component and unevaluated function calls, which if allevaluated would combine to give the full solution. Thesealgorithms are implemented in the Regular Chains Library in Maple which we use throughout this subsection.We apply LRT on the quantifier-free formula (5) evaluatedwith the parameter estimates for k , . . . , k given at thestart of Section 2, so we have one free parameter as in theprevious section. We need to choose a variable ordering: ouranalysis requires that k be the indeterminate consideredalone; the remaining variables are placed in lexicographi-cal order (the in-built heuristics to make the choice couldsuggest nothing better). The solutions must hence containconstraints in k , constraints in ( x , k ), in ( x , x , k )and so on. We define the main variable of a constraint tobe the highest one present in this ordering.LRT produces one solution component and 6 unevaluatedfunction calls in less than 3 seconds. In the evaluated com-ponent: for each of x , . . . , x there is a single equationwhich had this as the main variable. Further, these are alllinear in their main variable meaning they can be easily re-arranged into the solution formulae in Table 3.The constraints on ( x , k ) are that x > igure 2: All CAD trees for ψ (1) k , . . . , ψ (11) k . For positive k there are always either one or three positivesolutions for the corresponding x i . The break point from one to three solutions is the same in all trees. In thesecond but last row on the left hand side there is the tree for ψ (1) k , which is shown in more detail in Figure 1.able 2: Coefficients c i and d j of polynomials occurring in Equations (6) and (7), respectively. c = 351590934502740290936895033267017158736060313940693076650155371250411 c = − c = 25374851641220554774259605635053469432582109883965015804077119110958034090 c = 12972493018300022707027639267804259251235991618029852880330004508564391594000 c = − c = 2231098270337406450670301663172664333421440833875848621423683265663846533079600000 c = − c = 39262101548790869407057994985320156500968958361396178908180026842806643766783104000000 c = − c = 70978850735887473459176997186175978425873267246760023212940616924643171868478080000000000 c = − d = 16838105723097694257603469 d = − k + 7723967969644977896148686580 d = 8176202638735769127032169 k − k + 1232154357941338876156606812900 d = 1465408757440589841803452380 k − k + 83152655240002767729550477640000 d = 85462524901276846107251669400 k − k + 2556805354853318332197489636000000 d = 1631685649719702672282505500000 k − k + 28843755938318780823218400000000000 d = − k . polynomial equation of degree 6 be satisfied: f ( x , k ) = (cid:80) i =0 d i x i = 0 (7)where the coefficients d i are univariate polynomials in k of maximum degree 2 as given in Table 2.Finally, the constraints on k are that it be positive; itnot be a root of the polynomial in Equation (6); nor twoother polynomials as described in Table 4.Thus this solution component is valid for all positive val-ues of k excluding three points. As before, we could givethese as exact algebraic numbers but for brevity give floatapproximations: 409 . . . k . However, it does not identify where the number of realsolutions change: although the break point identified earlierhas been rediscovered there is no information from which wecan infer its significance; and there is no significance in ourapplication of the other two isolated points.To finish the analysis we need to decompose ( x , k )-space according to the real roots of f ( x , k ); and also x since the constraint x > f ( x , k ) divides the plane into 135cells in a few seconds. This CAD decomposes the k axisinto 11 cells, i.e. identifying five points which approximateto: − . − . . . k ∈ ]0 , . x , k ) plane is divided into 11 cells: three of which cover x > f ( x , k ) has a single positive real solution for such k . On the two cells for k ∈ ]409 . , . k ∈ ]25084 . , ∞ [ the cylinders above are divided into15 cells; seven of which cover x >
0. This indicates that f ( x , k ) has three positive real solutions for such k .At the end of this analysis we have rediscovered the breakpoint where the system moves from a single positive realsolution to three. We also have explicit solutions valid forall except three isolated k values. To obtain a solutionselect the k value of interest then identify the real rootsof f ( x , k ) (we know in advance how many depending onthe k value chosen); then for each x solution substituterecursively into the equations of Table 3; starting from thebottom and including the new variable solution discoveredfrom each substitution into the next. The solutions in Table1 may be easily rediscovered this way. Repeating the Process for Different Choices of theLone Free Parameter/Fixed Parameter Values.
We may repeat the approach described above for differ-ent choices of free parameter and different choices of fixedparameter values. For example: • With k set to 95 instead of 100 we find that thebreak point between 1 and 3 real positive solutionsmoves to k = 369 . k set to 105 it movesto k = 450 . • Allowing k to be free and fixing k = 200 we findthat there is only ever one positive real solution. • Allowing k to be free and fixing k = 500 we findthe number of positive real solutions moving from 1 to3 to 1 breaking at k = 85 .
988 and k = 110 . able 3: Triangular solution formulae valid for all positive k excluding three isolated points x = − x + 1600 (10 k − x − x + 10 x − x − x + 1600 ( − x + 27 x + 27 k − x − x + x + k − x = 1150 x ( x + x − x − k + x + 150) x = 118200 (69 x + 182 x )( x + x − x − k + x + 150) x = 15364 ( x + x − x − k + x + 150) x x = 50 − x x − x x = 2101 x x x = x + x − x − k + x + 150 x = 2525000101 x + 1000 x + 50500 x = − x − ( − k + 1101 x + 65650) x − (1000 x + ( − k + 200500) x − k + 5050000) x + 150000 x )101 x + (1000 x + 50500) x x = nd where n = 30625833064790009548991419920 x + ( − k + 37749979225487731805273686504663200) x + (14871210647782462053693235920 k − k + 6815925407229297763234036009365120000) x + (1538325448222983229930530049200 k − k + 279241219028720368578809336249748000000) x + (29370341694954648101085099000000 k − k + 3705960282117523242886769213700000000000) x − k d = 232763663752113237974029404420089 x + ( − k + 88646303215205075376308147029677220) x + (113024761399450186949390623074789 k − k + 11682465068391769796632986929072776500) x + (11455232309649034305597048791479020 k − k + 619147207587597001268026254404647600000) x + (290245997063001550130198026458525000 k − k + 14547288529581382252587071541494600000000) x − k − x has at most 6 solutions for a given value of k , according to Equation (7). Table 4: Constraints on k for solution formulae in Table 3 and Equation (7) to be valid k > (cid:54) = 023197989433419579994929 k − k − (cid:54) = 0505465566622475867655547880786544637953790406059982726185509 k − k + 1175510330915205241831243213231417517003037315562884193657451445400 k − k − (cid:54) = 0 Similarly, allowing k to be free and fixing k = 200we find there is only ever one positive real solution;but fixing k = 500 instead we find 3 real solutionsbetween k = 51 .
382 and 58 .
329 and 1 otherwise.The results above hint that there is a shape approximat-ing a narrow paraboloid in ( k , k , k )-space within whichbistability may occur; with bistability available for any k and k value but bounded from below in the k coordinate.We note that these additional experiments all produce re-sults which, as with the one described in detail, are invalidat a handful of isolated values of the free parameter. We use the three linear conservation constraint equations(3) to eliminate x , x , and x from system (1) and sym-bolically compute the Jacobian ˜ J of the obtained reducedsystem. We then numerically compute the eigenvalues of ˜ J for the instances arising from the substitution of the differentpositive fixed points for the variables and the correspondingparameter values.We have used the float approximations for the unique so-lution x (200) with k = 200 and the three solutions x (500)1 ,. . . , x (500)3 for k = 500 in Table 1. For the single positivefixed point x (200) the Jacobian ˜ J ( x (200) ) has eigenvalues withnegative real part only and hence can be shown to be stable.For k = 500 one of the three positive fixed points x (500)2 can be shown to be unstable, as ˜ J ( x (500)2 ) has one eigen-value with positive real part; the other seven had negativereal parts. In contrast x (500)1 and x (500)3 can be shown to bestable. Hence for k = 500 the system is indeed bistable.A verification of the stability of the fixed points usingexact real algebraic numbers by the well-known Routh–Hurwitz criterion is possible algorithmically [15], but seemsto be out of range of current methods for this example. No-tice that also in other studies on multistationarity of sig-nalling pathways [6, 14] the question of stability has notbeen addressed. Finally, we compare our symbolic results with numericalones obtained using the homotopy solver Bertini [2]. Bertinicomputes complex roots of polynomial systems using meth-ods from numerical algebraic geometry [22].For the parameter values as above and k = 500 we ob-tain six solutions, three of which are positive real solutions.For k = 200 we obtain a single positive solution. In bothcases the relevant solutions coincide with the ones obtainedwith our symbolic analyses up to the used numeric precision.However, for larger values of k Bertini produces incor-rect results due to numerical instability. For instance, wefalsely obtain exactly one positive real solution for k =6000 and no positive real solution for k = 10000.Figure 3 shows a Bertini-based grid sampling of parameterregions, varying k between 200 and 1000 and fixing one of k and k while varying the other among the default values(4). While this suffers from the discrete nature of samplingand potentially unreliable results as discussed, it is never-theless useful for the generation of hypothesis about the na-ture of the parameter regions. Figure 3 seems to identify aregion of bistability (in blue) within the parameter space,as hypothesised at the end of Section 2.1.2. The results ofBertini indicate holes in this region (the green dots within the blue). However, computation at these particular pointsreveals these to be the result of numerical errors: where aninsufficiently high precision causes what is actually a posi-tive real solution to appear to have a negative component.It seems there is scope for fruitful interplay between sym-bolic and numeric methods here; with numerics postulatinghypotheses for the symbolic methods to check and refine.
3. CONCLUSIONS AND FUTURE WORK
We have shown that the determination of multistationar-ity of an 11-dimensional MAPK network can be achieved bycombinations of currently available symbolic computationmethods for mixed equality/inequality systems if, for all butpotentially one parameter, numeric values are known. Theaspiration of a semi-algebraic description of the ranges forall parameters in the conservation laws (3) yielding multi-stationarity will now be pursued, with the present resultsdemonstrating that this aspiration may be within reach.As there are many very relevant systems having dimen-sions between 10 and 20 it seems to be worth the effort toenhance and improve the present algorithmic methods, andin particular their combination, to solve such important ap-plication problems for symbolic computation.
Acknowledgments
For our QE-related computations we used two great free soft-ware tools: GNU Parallel for distributing computations onseveral processors, and yEd for visualization of CAD trees.D. Grigoriev is grateful to the grant RSF 16-11-10075 and toMCCME for wonderful working conditions and inspiring at-mosphere. M. Koˇsta has been supported by the DFG/ANRProject STU 483/2-1 SMArT. H. Errami, A. Weber, andO. Radulescu thank the French-German Procope-DAADprogram for partial support of this research. J. Davenport,M. England and T. Sturm are grateful to EU CSA project712689 SC .
4. REFERENCES [1] D. S. Arnon, G. E. Collins, and S. McCallum.Cylindrical algebraic decomposition I: The basicalgorithm.
SIAM J. Comput. , 13(4):865–877, 1984.[2] D. J. Bates, J. D. Hauenstein, A. J. Sommese, andC. W. Wampler. Bertini: Software for numericalalgebraic geometry. doi:10.7274/R0H41PB5.[3] U. S. Bhalla and R. Iyengar. Emergent properties ofnetworks of biological signaling pathways.
Science ,283(5400):381–387, 1999.[4] C. Chen, J. Davenport, J. May, M. Moreno Maza,B. Xia, and R. Xiao. Triangular decomposition ofsemi-algebraic systems.
J. Symb. Comput. , 49:3–26,2013.[5] C. Chen, M. Moreno Maza, B. Xia, and L. Yang.Computing cylindrical algebraic decomposition viatriangular decomposition. In
Proceedings of the ISSAC2009 , pages 95–102. ACM, 2009.[6] C. Conradi, D. Flockerzi, and J. Raisch.Multistationarity in the activation of a MAPK:parametrizing the relevant region in parameter space.
Math. Biosci. , 211(1):105–31, 2008.[7] C. Conradi and M. Mincheva. Catalytic constantsenable the emergence of bistability in dual igure 3: Grid sampling two-parameter regions using Bertini. We combine k with k (left) and with k (right). Colors indicate the computed numbers of positive real fixed points: blue 3, green 2, yellow 1, red 0.The dotted lines indicate values of the parameters as in Equation (4). phosphorylation. Journal of The Royal SocietyInterface , 11(95), 2014.[8] G. Craciun, A. Dickenstein, A. Shiu, and B. Sturmfels.Toric dynamical systems.
J. Symb. Comput. ,44(11):1551–1565, 2009.[9] A. Dolzmann and T. Sturm. Redlog: Computeralgebra meets computer logic.
ACM SIGSAMBulletin , 31(2):2–9, 1997.[10] M. England, R. Bradford, and J. Davenport.Improving the use of equational constraints incylindrical algebraic decomposition. In
Proceedings ofthe ISSAC 2015 , pages 165–172. ACM, 2015.[11] M. England and J. Davenport. The complexity ofcylindrical algebraic decomposition with respect topolynomial degre. In
Proceedings of the CASC 2016 ,volume 9890 of
LNCS , pages 172–192. Springer, 2016.[12] M. Feinberg. Stability of complex isothermalreactors–I. the deficiency zero and deficiency onetheorems.
Chem. Eng. Sci. , 42(10):2229–2268, 1987.[13] D. Grigoriev and N. N. Vorobjov. Solving systems ofpolynomial inequalities in subexponential time.
J. Symb. Comput. , 5:37–64, 1988.[14] E. Gross, H. A. Harrington, Z. Rosen, andB. Sturmfels. Algebraic systems biology: A case studyfor the Wnt pathway.
Bull. Math. Biol. , 78(1):21–51,2016.[15] H. Hong, R. Liska, and S. Steinberg. Testing stabilityby quantifier elimination.
J. Symb. Comput. ,24(2):161–187, 1997.[16] M. D. Johnston. A note on “MAPK networks andtheir capacity for multistationarity due to toric steadystates”. arXiv:1407.5651, 2014.[17] B. Joshi and A. Shiu. A Survey of Methods forDeciding Whether a Reaction Network is Multistationary.
Math. Model. Nat. Phenom. ,10(5):47–67, 2015.[18] C. Li, M. Donizelli, N. Rodriguez, H. Dharuri,L. Endler, V. Chelliah, L. Li, E. He, A. Henry, M. I.Stefan, J. L. Snoep, M. Hucka, N. Le Nov`ere, andC. Laibe. BioModels database: An enhanced, curatedand annotated resource for published quantitativekinetic models.
BMC Systems Biology , 4:92, 2010.[19] N. I. Markevich, J. B. Hoek, and B. N. Kholodenko.Signaling switches and bistability arising frommultisite phosphorylation in protein kinase cascades.
J. Cell Biol. , 164(3):353–359, 2004.[20] M. P´erez Mill´an and A. G. Turjanski. MAPK’snetworks and their capacity for multistationarity dueto toric steady states.
Math. Biosci. , 262:125–37, 2015.[21] S. Schuster and T. H¨ofer. Determining all extremesemi-positive conservation relations in chemicalreaction systems: a test criterion for conservativity.
J.Chem. Soc. Faraday T. , 87(16):2561–2566, 1991.[22] A. J. Sommese, J. Verschelde, and C. W. Wampler.Introduction to numerical algebraic geometry. In
Solving Polynomial Equations: Foundations,Algorithms, and Applications , pages 301–337.Springer, 2005.[23] D. Wang.
Elimination Methods . Springer, 2000.[24] D. Wang and B. Xia. Stability analysis of biologicalsystems with real solution classification. In
Proceedingsof the ISSAC 2005 , pages 354–361. ACM, 2005.[25] V. Weispfenning. Quantifier elimination for realalgebra—the quadratic case and beyond.
Appl. Algebr.Eng. Comm. , 8(2):85–101, 1997.[26] G. Weng, U. S. Bhalla, and R. Iyengar. Complexity inbiological signaling systems.