Polynomial Superlevel Set Representation of the Multistationarity Region of Chemical Reaction Networks
PPolynomial Superlevel Set Representation of theMultistationarity Region of Chemical Reaction Networks
AmirHosein Sadeghimanesh ∗ March 18, 2020
Abstract
In this work a new representation of the multistationarity region of reaction net-works is introduced using the polynomial superlevel sets. The advantages of using thepolynomial superlevel set representation to the former existing representations such asCAD, the finite and the grid representations are discussed. And finally the algorithmsto compute this new representation are provided. The results are given in a generalmathematical formalism of parametric system of equations and therefore can be usedin other applied areas.
Many questions in the application can be modeled by a parametric polynomial system,and therefore to answer them, one needs to explore the properties of these systems. Onesuch case is when one wants to study the multistationarity behavior of a chemical reactionnetwork. Variables in a chemical reaction network are the concentration of the species,which vary as time goes on. Hence we have a dynamical system which is of polynomialtype when the kinetics is assumed to follow the mass action rules. The equilibriums ofthis dynamical system, therefore, are solutions to a system of the polynomial equations.However, the coefficients of the terms in the polynomials in our system involve someparameters. These parameters are usually the rates under which a reaction occurs and thetotal amounts (can be thought of dependency on the initial concentration of the species).These variables and parameters can only attain nonnegative real values. A network iscalled multistationary if there exists a choice of parameters for which the network hasmore than one equilibrium. There are many algorithms developed for answering theyes/no question of multistationarity [6, 8–11, 14, 19, 21]. The input of these algorithmsis a reaction network and the output is the confirmation or rejection of the possibility ofexhibiting multistationary behavior. Nevertheless, more importantly, it is to determinethe parameters where the network has this behavior. Unfortunately, there is less successin this direction.Reviewing the state of the art in the literature, we see that in some works, one isonly focused on a specific reaction network and do heuristic manual calculations to find asuitable parameter [3, 12]. In some other works the system of equations for finding equi-libriums are solved for many random points from the parameter space. This is one wayto approximate the region where the network is multistationary. Recently in [20, PaperIII] a new approach to get a description of the multistationarity region is proposed. In ∗ Bolyai Institute, University of Szeged, Szeged, Hungary; [email protected]. a r X i v : . [ q - b i o . M N ] M a r his method one does not need to solve the system of equations to count the number ofequilibriums. Instead one computes an integral to get the expected number of equilibriumswhen the parameters are equipped by a random distribution. Therefore by choosing theuniform distribution and computing the average number of equilibriums on subhyperrect-angles of the parameter space, one can approximate the multistationarity region as unionof subhyperrectangles. However, looking at a list of hyperrectangles, one may not getmuch information about the geometry of this region such as connectedness or convexivity.In this work, we propose using the polynomial superlevel sets to approximate the unionof these hyperrectangles as a set that can be described by the help of one polynomial.Polynomial superlevel sets are already employed to approximate semialgebraic sets inthe literature and used in control and robust filtering contexts, see [4, 5]. The polynomialsuperlevel set representation can be a more compact representation of the region instead oflisting many hyperrectangles, which each are described as a Cartesian product of intervals.Not only that, to check if a point belongs to the region, one can easily just evaluate thepolynomial in this point and check if the answer is greater than or equal to 1. There aremany other things which one can do with the polynomial superlevel set description of theregion more comfortable than by the union of hyperrectangles. These will be explored inSection 4.The organization of this paper is as follows. The mathematical framework of reactionnetworks and definition of the multistationarity region is given in Section 2. Section 3contains the notations regarding parametric functions and definitions of the finite and thegrid representations of the multistationarity region. In section 4 we define polynomialsuperlevel sets formally and describe how one can algorithmically find a polynomial su-perlevel set representation of a set using the finite and the grid representations. We useit to find the polynomial superlevel set representation of the multistationarity region ofa reaction network. Moreover, finally, in Section 5, we discuss methods that sometimescan speed up computation of the polynomial superlevel set representation by the help ofbisecting algorithms and whenever possible, algorithms of computing expected number ofsolutions independently of solving the system itself. Notations.
Cardinal of a set A is denoted by A ). Let x ∈ Z and n ∈ Z − { } . In this paperwe define x modulo n to be n instead of 0 whenever x is a multiple of n . For a function f : A → A and a point u ∈ A , the level set of f is denoted by L u ( f ) and is definedas { x ∈ A | f ( x ) = u } . For two points a = ( a , · · · , a n ) and b = ( b , · · · , b n ) in R n ,the notation [ a, b ] is used to show the hyperractangle (cid:81) ni =1 [ a i , b i ]. For a subset S of ahyperrectangle B ⊆ R n , let Vol( S ) denote the normalized volume of S with respect to B ,i.e. Vol( S ) = Vol( S )Vol( B ) . When a random vector X = ( X , · · · , X n ) is distributed by a uniform distribution on aset S ⊆ R n , we write X ∼ U ( S ). If X is distributed by a normal distribution with mean µ ∈ R n and variance σ ∈ R > , then we write X ∼ N ( µ, σ ) and we mean X , · · · , X n areidentically and independently distributed by N ( µ i , σ ). The expectation of g ( X ) when X is distributed by a probability distribution q is denoted by E (cid:0) g ( X ) | X ∼ q (cid:1) . In this section, we introduce the concepts of reaction network theory that are neededthroughout the paper with the help of a simple gene regulatory network. One can think2f a gene as a unit encoding information for the synthesis of a product such as a protein.First, a group of DNA binding proteins called transcription factors bind a region of thegene called promoter. Now an enzyme called RNA polymerase starts reading the geneand produces an RNA until it arrives in the terminator region of the gene. Until here iscalled the transcription step. After transcription got completed, the resulted RNA leavesthe nucleus (in eukaryotes) and reaches ribosomes. In ribosomes, the second step, calledtranslation, gets started. Ribosomes assemble a protein from amino acids using the manualguide written in the RNA. A gene encoding a protein recipe is said to be expressed whenit gets transcribed to an RNA, and the RNA translated to the protein. A gene is notalways expressed in a constant rate. There might be proteins that bind the transcriptionfactors or the promoter region and, as a result, inhibits the RNA polymerase starting thetranscription process. On the other hand, there might be other proteins in which theirbinding to the transcription factors or the promoter region enhances the transcription.Consider a simple example from [16, Figure 2], depicted here in Figure 1a. There are threegenes with proteins A , B, and C as their final products. Denote their concentrations attime t by [ A ]( t ), [ B ]( t ) and [ C ]( t ) respectively. The concentration of these proteins willnot remain constant all the time, and therefore we have an ODE describing the variationof the concentrations as time goes on, see Figure 1b. Each protein is degraded with a first-order kinetic with the reaction rate constant k A,d , k B,d and k C,d correspondingly. Protein A activates the expression of the second gene with a Michaelis-Menten kinetics with themaximum rate k B, max and the Michaelis constant k − B,A . The third gene gets activatedby both proteins A and B together with product of two Michaelis-Menten kinetics, withmaximum rate k C, max and Michaelis constants k − A,C and k − B,C . The first gene gets expressedby the rate k A, max in the absence of protein C, and protein C has an inhibitory effect onthe expression of the first gene captured by the denominator (1 + k C,A [ C ]( t )) in the rateexpression. Gene 1 Gene 2 Gene 3 (a) d [ A ]( t ) dt = k A, max ·
11 + k A,C [ C ]( t ) − k A,d [ A ]( t ) d [ B ]( t ) dt = k B, max · k B,A [ A ]( t )1 + k B,A [ A ]( t ) − k B,d [ B ]( t ) d [ C ]( t ) dt = k C, max · k C,A [ A ]( t )1 + k C,A [ A ]( t ) · k C,B [ B ]( t )1 + k B,C [ B ]( t ) − k C,d [ C ]( t ) (b) Figure 1: A regulatory network of 3 genes. (a) This graph shows the relations betweenexpressions of the genes. An inhibitory relation is shown by (cid:173) and a positive relation isshown by → . (b) The system of ODE equations of the network.A solution to the system obtained by letting d [ X i ]( t ) dt = 0 (here X i ’s are A , B and C ) is called an equilibrium of the ODE system. Since the concentration of the proteinscan only be nonnegative real numbers, the complex or negative real solutions are ignored.Sometimes we may only consider the positive solutions, for example, if a total consumptionof a protein is not possible or of interest. Therefore by steady states we mean positivesolutions to the system of equations d [ X i ]( t ) dt = 0. The equations in this system are called thesteady state equations. Now we are ready to define a reaction network formally. A reactionnetwork or a network for short is an ordered pair, N = ( S , R ) where S and R are twofinite sets called the set of species and the set of reactions. In our example, S = { A, B, C } and R contains six reactions; three gene expressions and three protein degradation. To3ach network, an ODE is attached with a concentration of the species as its variables andthe constants of the reaction rate expressions as its parameters. In our example, we have3 variables and 10 parameters. To fix the notation assume S = { X , · · · , X n } and thereare r constants involved in the reaction rate expressions. Then we use x i instead of [ X i ]( t )and k i for the i -th parameter. Denote by f i,k ( x ) the i -th steady state equation where x = ( x , · · · , x n ) and k = ( k , · · · , k r ).The network in Figure 1a is an open network because of the presence of the degradationreactions. A network can also be fully or partially conserved. Consider the simple singlereaction network depicted in Figure 2a. It’s easy to see that x + 2 x = T , x + 2 x = T and x + 2 x = T for real numbers T , T and T determined by the initial conditionsof the ODE system. Therefore three of the linearly redundant steady state equations canbe replaced by these three linear invariants which are called conservation laws in CRN(Chemical Reaction Network theory). The linear subspcae determined by the conservationlaws is called the stoichiometric compatibility class . One should note that the trajectoriesof the ODE system are confined to stoichiometric compatibility classes. In this case oneonly care about the steady states in one stoichiometric compatibility class. –2 + 2 H + k −−→ O + H O (a) dx dt = − kx x , dx dt = kx x dx dt = − kx x , dx dt = kx x (b) Figure 2: A simple example of a closed network consisting of one reaction. (a) Twomolecules of superoxide and two hydron atoms react to each other and produce onemolecule of dioxygen and a molecule of hydrogen perixide. The reaction rate here fol-lows the mass-action kinetics with the reaction rate constant k . (b) The system of ODEequations of the network. The concentrations of O –2 , H + , O and H O are denoted by x , x , x and x respectively.Now we are ready to define the main concept of interest, multistationarity. Definition 2.1.
Consider a network with n species. Replace redundant steady stateequations by conservation laws if there exists any. Let k stands for the vector of constantsof both the reaction rates and conservation laws and of the size r . A network is calledmultistationary over B ⊆ R r if there exists a k ∈ B such that f k ( x ) = 0 has more thanone solution in R n> . Remark 2.2. i) One may also consider non-linear invariants such as first integrals as defined in [18,Definition 11].ii) Note that we are not concerned with the choice of the kinetics such as mass-action,Michaelis-Menten, Hill function, power-law kinetics and s-systems, or the form ofthe steady state equations such as polynomial or rational functions. Therefore theresults of this paper will remain valid and practical for a general reaction network.To answer the question of whether a network is multistationary or not one can useone of many algorithms available in the literature, see [14] and [21] for a few examples.However, to partition the parameter space to two subsets, one consisting of the choices4f parameters for which f k ( x ) has more than one solution and the other comprising thoseparameter choices for which f k ( x ) has at most one solution is a more laborious task. Definition 2.3.
Consider a reaction network with the setting and notation of Definition2.1. The set { k ∈ B | (cid:0) f − k (0) ∩ R n> (cid:1) ≥ } is called the multistationarity region of thenetwork.The region B in Definition 2.3 is usually a hyperrectangle made by the inequality re-strictions of the form k i, min < k i < k i, max for the parameters. For example the rate ofexpression of a gene can not be any arbitrary positive number or the constant of conser-vation laws may be limited from the above due to the limitation of the materials in thelab. Let f k : R n → R m be a parametric function with B ⊆ R r as its parameter region and u apoint in R m . For each choice of the parameters k (cid:63) ∈ B , the system f k (cid:63) ( x ) = u is a non-parametric system of equations. One can solve this system and look at the cardinal of thesolution set. For different choices of k (cid:63) , this number can be different. Therefore we definea new function Φ uf : B → Z ≥ ∪ {∞} sending k ∈ B to (cid:0) L u ( f ) (cid:1) . Now one can partition B to the union of level sets of the map Φ uf . For a general form of f k ( x ), finding L i (Φ uf ) is ahard question. In the case where f k ( x ) ∈ (cid:0) R ( k )[ x ] (cid:1) m and A and B are semialgebraic setsone can use Cylindrical Algebraic Decomposition (CAD). CAD decomposes B to a finitenumber of connected sets called cells. Each cell has intersection with only one L i (Φ uf ) andtherefore L i (Φ uf ) can be expressed as union of finite number of cells with exact descriptionof their boundaries. The problem with this method is that the number of cells growsdoubly exponential on the total number of variables and parameters of f k ( x ). This makesCAD impractical for studying parametric systems of polynomial equations with more thana few variables and parameters. Another approach adopted by scientists is to solve thesystem f k ( x ) = u for many choices of k ∈ B . Mathematically speaking, B is replaced bya finite set. Then each L i (Φ uf ) is expressed as a subset of this finite set. This approachhereafter is referred as the finite representation approach.Since we are motivated from the application, we should note that in a lab, it is usuallynot possible to design the experiment so that the parameter values are exactly the numbersthat we decide. Therefore when the experiment is designed to have k = k (cid:63) , what hapeensis that k is a point in a neighborhood of k (cid:63) and not necessarily k (cid:63) itself. This can happenfor example because of errors coming from the measurement tools or the noise from theenvironment or any other reason depending on the context. In such cases picking upa point close to the boundaries of L i (Φ uf ) can led to a different result than what theexperimentalist expects to see. A different discretization of B can be done using a gridinstead of a finite subset. For example if B is a hyperrectangle [ a, b ] then a grid on B isachieved by dividing B along each axis to equal parts. Then for each subhyperrectanglesof B in this grid we assign the average of the number of solutions of f k ( x ) = u for severalchoices of k coming from the subhyperrectangle. This approach hereafter is referred asthe grid representation approach. See Figure 4 to compare the three approaches visually. Example 3.1.
Consider the gene regulatory network in [22, Figure 3B], depicted herein Figure 3a with the ODE system in Figure 3b. This network has one conservationlaw x + x = k . Therefore we consider the system of equations obtained by the firstthree steady state equations in the ODE system and the conservation law to study the5ultistationarity of this network. We fix the values of all parameters other than k and k to the values given at [22, Figure 4] which are listed below. k = 2 . , k = 1 , k = 0 . , k = 2 . , k = 1 . , k = 46 . . (1)Using the “RootFinding[Parametric]” package of Maple [13] we get the exact description ofthe multistationarity region of the network in the hyperrectangle made by the constraints0 . < k < .
001 and 0 < k <
2. This algorithm uses CAD. The result is depicted inFigure 4a.A finite representation of the same area is found by solving the system of the equationsfor 1000 points ( k , k ) uniformly sampled from [(0 . , , (0 . , . , , (0 . , k , k ) uniformly sampled fromeach subrectangles. Then the subrectangles are colored with respect to the average numberof solutions. This computation also was done by Matlab and took 362 seconds. See Figure4c. X k −−→ X + PP k −−→ P k −− (cid:42)(cid:41) −− k P PX + P P k −− (cid:42)(cid:41) −− k XP PXP P k −−→ XP P + P (a) dx dt = − k x x + k x ,dx dt = k x − k x − k x + 2 k x + k x ,dx dt = k x − k x − k x x + k x ,dx dt = k x x − k x (b) Figure 3: A bistable autoregulatory motif presented in [22, Figure 3B]. (a) X is a gene, P is a protein that can form a dimer P P and then binding to X . The gene X will getexpressed and produce P in both forms X and XP P . And finally there is a degradationof P . (b) The ODE system of the gene regulatory network in part (a). The variables x , x , x and x are standing for the concentration of the species X , P , P P and
XP P respectively.
We open this section with formally defining a polynomial superlevel set.
Definition 4.1.
Consider an arbitrary function f : R n → R . For a given u ∈ R a superlevelset of f is the set of the form U u ( f ) = { x ∈ R n | f ( x ) ≥ u } . When u = 1 we drop the subindex and only write U ( f ). Naturally, a polynomial superlevelset is a superlevel set of a polynomial.Polynomial sublevel sets are defined similarly as in Definition 4.1 with the only differ-ence of the direction of the inequality. However, in this paper, we only focus on superlevelsets. For d ∈ Z ≥ let P d denote the set of polynomials of degree at most d . A sum of squares(SOS) polynomial of degree 2 d is a polynomial p ∈ P d that there exist p , · · · , p m ∈ P d such that p = (cid:80) mi =1 p i . Denote the set of SOS polynomials of degree at most 2 d with Σ d .6 a) k k (b) k k (c) Figure 4: Three representations of the multistationarity region of the network in 3aafter fixing all parameter values other than k and k to the values in equation (1).(a) CAD gives the exact boundary of L (Φ f ) and L (Φ f ), the first one is colored bywhite and the later one with yellow. (b) Finite representation of the parameter region B = [(0 . , , (0 . , L (Φ f ) and are colored by yellow. The 886 other points belong to L (Φ f ) andare colored by sky-blue. (c) The grid representation of B . Each subrectangle is coloredwith respect to the average number of solutions for 10 random points sampled uniformlyfrom the subrectangle. The color bar of the figure is in the right side. Theorem 4.2 ([5, Theorem 2]) . Let B ⊆ R n be a compact set and K a closed subset of B . For d ∈ N define S d = { p ∈ P d | p ≥ on B, p ≥ on K } . Then there exists a polynomial p d ∈ S d such that (cid:90) B p d ( x ) dx = inf { (cid:90) B p ( x ) dx | p ∈ S d } . Furthermore lim d →∞ Vol( U ( p d ) − K ) = 0 . Given a pair (
B, K ) where B ⊆ R n is a compact set and K ⊆ B a closed set and d ∈ N , the superlevel set U ( p ) with p being the polynomial p d ∈ S d found in Theorem4.2 is called the PSS representation of K ⊆ B of degree d . When K is a semialgebraicset, one can find p d numerically using a minimization problem subject to some positivityconstraints [4, Equation 13]. Let B = [ a B , b B ] and K i = [ a K i , b K i ], i = 1 , · · · , m be somehyperrectangles in R n such that K := ∪ mi =1 K i ⊆ B . Using a similar optimization problemit is possible to find the PSS representation of K ⊆ B . Let d ∈ N , the goal is to find thecoefficients of a polynomial of degree d such that (cid:82) B p ( x ) dx becomes minimum subjectto some conditions. Before presenting the constraints, let us look at the target function.A polynomial p ( x ) of degree d can be written as (cid:80) α ∈ N nd c α x α . Where N nd is the set of α = ( α , · · · , α n ) ∈ Z n ≥ such that (cid:80) ni =1 α i = d . Now the integral can be simplified as inbelow. (cid:90) B p ( x ) dx = (cid:90) B ( (cid:88) α ∈ N nd c α x α ) dx = (cid:88) α ∈ N nd c α (cid:90) B x α dx = (cid:88) α ∈ N nd ( (cid:90) B x α dx ) c α . Since (cid:82) B x α dx are constant real numbers independent of the coefficients of the polynomial,the target function is a linear function on the coefficients of p ( x ) which are the main7ariables of the optimization problem. Now looking at the constraints. First of all p ( x )has to be nonnegative on B . This can be enforced by letting p ( x ) − n (cid:88) j =1 s B,j ( x ) (cid:0) x j − a B,j (cid:1)(cid:0) b B,j − x j (cid:1) ∈ Σ r , s B,j ∈ Σ r − , j = 1 , · · · , n, (2)where r = (cid:98) d (cid:99) the largest integer less than or equal to d . Secondly we need p ( x ) ≥ K or in another word p ( x ) − ≥ K . This holds if and only if p ( x ) − ≥ K i .Therefore for every i = 1 , · · · , m one more constraints of the shape (2) has to be added. p ( x ) − − n (cid:88) j =1 s K i ,j ( x ) (cid:0) x j − a K i ,j (cid:1)(cid:0) b K i ,j − x j (cid:1) ∈ Σ r , s K i ,j ∈ Σ r − , j = 1 , · · · , n. Recall Definition 2.3. The multistationarity region of a network is in fact a superlevelset, U (Φ f ). The goal is to find a PSS representation of the set U (Φ f ). One way to ac-complish this goal is to find a grid representation of B and then using the above mentionedSOS optimization problem. The following example illustrates this idea. Example 4.3. (Continued from Example 3.1) Consider the grid representation of themultistationarity region of Example 3.1 given at Figure 4c. To find the PSS representationof this set, we let B = [(0 , , (0 . , . K to be the union of rectangles K i ’s that theirassociated number is greater than 1.8. From the total 100 subrectangles of B , 12 ofthem satisfy this condition. These subrectangles are colored with blue in Figure 5a. Weuse YALMIP and SeDuMi packages of Matlab [7, 17, 23] to solve the SOS optimizationdiscussed before this example. To report the computation time we add the two timesreported in the output of YALMIP which includes the “yalmiptime” and “solvertime”. Ittakes about 1 second to get the coefficients of the polynomial p of the PSS representationof degree 2. Figure 5a shows the plot of U ( p ). Unfortunately for some reason possiblythe numerical issues, the combination of YALMIP and SeDuMi does not give a betterapproximation of p for the PSS representation of this example for higher degrees. Afterrunning the codes for d = 4 , ,
8, where d is degree of p , Matlab plots the same figure asFigure 5a.Consider another gene regulatory example from [20, Chapter 2]. To avoid lengtheningthe text, we only bring the system of equations needed to study the multistationarity ofthe network in equation (3). k x x − k x = 0 k x x − k x = 0 k x − k x = 0 k x − k x = 0 k x x − k x = 0 k x x − k x = 0 k x x − k x = 0 k x x − k x = 0 x = k x = k x + x + x = k x + x + x = k . (3)Fix all parameters other than k and k to the following values coming from equation(2.10) of [20, Chapter 2]. k = k = k = k = 1 , k = 0 . , k = 0 . , k = k = 0 . , k = k = 10000 ,k = 2 , k = 25 , k = 1 , k = 9 , k = k = k = 1 , k = 4 . (4)We reproduced the grid representation of the multistationarity region of this network byMatlab, see Figure 5b. From the total 100 subrectangles of the grid, on 27 of them the8verage number of steady states is greater than 2. Using YALMIP and SeDuMi it tookbetween 3 and 17 seconds to get the polynomial of the PSS representation of degree 2, 4and 6 represented in Figure 5c, 5d and 5e respectively. Luckily for this example, the PSSapproximation of degree 4 gets better than of degree 2 via the code that we wrote. Butfor of degree 6 only a slightly change happens at the bottom of the plot. k k (a) k k (b) k k (c) k k (d) k k (e) Figure 5: PSS representation of different degrees of the mulistationarity region of two generegulatory networks using the information we got from the grid representation. The unionof blue colored subrectangles is considered as the approximation of the multistationarityregion obtained by the grid representation and chosen as the set K . The yellow colored areais the difference of U ( p ) − K . (a) The PSS approximation of K for the network in Figure3a of degree 2 obtained by the information of Figure 4c. (b) The grid representation ofmultistationarity region of the network with the system of equations given in equation (3)and some parameters being fixed by the values in equation (4). (c-e) SPP representations ofthe multistationarity region of degrees 2, 4 and 6 respectively obtained by the informationof Figure 5b. Remark 4.4.
One may ask why should one find a PSS representation of the multistation-arity region using the grid representation if he already has a representation. Let B ⊆ R r be the parameter region of the form of a hyperrectangle, and K ⊆ B be the multistation-arity region. In the grid representation we have K (cid:39) ∪ mi =1 K i where K i = [ a K i , b K i ] arehyperrectangles. In the PSS representation we have K (cid:39) U ( p ) where p is a polynomial ofdegree d .1- When r ≥
4, plotting K is impossible. In order to save or show the grid representa-tion one needs to use a matrix of the size ( m ) × (2 r ), where each row stands for one9 i and the first r columns have the coordinates of the point a K i and the second r columns correspond to the coordinates of the point b K i . However for the PSS rep-resentation one needs to use a vector of the size (cid:80) di =0 (cid:0) r − ir − (cid:1) , where (cid:0) r − ir − (cid:1) entriesare coefficients of the terms of degree i . The terms are ordered from smaller totaldegree to larger and for the terms of the same total degree we use the lexicographicorder.2- To test if a point k (cid:63) ∈ B belongs to K , using the grid representation one shouldcheck m conditions of the form k (cid:63) ∈ K i which means verifying inequality on eachcoordinate of the point i.e. a K i ,j ≤ k (cid:63)j ≤ b K i ,j . If one of the conditions k (cid:63) ∈ K i ispositive, then there is no need to check the rest, otherwise all should fail to concludethat k (cid:63) (cid:54)∈ K . However, using the SPP representation one needs to check only onecondition of an evaluation form, i.e. p ( k (cid:63) ) ≥ k (cid:63) ∈ B to the boundaries of K , using the gridrepresentation one should find distance of k (cid:63) from each K i and then taking theminimum. However using the SPP representation, one just need to find the distanceof k (cid:63) from the algebraic set defined by p ( k ) = 0, for example by the help of Lagrangemultipliers. See Example 4.5. Example 4.5. (Continued from Example 4.3) To illustrate how to approximate distanceof a parameter point from the boundaries of the multistationarity region using a PSSrepresentation consider the network in Example 4.3. Let p be the polynomial of degree4 in two variables k and k corresponding to U ( p ) in Figure 5d. We will approximatedistance of the point k (cid:63) = (0 . , .
02) from the boundary of the multistationarity region bythe distance of k (cid:63) from the algebraic set p ( k , k ) − k , k ) from the point k (cid:63) subject tothe constraint ( k , k ) ∈ L ( p ). The target function is (cid:112) ( k − . + ( k − . whichgets minimized if and only if ( k − . + ( k − . gets minimized. An elementaryway to solve this minimization problem is to use the Lagrangian multipliers. Define F ( k , k , λ ) = ( k − . + ( k − . + λ (cid:0) p ( k , k ) − (cid:1) . Now find the critical points of F ( k , k , λ ). So we should solve the system of equationsobtained by ∂F∂k = ∂F∂k = ∂F∂λ = 0. It takes 1.85 second to solve this system of equationsby Maple. It has 20 solutions, from which 12 of them are complex. From the remainedreal solutions, only one belongs to the rectangle B = [(0 , , (0 . , . . , . k (cid:63) is 0 . B = [ a B , b B ] be a hyperrectangle and K = { a (1) , · · · , a ( m ) } a finite set. Let d ∈ N , thegoal is to find the coefficients of a polynomial of degree d such that (cid:82) B p ( x ) dx becomesminimum subject to some conditions. We already saw that the target function is linear.The constraint p ≥ K can be enforced by p ( a ( i ) ) ≥ i , which are linearconstraints. The positivity of p on B can be enforced by equation (2) or by adding enoughmore number of random points from B and putting the constraints p ( a ) >
0. The lateridea makes the question to be solvable by any common linear programming tool. Howeverhere we still use the equation (2). The following example illustrates how to find the PSSrepresentation via a finite representation. 10 xample 4.6. (Continued from Example 3.1) Consider the finite representation of themultistationarity region of the network of Example 3.1 given at Figure 4b. To find thePSS representation of this set, we let B = [(0 , , (0 . , . K to be the set of pointsfor which the system f k ( x ) = 0 had more than one positive solution. There are 1000points from which 114 of them are parameter choices where the network has three steadystates. Using YALMIP package of Matlab, it takes between 1 and 2 seconds to get thecoefficients of the polynomial p of the PSS representation of degree 2, 6 and 10. Figure 6shows the plot of U ( p ). In contrast to of SPP via grid, SPP via finite gets improved forthis example using Matlab, YALMIP and SeDuMi as we increase d the degree of p . k k (a) k k (b) k k (c) Figure 6: PSS representation of different degrees of the mulistationarity region of thenetwork of Example 3.1 inside the hyperrectangle B = [(0 , , (0 . , . K . The yellow colored area is the difference of U ( p ) − K . As it can beseen, this difference is getting smaller as the degree increases. (a) degree 2. (b) degree 6.(c) degree 10. Remark 4.7.
The same question as in Remark 4.4 can be asked here. Why should onefind a PSS representation of the multistationarity region using the finite representationif he already has a representation. Let B ⊆ R r be the parameter region of the form ofa hyperrectangle, K ⊆ B be the multistationarity region. In the finite representation wehave { a (1) , · · · , a ( m ) } ⊆ K . In the PSS representation we have K (cid:39) U ( p ) where p is apolynomial of degree d .1- When r ≥
4, plotting K is impossible. In order to save or show the finite represen-tation one needs to use a matrix of the size ( m ) × ( r ), where each row stands forone point a ( i ) and the columns correspond to the coordinates of the points. Howeverfor the PSS representation one needs to use a vector of the size (cid:80) di =0 (cid:0) n − in − (cid:1) asexplained in Remark 4.4 item 1.2- To test if a point k (cid:63) ∈ B belongs to K , using the finite representation is not astraightforward task. However, using the SPP representation one needs to verifyonly one condition of the evaluation form, p ( k (cid:63) ) ≥ k (cid:63) ∈ B to the boundaries of K , using the finiterepresentation, if k (cid:63) (cid:54)∈ K , one should compute the distance of k (cid:63) from each point inthe finite representation of K and then take the minimum. However using the SPPrepresentation, whether k (cid:54)∈ K or not, one just need to find the distance of k (cid:63) fromthe algebraic set defined by p ( k ) = 0. 11 Bisect search
If one needs to solve the system in random points and take average in order to get thegrid representation, then using the finite representation and then getting the PSS repre-sentation from the finite representation is a better and faster idea. But in some cases it ispossible to compute the average number of the solutions without solving the system. Onesuch cases is introduced in [20]. Instead of solving the system for many points, it is enoughto compute one integral called Kac-Rice integral. In this situation if the computation ofthe integral is possible and faster than solving the system for many random points, thenthe grid representation can be preferred to the finite representation. However using a gridstill needs a computation per each subhyperrectangle of the grid. This number can growby the number of parameters. If B ⊆ R n and we divide it along each axis to m equal parts,then the number of subhyperrectangles in the grid becomes m n . In this section we in-troduce different decomposition which usually contains less number of subhyperrectangles(not necessarily of equal volume).Let simplify the question. There is a hyperrectangle B ⊆ R r and a function g : B → Z ≥ ∪ {∞} which in our case is Φ f associated to a parametric function f k ( x ). Let B bethe set containing all subhyperrectangles of B . The goal is to express L i ( g ) or U ( g ) asunion of subhyperrectangles of B . One of the common shapes of multistationary networksare bistable networks with folding type of bifurcation. In our settings these networks haveone steady state for some choices of parameters and three steady states for some otherchoices of the parameters and for a zero measure set of parameters in the boundary of thetwo regions it has two steady states. The networks in Examples 3.1 and 4.3 are examplesof such networks. In such cases Φ f almost always takes one of the two values 1 or 3.Going back to our question, motivated from application, assume Im( g ) = { n , n } where n (cid:8) n . In this case for each K ∈ B one of the followings occurs.i) E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = n . This can happen if and only if for almost every k ∈ K , g ( k ) = n .ii) E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = n This can happen if and only if for almost every k ∈ K , g ( k ) = n .iii) E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = α , n (cid:8) α (cid:8) n .This can happen if and only if K ∩ L n ( g ) and K ∩ L n ( g ) both are nonempty andof non-zero measure.The proof is straighforward by notting that E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = n Vol (cid:0) K ∩ L n ( g ) (cid:1) + n Vol (cid:0) K ∩ L n ( g ) (cid:1) . Therefore one can compute E (cid:0) g ( k ) | k ∼ U ( B ) (cid:1) , then if the answer is n concluding thatalmost the whole B is subset of L n ( g ). If the answer is n , concluding that the whole B is subset of L n ( g ). Otherwise dividing B along only one axis to two equal subhyperrect-angles. Continue in this fashion until each subhyperrectangle is inside L n ( g ) or L n ( g )or a termination condition on the length of the edges of the subhyperrectangles happen.When the termination condition on the edges is obtained, put the subhyperrectangle in L n i ( g ) if E ( g ( k )) on this subhyperrectangle is closer to n i . We refer to this approach asthe two-value bisect search hereafter. Two stable and on nonstable which we do not mention stability of the steady states in this paper. lgorithm 5.1. Input: B = [ a, b ] ⊆ R r , g : B → { n , n } ⊆ Z ≥ , n (cid:8) n , (cid:15) ∈ R > . Output: L n ( g ) (cid:39) ∪ m i =1 K n ,i , L n ( g ) (cid:39) ∪ m i =1 K n ,i , ∀ i, j : K n i ,j ∈ B and the minimumof the length of the edges of K n i ,j > (cid:15) . Procedure:Initializing step: L = {} , L = {} , S = { ( B, } . If S (cid:54) = ∅ , choose the first element of S and call its first element by K andits second element by i and remove ( K, i ) from S . Otherwise terminate thealgorithm. Compute α = E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) . If α = n , add K to L . If α = n , add K to L . Otherwise define β to bethe minimum of the length of the edges of K . If β ≤ (cid:15) then add K to L if α ≤ n + n and add K to L if α > n + n .Otherwise define K and K by dividing K along the i -th axis to two equalsubhyperrectangles. Replace i by ( i + 1) mod n and add ( K , i ) and ( K , i ) to S . Return to step 1.If the length of edges of B are of different scales, then it is better to replace the terminationcondition with the following.min { b K,j − a K,j b B,j − a B,j | ≤ j ≤ r } ≤ (cid:15), where K = [ a K , b K ], B = [ a B , b B ] and 0 < (cid:15) < E (cid:0) Φ f ( k ) | k ∼ q (cid:1) where q is a distribution on K and then find a representationusing two-valued bisect search and afterwards a PSS representation. To illustrate thismethod we use the example 2.1 of [20, Paper III] for which the Kac-Rice integral isalready written there. Example 5.2.
Consider the network in Figure 7a taken from [20, Paper III, Example1.2]. The system of equations for studying multistationarity is given at 7b. Fixing allvalues of parameters other than k and k to the following values, we want to find themultistationarity region of the network in the rectangle [(0 , , (5 , k = 0 . , k = 100 , k = 73 . , k = 50 , k = 100 , k = 5 . Figure 7c shows the exact region computed by CAD. Using the Algorithm 5.1 we get theapproximation of this region represented in Figures 7d-7f. As it can be seen, by decreasingthe (cid:15) of the termination condition, the approximation gets improved. Furthermore usingthe Kac-Rice integral given in [20, Paper III] it takes 2.39, 6.88 and 14.61 seconds for ourcode written in Julia to compute the approximations in Figures 7d, 7e and 7f respectively,13hile solving the system in 1000 points to get the finite representation takes 111.70 sec-onds. Figure 7g shows the PSS approximation of degree 4 achieved from Figure 7f. Toavoid numerical/software issues of finding PSS representation via a union of rectangles ashappened in Example 4.3, one can generate random points from rectangles and find thePSS representation from this finite approximation. Adding the time of using Kac-Riceintegral, two-valued bisect search, generating random point and computing PSS represen-tation, all together for this example is still very less than finding a grid representation bysolving the system in many points. The result is shown in Figures 7h and 7i. Remark 5.3.
The reader should note that having less number of hyperrectangles in therepresentation obtained by the two-valued bisect search does not guarantee a faster speedthan a grid search algorithm. Consider the setting in Algorithm 5.1. Assume length of alledges of B ⊆ R n are the same and equal to 2 m . Let (cid:15) = 1 and assume E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) is not getting enough close to n or n for any K in the process of this Algorithm. Thenthe total number of expectations that are needed to be computed until the terminationof this Algorithm is equal to (cid:80) mni =0 i . On the other hand in a grid search by dividing B along each axis to 2 m equal parts, the number of needed expectations to be computed is2 mn .Now consider a more general case where Im( g ) = { n , · · · , n s } ⊆ Z ≥ . In this casewe can not judge about K ∩ L n i ( g ) just by looking at E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) . For exampleif Im( g ) = { , , } and we receive E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = 3, it is not clear that K isalmost subset of L ( g ) or almost half of it is inside L ( g ) and the other half in L ( g ).So now the goal is to find a way to figure out how to decide when to add K to L i when E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = i and when to not add it to L i and instead bisecting it to twosubhyperrectangles in Algorithm 5.1.Note that E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = n Vol (cid:0) K ∩ L n ( g ) (cid:1) + · · · + n s Vol (cid:0) K ∩ L n s ( g ) (cid:1) . Assume { n α , · · · , n α t } ⊆ { n , · · · , n s } such that Vol (cid:0) K ∩ L i ( g ) (cid:1) (cid:54) = 0 only for i ∈ { n α , · · · , n α t } .In that case for any distribution on K which has the same zero measure sets as theLebesgue’s measure we have (cid:82) K ∩ L i ( g ) q ( x ) dx = 0 ; i (cid:54)∈ { n α , · · · , n α t } , (cid:82) K ∩ L i ( g ) q ( x ) dx = (cid:82) K ∩ Li ( g ) q ( x ) dx Vol (cid:0) K ∩ L i ( g ) (cid:1) Vol (cid:0) K ∩ L i ( g ) (cid:1) ; i ∈ { n α , · · · , n α t } , where (cid:82) K ∩ Lni ( g ) q ( x ) dx Vol (cid:0) K ∩ L ni ( g ) (cid:1) ∈ R > is a constant number for a fixed probability density function q .Denote this constant by β q,i when i ∈ { n α , · · · , n α t } and let β q,i = 0 if i (cid:54)∈ { n α , · · · , n α t } .Returning to our goal assume E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) = n i for some i ∈ { , · · · , s } . Iffor every j (cid:54) = i , Vol (cid:0) K ∩ L n j ( g ) (cid:1) = 0, then for any other distribution q on K we have β q,i = 1 and β q,j = 0 for j (cid:54) = i . Therefore E (cid:0) g ( k ) | k ∼ q (cid:1) = n i . Now again assume that { n α , · · · , n α t } ⊆ { n , · · · , n s } such that Vol (cid:0) K ∩ L i ( g ) (cid:1) (cid:54) = 0 only for i ∈ { n α , · · · , n α t } .This time let t ≥
2. Define the following two sets. T = { ( x , · · · , x t ) ∈ (0 , t | x + · · · + x t = 1 } ,T = { ( x , · · · , x t ) ∈ (0 , t | x + · · · + x t = 1 , n α x + · · · + n α t x t = 0 } . To get the grid representation by 100 equal subrectangles and solving for 10 points in each subrectangle,it is again necessary to solve the system for 1000 points. k −−→ X k −−→ X k −−→ X X + X k −−→ X + X X + X k −−→ X + X X k −−→ X (a) k x x − k x = 0 k x x + k x − k x = 0 − k x x + k x − k x = 0 − k x x − k x x + k x = 0 x + x + x + x − k = 0 x + x − k = 0 (b) (c) k k (d) k k (e) k k (f) T T (g) k k (h) k k (i) Figure 7: Using Kac-Rice integral and two-valued bisect search to get the PSS repre-sentation of the multistationarity region of a reaction network presented in (a). Thisnetwork has two conservation laws, therefore two of the six steady state equations are re-placed by these linear invariants. The system of equations for studying multistationarityis given in part (b). (c) The CAD representation of the multistationarity region. (d-f)Approximations of the multistationarity region computed by Algorithm 5.1. (g) The PSSrepresentation of degree 4 computed by considering union of yellow rectangles in (f) as K.(h-i) PSS representations of degree 2 and 4 computed by considering 140 sample pointsfrom the union of yellow rectangles in (f) as approximation of K.15ote that T is a set of one dimension lower than dimension of T . By varying q one canattain any point in T by (cid:0) β q,α n Vol( K ) , · · · , β q,α nt Vol( K ) (cid:1) and E (cid:0) g ( k ) | k ∼ q (cid:1) = n i ifand only if this point belongs to T . Therefore by probability one for randomly chosendistribution q , we will not get E (cid:0) g ( k ) | k ∼ q (cid:1) = n i . Hence we proved the following lemma. Lemma 5.4.
Let B ⊆ R r be a hyperrectangle and g : B → { n , · · · , n s } ⊆ Z ≥ . Assumethat E (cid:0) g ( k ) | k ∼ U ( B ) (cid:1) = n i for some i ∈ { , · · · , s } . Then with probability one we havethat B is almost subset of L n i ( g ) if and only if E (cid:0) g ( k ) | k ∼ q (cid:1) = n i for a randomly chosendistribution q on B with the same zero measure sets as the Lebesgue’s measure. Noting that in [20, Paper III] it is mentioned that the Kac-Rice integral can also beused to compute the expected number of steady states when the parameters are equippedby normal distributions we get the following algorithm which we call it two-step bisectsearch . Algorithm 5.5.
Input: B = [ a, b ] ⊆ R r , g : B → { n , · · · , n s } ⊆ Z ≥ , n (cid:8) · · · (cid:8) n s , (cid:15) ∈ R > . Output: L n ( g ) (cid:39) ∪ m i =1 K n ,i , · · · , L n s ( g ) (cid:39) ∪ m s i =1 K n s ,i , ∀ i, j : K n i ,j ∈ B and theminimum of the length of the edges of K n i ,j > (cid:15) . Procedure:Initializing step: L = {} , · · · , L s = {} , S = { ( B, } . If S = ∅ , then terminate the algorithm. Otherwise choose the first elementof S and denote its first element by K and its second element by i . Remove( K, i ) from S . Compute α = E (cid:0) g ( k ) | k ∼ U ( K ) (cid:1) . Define β to be the minimum of the length of the edges of K . If β ≤ (cid:15) , thenadd K to L j if α is closer to n j and return to step 1. If α (cid:54)∈ { n , · · · , n s } , replace i with ( i + 1) mod n . Define K and K bydividing K along the i -th axis to two equal subhyperrectangles. Add ( K , i )and ( K , i ) to S and return to step 1. Choose a distribution on K randomly and call it q , for example choosea random point from K and call it k (cid:63) , then let q be the truncated normaldistribution on K with mean being k (cid:63) and let the variance to be a number nottoo small and not too large comparing to the length of the edges of K . Compute γ = E (cid:0) g ( k ) | k ∼ q (cid:1) . If γ = n j , then add K to L j , otherwise replace i with ( i + 1) mod n . Define K and K by dividing K along the i -th axis to two equal subhyperrectangles.Add ( K , i ) and ( K , i ) to S . Computer information.
All computations of the examples of this paper are doneon a computer with the following information. It should be noted that the reportedcomputation times may differ on a computer with different properties than the ones inbelow.
Processor: Intel(R) Core(TM) i7-2670QM CPU @2.20GHz 2.20GHz,Installed memory (RAM): 6.00 GB,System type: 64-bit Operating System, x64-based processor.
Python 3.7.4, Julia 1.1.1, Maple 2019, Matlab R2019b, YALMIP, SeDuMi 1.3.
Acknowledgements.
Thanks to Hamed Baghal Ghafari for reading and discussing allthe earlier drafts of the paper and helpful comments about Matlab and L A TEX.This work is supported by NKFIH KKP 129877.
References [1] Robert J. Adler and Jonathan E. Taylor.
Random Fields and Geometry . Springer-Verlag New York, 1st edition, 2007.[2] Afrim Bojnik.
Statistics of real roots of random polynomials . Master’s thesis, Sa-bancıUniversity, 2019.[3] Russell Bradford, James H. Davenport, Matthew England, Hassan Errami, VladimirGerdt, Dima Grigoriev, Charles Hoyt, Marek Koˇsta, Ovidiu Radulescu, ThomasSturm, and G. Andreas Weber.
Identifying the parametric occurrence of multiplesteady states for some biological networks . J. Symb. Comput., 98(Special Issue onSymbolic and Algebraic Computation: ISSAC 2017):84 – 119, 2020.[4] Fabrizio Dabbene and Didier Henrion.
Set approximation via minimum-volume poly-nomial sublevel sets . In European Control Conference (ECC), page 11p., Zurich,Switzerland, July 2013.[5] Fabrizio Dabbene, Didier Henrion, and Constantino Lagoa.
Simple Approximationsof Semialgebraic Sets and their Applications to Control . Automatica, 78:110 – 118,2017.[6] Pete Donnell, Murad Banaji, Anca Marginean, and Casian Pantea.
CoNtRol: anopen source framework for the analysis of chemical reaction networks . Bioinformatics,30(11):1633–1634, 2014.[7] Johan Efberg.
YALMIP : a toolbox for modeling and optimization in MATLAB .In 2004 IEEE International Conference on Robotics and Automation (IEEE Cat.No.04CH37508), pages 284–289, Sep. 2004.[8] Phillipp Raymond Ellison.
The Advanced Deficiency Algorithm and Its Applicationsto Mechanism Discrimination . PhD thesis, The University of Rochester, EastmanSchool of Music, 1998. AAI9905373.[9] Martin Feinberg.
Chemical reaction network structure and the stability of complexisothermal reactors–I. The deficiency zero and deficiency one theorems . Chem. Eng.Sci., 42(10):2229–2268, 1987.[10] Martin Feinberg.
Chemical reaction network structure and the stability of complexisothermal reactors–II. Multiple steady states for networks of deficiency one . Chem.Eng. Sci., 43(1):1–25, 1988.[11] Martin Feinberg.
Multiple steady states for chemical reaction networks of deficiencyone . Arch. Rational Mech. Anal., 132(4):371–406, 1995.1712] Dietrich Flockerzi, Katharina Holstein, and Carsten Conradi.
N-site PhosphorylationSystems with 2N-1 Steady States . Bull. Math. Biol., 76(8):1892–1916, 2014.[13] J¨urgen Gerhard, D. J. Jeffrey, and Guillaume Moroz.
A Package for Solving Para-metric Polynomial Systems . ACM Commun. Comput. Algebra, 43(3/4):61–72, 2010.[14] Haixia Ji, Phillipp Ellison, Daniel Knight, and Martin Feinberg.
CRNToolbox
Version 2-3 — The Chemical Reaction Toolbox . http://crnt.osu.edu/CRNTWin ,2015.[15] Mark Kac. On the average number of real roots of a random algebraic equation . Bull.Amer. Math. Soc., 49(4):314–320, 1943.[16] Guy Karlebach and Ron Shamir.
Modelling and analysis of gene regulatory networks .Nat. Rev. Mol. Cell Biol., 9(10):770–780, 2008.[17] Johan L¨ofberg.
Pre- and post-processing sum-of-squares programs in practice . IEEET. Automat. Contr., 54(5):1007–1011, 2009.[18] Adam Mahdi, Antoni Ferragut, Claudia Valls, and Carsten Wiuf.
Conservation Lawsin Biochemical Reaction Networks . SIAM J. Appl. Dyn. Syst., 16(4):2213–2232, 2017.[19] Stefan M¨uller, Elisenda Feliu, Georg Regensburger, Carsten Conradi, Anne Shiu,and Alicia Dickenstein.
Sign conditions for injectivity of generalized polynomial mapswith applications to chemical reaction networks and real algebraic geometry . Found.Comput. Math., 16(1):69–97, 2016.[20] AmirHosein Sadeghimanesh.
Algebraic tools in the study of Multistationarity of Chem-ical Reaction Networks . PhD thesis, University of Copenhaegn, 10 2018.[21] AmirHosein Sadeghimanesh and Elisenda Feliu.
The multistationarity structure ofnetworks with intermediates and a binomial core network . Bull. Math. Biol., 2019.[22] Dan Siegal-Gaskins, Erich Grotewold, and Gregory D. Smith.
The capacity for mul-tistability in small gene regulatory networks . BMC Syst. Biol., 96(3), 2009.[23] Jos F. Sturm.