# Routes to Multiple Equilibria for Mass-Action Kinetic Systems

Authors: Antonio A. Alonso; Irene Otero-Muras; Manuel Pájaro
Journal: Complexity
 (2018)

---

## Abstract

In this work we explore two potential mechanisms inducing multiple equilibria for weakly reversible networks with mass-action kinetics. The study is performed on a class of polynomial dynamic systems that, under some mild assumptions, are able to accommodate in their state-space form weakly reversible mass-action kinetic networks. The contribution is twofold. We provide an explicit representation of the set of all positive equilibria attained by the system class in terms of a set of (positive parameter dependent) algebraic relations. With this in hand, we prove that deficiency-one networks can only admit multiple equilibria via folding of the equilibrium manifold, whereas a bifurcation leading to multiple branches is only possible in networks with deficiencies larger than one. Interestingly, some kinetic networks within this latter class are capable of sustaining multiple equilibria for any reaction simplex, as we illustrate with one example.


---

## Body

## 1. Introduction

The dynamics of chemical reaction systems is typically described by means of mass-action kinetics. These models constitute an essential tool not only for chemical reactor design and control [1], but also for the analysis of complex biochemical networks involved in the functioning of cells. Mass-action kinetic systems exhibit particularly interesting features related to bifurcation phenomena such as multiple equilibria, oscillations, or chaos [2]. They belong to the class of nonnegative systems [3] and, as such, many of their properties are of interest to a range of applications that go beyond chemistry or chemical engineering reaching fields as diverse as mechanics, economy, biology, or nanoscale communication [4–7].The central objective of this note is to identify the mechanisms leading to the emergence of multiple equilibria in mass-action kinetic systems. Two main routes are linked to multiple equilibria in dynamic systems [8]. One is the appearance of turning points via folding of the equilibrium manifold, also known in the literature as fold, saddle node, or limit point bifurcations. Alternatively, multiple equilibria might emerge from a branching point, in which two or more equilibrium curves with different tangents intersect.In the context of the Chemical Reaction Network Theory (CRNT), the question of existence and uniqueness (or not) of equilibria has attracted most of the attention, although little has been said about the bifurcation mechanisms underlying multiple equilibria. Existence and uniqueness of equilibria under a so-called complex balance condition have been discussed in [9–11]. Some particularly sound results are stated in the deficiency-zero and deficiency-one theorems [12, 13]. Explicit conditions for multiple equilibria were given in [14–17] whereas graph-oriented methods have been proposed in [18]. In connection with the stability of equilibria, a thermodynamic based stability criterion for biochemical networks under uncertainty has been proposed in [19]. In [20], results from Stoichiometric Network Analysis (SNA) have been exploited to detect bistability and oscillations in metabolic networks. A geometric perspective has been adopted in [21] to extend the stability analysis to reaction-diffusion systems. In [22], a robust stability method is proposed to study changes in the qualitative behavior of biochemical networks. The method makes use of a feedback loop breaking approach [23] to establish relationships between equilibrium states, parameters, and the eigenvalues of the corresponding linearized system. Given a biochemical reaction network structure and stoichiometry, necessary and sufficient algebraic conditions for multiple equilibria are proposed in [24, 25]. In [25] an algorithm has been developed to identify equilibrium multiplicity as a function of network parameters. Recent applications in the context of biochemical networks include [26, 27].In this contribution, we explore the bifurcation phenomena leading to multiple equilibria on a class of polynomial dynamic systems that, under some constraints on the structure of the state-space representation, can accommodate a large class of weakly reversible networks obeying mass-action kinetics. We derive a set of algebraic relations that describe the set of all positive equilibria attained by the system class, as a function of some positive parameters. With this result in hand, we find out that network deficiency, as defined in CRNT [12], plays a central role in the way multiplicity arises. Remarkably, deficiency-one networks can only accept multiple equilibria via folding of the equilibrium manifold, thus precluding branching phenomena. Equilibrium branching might therefore only be possible in networks with deficiencies larger than one. Interestingly, we find some kinetic networks with deficiency higher than one capable of sustaining multiple equilibrium for any reaction simplex.The paper is organized as follows. In the Methods we formally describe the class of polynomial dynamic systems under study and illustrate its connections with weakly reversible chemical reaction networks (CRN). In addition, we elaborate on some results related to the notion of feasible equilibria as developed in [28] to be employed in the Results and Discussion to characterize the set of positive equilibria, on the one hand, and to discuss on potential mechanisms leading to multiple equilibria, on the other hand. We complete the section with some examples representative of multiple equilibria via folding and branching. The paper ends up with some conclusions and perspectives for future work.


## 2. Methods

### 2.1. A Class of Dynamic Polynomial Systems

The class of dynamic systems we consider in this note is of the form:(1)c˙=∑λSλϕλc,where the state vector c∈R≥0m is nonnegative and usually denotes species concentrations. The summation extends over λ=1,…,l, and for each λ we have Sλ∈Zm×Nλ and ϕλ:R>0m↦RNλ of the form:(2)ϕλc=ψrλcaλ+Eλfλc,where Eλ∈RNλ×Nλ are a class of invertible Metzler matrices [29], aλ∈RNλ are nonnegative vectors, and both Eλ,aλ contain system parameters.Assumption 1.
In this contribution, we specialize on pairs(aλ,Eλ) such that the resulting vector -Eλ-1aλ>0 (i.e., strictly positive). (As it is extensively discussed in [28],any weakly reversible chemical reaction network satisfies this assumption.)Each functionψrλ:R>0m↦R>0 in expression (2) is a component of the vector function ψr:R>0m↦R>0l, which together with fλ:R>0m↦R>0Nλ takes the form:(3)ψrc=ψlcexp⁡LTln⁡c1,(4)fλc=exp⁡SλTln⁡c,where ψl(c)=exp⁡(yTlnc), y∈Zm is a (positive) coefficient vector, and L∈Zm×(l-1). We refer to (1) as a dynamic polynomial system, which can be written in compact form:(5)c˙=Sϕc,where S∈Zm×(n-l) (with n=∑λ=1lNλ+l) reads S=S1,…,Sλ,…,Sl and ϕ:R>0m↦Rn-l is of the form ϕT=ϕ1T,…,ϕλT,…,ϕlT with subvectors ϕλ as in (2).Definition 2.
We consider two subspaces associated with matricesL in (3) and S in (4)-(5) with the following dimensions:(6)η=dim⁡Im⁡L∩Im⁡S,(7)μ=dim⁡NullL.Remark 3.
Whenl=1, expression (3) reduces to ψr(c)(=ψ1(c))=exp⁡(yTlnc).Although not restricted to them, fory≥0 and appropriate matrices S and L, system (5) can accommodate weakly reversible CRN with mass-action kinetics. This will be illustrated next on the basis of the example in Figure 1.Figure 1
Example1: a reversible chemical reaction network describing a Michaelis-Menten mechanism. The network involves five chemical species [E,S,ES,ESS,P], where S and P represent substrate and product, respectively, transformed under the action of the enzyme E via intermediates ES and ESS. Each maximal set of connected complexes is a linkage class. Following the representation proposed in [28], with each linkage class we associate a reference complex (in squares) that serves to construct the stoichiometric subspace in the first place and will lead to the dynamic structure (1). Assuming mass-action law, we associate with each arrow a chemical reaction coefficient kij. The zero complex (complex number 3 in the graph) represents the environment which in this case exchanges substrate and product with the system.


### 2.2. Invariance and Equilibrium

All possible equilibrium solutionsc∗ of system (5) necessarily correspond to ϕ(c∗) either being a zero vector or lying in the null space of S. Let s=rank(S); then by the rank-nullity theorem, the dimension of the null space of S is δ=n-l-s. Such number is known in CRNT as the deficiency of the chemical reaction network, so we will refer to the null space of S as the deficiency subspace. Let gq∣q=1,…,δ be a basis of the null space of S, and express each vector gq in terms of l subvectors gλq∈RNλ such that(8)gqT=g1qT⋯gλqT⋯glqT.From the above description, the equality Sgq=0, for every q, can be expressed as(9)∑λSλgλq=0forq=1,…,δ.Note that by using matricesGλ∈RNλ×δ of the form:(10)Gλ=gλ1,…,gλq,…,gλδ,expression (9) can be rewritten as ∑λSλGλ=0. On occasions, we may use the matrix GT=G1T,…,GλT,…,GlT.Assumption 4.
We assume in what follows that system (1) is such that Gl is full rank.Starting from an arbitrary reference concentration vectorc0, we define a set(11)Ωc0=c∈Rm∣c≥0,BTc-c0=0,where the columns of B∈Rm×(m-s) are a basis of ImS⊥. Note that σ(c)=BTc is an invariant of motion for system (1), since σ˙≡BTc˙=0, so that any trajectory that starts in Ω(c0) remains there. (Such vector function σ is known in CRNT as the vector of mass conservation constants [26].) In the context of chemical reaction networks, the set Ω(c0) is known as the stoichiometric compatibility class at c0.


### 2.3. Representing CRN Dynamics

A CRN involvingn chemical species is typically represented in CRNT by a graph indicating how some sets of chemicals (reactants) transform into reaction products. Both reactants and products are known as reaction complexes which relate to each other by reaction steps. Reaction steps and complexes correspond in the graph structure with the nodes and edges, respectively. Each maximal set of connected complexes is, according to standard CRNT notation [12], a linkage class of the network. Two complexes are strongly linked if they can reach each other by a directed path. A linkage class is said to be weakly reversible if all the complexes belonging to the linkage class are strongly linked. Weakly reversible networks are those for which all the linkage classes are weakly reversible. Figure 1 shows as an example a typical reversible network describing a Michaelis-Menten mechanism that takes place on a continuous stirred tank reactor (CSTR) under constant volume and input/output fluxes. The network involves the five chemical species E, S, ES, ESS, and P which will be numbered as 1, 2, 3, 4, and 5, respectively. They are distributed in eight complexes organized into three linkage classes. The reactor is such that it only allows S and P to be exchanged with the environment. We collect in columns the stoichiometry of each complex to obtain the molecularity matrix which for this example, following the ordering of species and complexes in Figure 1, becomes as follows.(12)Y=0000011010010100000100010000100001000010Note that the third column is the zero vector (the complex 3 corresponds to the environment). Let Rij be the rate at which a set of reactants in complex i is transformed into a set of products in complex j, and assume that it obeys the mass-action law so that(13)Rij=kijψicwithψic=exp⁡yiTlnc.Under constant volume and flux and taking into account the fact that only S and P are exchanged with the environment, balance equations for the five chemical species become (14)dEdt=-R68+R86+R87-R78dSdt=1γS0-S+R86-R68-R45+R54dESdt=-R45+R54+R68+R78-R86-R87dESSdt=R45-R54dPdt=1γP0-P+R87-R78where γ represents the residence time of the system. At this point it is worth remarking that the above balance equations can be accommodated into the structure given by (1).The reference complexes indicated in Figure1 correspond with columns 1, 4, and 6 and are identified as yr1, yr2, and yr3, respectively. They are employed to construct the stoichiometric subspace as discussed in [28], thus leading to (1), where the Sλ matrices are of the form:(15)S1=00-1-1000010,S2=0-1-110,S3=0-1-1-1010010.Vectors aλ and matrices Eλ in (2) are as follows:(16)a1=01γ,E1=-1γαγ1γ-α+βγ,a2=2,E2=-1,a3=020,E3=-111-2,where α, β represent the input product and substrate concentrations, respectively, and γ the residence time of the system. Matrices E1, E2, and E3 are C-Metzler (see Lemma A.3 in the appendix) so their inverses are nonpositive:(17)E1-1=-γα+ββγαβγβγβ,E2-1=-1,E3-1=-2111.As in Assumption 1, the network is reversible (thus weakly reversible) and satisfies that -Eλ-1aλ>0, for λ=1,2,3.We obtain expression in (3) by constructing L and y from the reference complexes as follows:(18)L=yr1-yr3yr2-yr3=-1-100010000,andy=11000,where L=(yr1-yrl,…,yrl-1-yrl) and y=yrl. Note that, in this way, any weakly reversible mass-action kinetic network can be accommodated into the general expression (1), equivalently (5).The invariant set defined by (11) corresponds with B=10110T. Since the rank of S in (5) is 4 and N1=2, N2=1, and N3=2 so that n=N1+N2+N3+l=8, the dimension of the null space is δ=8-3-4=1. A basis for the null space is given by the vector G=-10010T, so that g1=-10T, g2=0T, and g3=10T.


### 2.4. The Family of Solutions

In this section, we define some auxiliary functions and study a few properties employed later to prove the main result of the paper. With some abuse of notation, let us define, for eachλ=1,…,l, a vector function fλ:R×Rδ↦RNλ as(19)fλxλ;χ=fλ∗+xλhλχ,where fλ∗=-Eλ-1aλ and hλ(χ)=Eλ-1Gλχ, with χ being any element of the set:(20)U=χ∈Rδ∣χ=1.For any χ∈U, there exists an interval Xλ(χ) where fλ(xλ;χ) remains positive. Note that, since by Assumption 1fλ∗>0, the interval(21)Xλχ=xλ∈R∣fλxλ;χ>0is never empty and contains the zero. We define an interval X(χ)=X1(χ)⊗⋯⊗Xλ(χ)⊗⋯⊗Xl(χ) and consider the vector function f:Rl×U↦Rn-l as(22)fTx;χ=f1Tx1;χ,…,fλTxλ;χ,…,flTxl;χ,with x=(x1,…,xλ,…,xl)T. Note that for every x∈X(χ) we have that f(x;χ)>0 (i.e., strictly positive). For convenience, let us define the following matrices:(23)f0=f1∗⋮fλ∗⋮fl∗andHx=x1·E1-1G1⋮xλ·Eλ-1Gλ⋮xl·El-1Gl,to rewrite the vector function (22) as(24)fx;χ=f0+Hx·χ.


#### 2.4.1. A Class of Strictly Monotonous Functions

Let us consider an open intervalX on the real line and define a vector function f(x):R↦Rn as(25)fx=E-1-a+xg,where E∈Rn×n is C-Metzler (thus invertible, see Appendix A), and a,g∈Rn, with g≠0 and a≥0 with at least one positive component. Moreover, let E-1a<0, such that the values of the scalar x for which the vector function f(x) remains strictly positive lie in the interval X=(L-,L+), which includes the zero. Note that since E-1g≠0, X⊂R. Moreover, at least one of the extremes L-,L+ is lower (respectively, upper) bounded, depending on whether the nonzero elements of E-1g have equal or distinct signs (see proofs of Proposition 4.1 and Theorem 5.1 in [28]).Associated withf(x), we have a scalar function F(x):X↦R of the form:(26)Fx=gTlnfx.As it is shown in [28], this function is monotonous decreasing. Such result on monotonicity is formally stated in the lemma below.Lemma 5 (Theorem 5.1 in [28]).
F ( x ) in (26) is monotonous decreasing with a continuous lower bounded derivative on the domain X. Moreover,(27)limx→L+⁡Fx=-∞andlimx→L-⁡Fx=+∞.Remark 6.
One main consequence of Lemma5 is that F(x) defines a bijective map between the set X (function domain) and its range (codomain) R.Forl>1 and λ=1,…,l, let us consider nonzero vectors gλ∈RNλ, C-Metzler matrices Eλ∈RNλ×Nλ, and nonnegative vectors aλ∈RNλ, such that -Eλ-1aλ>0, and construct the set of functions Fλ(xλ):Xλ↦R of the form (26). In addition, let us consider a domain Xρ=X1⊗⋯⊗Xλ⊗⋯⊗Xl∈Rl and construct a function ρ(x):Xρ↦R as(28)ρx=∑λFλxλ,with x=x1,…,xlT. First note that since -Eλ-1aλ>0 for every λ, 0∈Xρ. From Lemma 5 (in particular Remark 6), we can always find elements x∈Xρ such that ρ(x)=0. From Lemma 5 it also follows that each Fλ(xλ) has a continuous negative derivative. Hence, we can conclude that ρ(x) satisfies the conditions of the implicit function theorem (see, e.g., [30]).In this way, for everyλ=1,…,l, there exists a domain X¯λ=X1⊗⋯⊗Xλ-1⊗Xλ+1⊗⋯⊗Xl∈Rl-1, a vector x¯λ=(x1,…,xλ-1,xλ+1,…,xl)T, and a well-defined Lipschitz continuous function Zλ:X¯λ↦R, such that ρ(Zλ(x¯λ),x¯λ)=0. Moreover, since Fλ(xλ) is monotonous decreasing and its derivative is lower bounded, we have that(29)∂Zλ∂xj=-∂Fj/∂xj∂Fλ/∂xλ<0,for  everyλandj≠λ,which implies that each Zλ(x¯λ) has negative (and lower bounded) partial derivatives. Before presenting the result below, let us define for convenience the set Dρ+=Xρ∩R>0l (respectively, Dρ-=Xρ∩R<0l).Theorem 7.
Letz=z1,…,zl-1T be a strictly positive vector and Zl(0)>0 (respectively, <0). Then, there exists a function X:R>0l-1↦Dρ+ (respectively, Dρ-), well defined and Lipschitz continuous, of the form:(30)Xz=Θz·z1,where Θ:R>0l-1↦Xl∩R+ (respectively, Xl∩R-), that solves the following system.(31)ρX=0X1-z1Xl=0⋯=0Xl-1-zl-1Xl=0Proof.
Since0∈Xρ, there are vectors x∗ in Xρ∩R>0l that solve ρ(x∗)=0. To prove the above assertion, note that xl∗=Zl(0)(∈Xl∩R+) satisfies ρ(xl∗,0)=0. By assumption xl∗>0, because of (29), as each xλ increases above zero, xl decreases below xl∗. Thus, repeating this for every λ, by continuity one finds such positive vectors x∗. A similar argument can be employed for Zl(0)∈Xl∩R-. In fact, such vectors are also solutions of system (31) for z∗=z1∗,…,zl-1∗T, with zλ∗=xλ∗/xl∗.
Using a vector functionH(x,z):Xρ×R>0l-1→Rl to represent the left hand side of (31), we then have that H(x∗,z∗)=0 for some z∗∈R>0l-1 and x∗∈Dρ+. In addition, H(x,z) is Lipschitz continuously differentiable everywhere in Dρ+ and R>0l-1. Moreover,(32)∇xHx,z=∂F1∂x1⋯∂Fλ∂xλ⋯∂Fl∂xl10⋯0-z101⋯0-z2⋯⋯⋯⋯⋯00⋯1-zl-1is nonsingular. Hence the implicit function theorem [30] ensures the existence of a function X(z):R>0l-1→Dρ+ well defined and Lipschitz continuous, such that H(X(z),z)=0. Let Θ(z)≡Xl(z); then from (31) we have for each λ=1,…,l-1 that Xλ(z)=zλΘ(z), which leads to expression (30).Corollary 8.
IfZl(0)=0, then Θ(z) is the identically zero function.Proof.
ForZl(0)=0, we have that ρ(0)=0, so for every z>0, X=0 is the only solution of (31), and therefore Θ(z)=0.Remark 9.
Forl=1, (28) reduces to ρ(x)=F1(x1), and because of Lemma 5, there exists only one x1∈X1 such that F1(x1)=0.


#### 2.4.2. The Feasibility Condition

Provided that the logarithm off(x;χ) in (22) lies in the image of ST, any fλ in (19) is compatible with (4). We then say that the family expressed in (22) is feasible.Lemma 10 (feasibility condition).
f ( x ; χ ) in (22) is feasible for any χ∈U and x∈X(χ) such that (gq)Tln⁡f(x;χ)=0, for every q=1,…,δ.Proof.
Since feasibility implies thatln⁡f(x;χ) is in the image of ST, it must be orthogonal to all elements of the basis gq∣q=1,…,δ, which proves the statement.In preparing to explore the consequences of feasibility, for eachχ∈U we define a function ρχ:X(χ)→R, of the form:(33)ρχx=∑λFλxλ;χ,with Fλ≡gλT(χ)lnfλ(xλ;χ), and gλ(χ)=Gλχ. The following lemma describes possible solutions for ρχ(x)=0.Lemma 11.
For everyχ∈U, there exists a well-defined and Lipschitz continuous function Θχ:R>0l-1→R such that(34)Xχz=Θχz·z1,with Xχ(z)∈X(χ), that solves ρχ(x)=0. Moreover, Θχ(z) is strictly positive (respectively, negative) if and only if lnf(0;χ)∉Im(ST). Otherwise, it is the identically zero function. In addition, f(Xχ(z);χ)=f(X-χ(z);-χ).Proof.
l n  f ( 0 ; χ ) ∉ Im ⁡ ( S T ) implies that ρχ(0)≠0, so by Assumption 4, we can find a Zl(0)>0 (respectively, <0) so that Theorem 7 applies, and expression (34) follows with Θχ(z) strictly positive (respectively, negative).
Iflnf(0;χ)∈Im⁡(ST), we necessarily have that ρχ(0)=0, so by Corollary 8 we conclude that Θχ(z) is identically zero.
Concerning the final part, note that for everyχ∈U we have that ρχ(x)=-ρ-χ(-x), so the left hand side of system (31) in the proof of Theorem 7 becomes H(-x,z), which leads to Θχ(z)=-Θ-χ(z). From (34), then we have that Xχ(z)=-X-χ(z). Substituting the above equivalences in (22) we get that f(Xχ(z);χ)=f(-X-χ(z);χ)=f(X-χ(z);-χ) which concludes the proof.Lemma 12.
For everyz∈R>0l-1, there exists at least one vector χ∗∈U for which f(Xχ∗(z);χ∗) is feasible.Proof.
By PropositionB.2, we have that for every z∈R>0l-1 there exists a vector γ∈Rδ (of the form (B.12)) which makes (B.8) satisfy (B.9). This is equivalent to saying that (B.8) is feasible (see Lemma 10). The proof then follows by noting that (B.8) is of the form f(Xχ∗(z);χ∗), with χ∗=γ-1γ.Corollary 13.
Forl=1 there is at least one vector χ∗∈U for which (24) is feasible.Proof.
Ifl=1, then for every χ∈U there exists one scalar xχ∈X(χ) that solves ρχ(x)=0 (Remark 9). Thus for l=1, Xχ(z) in (34) reduces to a constant that depends on χ. Taking this into account, the proof follows the same argument as in Lemma 12.


## 2.1. A Class of Dynamic Polynomial Systems

The class of dynamic systems we consider in this note is of the form:(1)c˙=∑λSλϕλc,where the state vector c∈R≥0m is nonnegative and usually denotes species concentrations. The summation extends over λ=1,…,l, and for each λ we have Sλ∈Zm×Nλ and ϕλ:R>0m↦RNλ of the form:(2)ϕλc=ψrλcaλ+Eλfλc,where Eλ∈RNλ×Nλ are a class of invertible Metzler matrices [29], aλ∈RNλ are nonnegative vectors, and both Eλ,aλ contain system parameters.Assumption 1.
In this contribution, we specialize on pairs(aλ,Eλ) such that the resulting vector -Eλ-1aλ>0 (i.e., strictly positive). (As it is extensively discussed in [28],any weakly reversible chemical reaction network satisfies this assumption.)Each functionψrλ:R>0m↦R>0 in expression (2) is a component of the vector function ψr:R>0m↦R>0l, which together with fλ:R>0m↦R>0Nλ takes the form:(3)ψrc=ψlcexp⁡LTln⁡c1,(4)fλc=exp⁡SλTln⁡c,where ψl(c)=exp⁡(yTlnc), y∈Zm is a (positive) coefficient vector, and L∈Zm×(l-1). We refer to (1) as a dynamic polynomial system, which can be written in compact form:(5)c˙=Sϕc,where S∈Zm×(n-l) (with n=∑λ=1lNλ+l) reads S=S1,…,Sλ,…,Sl and ϕ:R>0m↦Rn-l is of the form ϕT=ϕ1T,…,ϕλT,…,ϕlT with subvectors ϕλ as in (2).Definition 2.
We consider two subspaces associated with matricesL in (3) and S in (4)-(5) with the following dimensions:(6)η=dim⁡Im⁡L∩Im⁡S,(7)μ=dim⁡NullL.Remark 3.
Whenl=1, expression (3) reduces to ψr(c)(=ψ1(c))=exp⁡(yTlnc).Although not restricted to them, fory≥0 and appropriate matrices S and L, system (5) can accommodate weakly reversible CRN with mass-action kinetics. This will be illustrated next on the basis of the example in Figure 1.Figure 1
Example1: a reversible chemical reaction network describing a Michaelis-Menten mechanism. The network involves five chemical species [E,S,ES,ESS,P], where S and P represent substrate and product, respectively, transformed under the action of the enzyme E via intermediates ES and ESS. Each maximal set of connected complexes is a linkage class. Following the representation proposed in [28], with each linkage class we associate a reference complex (in squares) that serves to construct the stoichiometric subspace in the first place and will lead to the dynamic structure (1). Assuming mass-action law, we associate with each arrow a chemical reaction coefficient kij. The zero complex (complex number 3 in the graph) represents the environment which in this case exchanges substrate and product with the system.


## 2.2. Invariance and Equilibrium

All possible equilibrium solutionsc∗ of system (5) necessarily correspond to ϕ(c∗) either being a zero vector or lying in the null space of S. Let s=rank(S); then by the rank-nullity theorem, the dimension of the null space of S is δ=n-l-s. Such number is known in CRNT as the deficiency of the chemical reaction network, so we will refer to the null space of S as the deficiency subspace. Let gq∣q=1,…,δ be a basis of the null space of S, and express each vector gq in terms of l subvectors gλq∈RNλ such that(8)gqT=g1qT⋯gλqT⋯glqT.From the above description, the equality Sgq=0, for every q, can be expressed as(9)∑λSλgλq=0forq=1,…,δ.Note that by using matricesGλ∈RNλ×δ of the form:(10)Gλ=gλ1,…,gλq,…,gλδ,expression (9) can be rewritten as ∑λSλGλ=0. On occasions, we may use the matrix GT=G1T,…,GλT,…,GlT.Assumption 4.
We assume in what follows that system (1) is such that Gl is full rank.Starting from an arbitrary reference concentration vectorc0, we define a set(11)Ωc0=c∈Rm∣c≥0,BTc-c0=0,where the columns of B∈Rm×(m-s) are a basis of ImS⊥. Note that σ(c)=BTc is an invariant of motion for system (1), since σ˙≡BTc˙=0, so that any trajectory that starts in Ω(c0) remains there. (Such vector function σ is known in CRNT as the vector of mass conservation constants [26].) In the context of chemical reaction networks, the set Ω(c0) is known as the stoichiometric compatibility class at c0.


## 2.3. Representing CRN Dynamics

A CRN involvingn chemical species is typically represented in CRNT by a graph indicating how some sets of chemicals (reactants) transform into reaction products. Both reactants and products are known as reaction complexes which relate to each other by reaction steps. Reaction steps and complexes correspond in the graph structure with the nodes and edges, respectively. Each maximal set of connected complexes is, according to standard CRNT notation [12], a linkage class of the network. Two complexes are strongly linked if they can reach each other by a directed path. A linkage class is said to be weakly reversible if all the complexes belonging to the linkage class are strongly linked. Weakly reversible networks are those for which all the linkage classes are weakly reversible. Figure 1 shows as an example a typical reversible network describing a Michaelis-Menten mechanism that takes place on a continuous stirred tank reactor (CSTR) under constant volume and input/output fluxes. The network involves the five chemical species E, S, ES, ESS, and P which will be numbered as 1, 2, 3, 4, and 5, respectively. They are distributed in eight complexes organized into three linkage classes. The reactor is such that it only allows S and P to be exchanged with the environment. We collect in columns the stoichiometry of each complex to obtain the molecularity matrix which for this example, following the ordering of species and complexes in Figure 1, becomes as follows.(12)Y=0000011010010100000100010000100001000010Note that the third column is the zero vector (the complex 3 corresponds to the environment). Let Rij be the rate at which a set of reactants in complex i is transformed into a set of products in complex j, and assume that it obeys the mass-action law so that(13)Rij=kijψicwithψic=exp⁡yiTlnc.Under constant volume and flux and taking into account the fact that only S and P are exchanged with the environment, balance equations for the five chemical species become (14)dEdt=-R68+R86+R87-R78dSdt=1γS0-S+R86-R68-R45+R54dESdt=-R45+R54+R68+R78-R86-R87dESSdt=R45-R54dPdt=1γP0-P+R87-R78where γ represents the residence time of the system. At this point it is worth remarking that the above balance equations can be accommodated into the structure given by (1).The reference complexes indicated in Figure1 correspond with columns 1, 4, and 6 and are identified as yr1, yr2, and yr3, respectively. They are employed to construct the stoichiometric subspace as discussed in [28], thus leading to (1), where the Sλ matrices are of the form:(15)S1=00-1-1000010,S2=0-1-110,S3=0-1-1-1010010.Vectors aλ and matrices Eλ in (2) are as follows:(16)a1=01γ,E1=-1γαγ1γ-α+βγ,a2=2,E2=-1,a3=020,E3=-111-2,where α, β represent the input product and substrate concentrations, respectively, and γ the residence time of the system. Matrices E1, E2, and E3 are C-Metzler (see Lemma A.3 in the appendix) so their inverses are nonpositive:(17)E1-1=-γα+ββγαβγβγβ,E2-1=-1,E3-1=-2111.As in Assumption 1, the network is reversible (thus weakly reversible) and satisfies that -Eλ-1aλ>0, for λ=1,2,3.We obtain expression in (3) by constructing L and y from the reference complexes as follows:(18)L=yr1-yr3yr2-yr3=-1-100010000,andy=11000,where L=(yr1-yrl,…,yrl-1-yrl) and y=yrl. Note that, in this way, any weakly reversible mass-action kinetic network can be accommodated into the general expression (1), equivalently (5).The invariant set defined by (11) corresponds with B=10110T. Since the rank of S in (5) is 4 and N1=2, N2=1, and N3=2 so that n=N1+N2+N3+l=8, the dimension of the null space is δ=8-3-4=1. A basis for the null space is given by the vector G=-10010T, so that g1=-10T, g2=0T, and g3=10T.


## 2.4. The Family of Solutions

In this section, we define some auxiliary functions and study a few properties employed later to prove the main result of the paper. With some abuse of notation, let us define, for eachλ=1,…,l, a vector function fλ:R×Rδ↦RNλ as(19)fλxλ;χ=fλ∗+xλhλχ,where fλ∗=-Eλ-1aλ and hλ(χ)=Eλ-1Gλχ, with χ being any element of the set:(20)U=χ∈Rδ∣χ=1.For any χ∈U, there exists an interval Xλ(χ) where fλ(xλ;χ) remains positive. Note that, since by Assumption 1fλ∗>0, the interval(21)Xλχ=xλ∈R∣fλxλ;χ>0is never empty and contains the zero. We define an interval X(χ)=X1(χ)⊗⋯⊗Xλ(χ)⊗⋯⊗Xl(χ) and consider the vector function f:Rl×U↦Rn-l as(22)fTx;χ=f1Tx1;χ,…,fλTxλ;χ,…,flTxl;χ,with x=(x1,…,xλ,…,xl)T. Note that for every x∈X(χ) we have that f(x;χ)>0 (i.e., strictly positive). For convenience, let us define the following matrices:(23)f0=f1∗⋮fλ∗⋮fl∗andHx=x1·E1-1G1⋮xλ·Eλ-1Gλ⋮xl·El-1Gl,to rewrite the vector function (22) as(24)fx;χ=f0+Hx·χ.


### 2.4.1. A Class of Strictly Monotonous Functions

Let us consider an open intervalX on the real line and define a vector function f(x):R↦Rn as(25)fx=E-1-a+xg,where E∈Rn×n is C-Metzler (thus invertible, see Appendix A), and a,g∈Rn, with g≠0 and a≥0 with at least one positive component. Moreover, let E-1a<0, such that the values of the scalar x for which the vector function f(x) remains strictly positive lie in the interval X=(L-,L+), which includes the zero. Note that since E-1g≠0, X⊂R. Moreover, at least one of the extremes L-,L+ is lower (respectively, upper) bounded, depending on whether the nonzero elements of E-1g have equal or distinct signs (see proofs of Proposition 4.1 and Theorem 5.1 in [28]).Associated withf(x), we have a scalar function F(x):X↦R of the form:(26)Fx=gTlnfx.As it is shown in [28], this function is monotonous decreasing. Such result on monotonicity is formally stated in the lemma below.Lemma 5 (Theorem 5.1 in [28]).
F ( x ) in (26) is monotonous decreasing with a continuous lower bounded derivative on the domain X. Moreover,(27)limx→L+⁡Fx=-∞andlimx→L-⁡Fx=+∞.Remark 6.
One main consequence of Lemma5 is that F(x) defines a bijective map between the set X (function domain) and its range (codomain) R.Forl>1 and λ=1,…,l, let us consider nonzero vectors gλ∈RNλ, C-Metzler matrices Eλ∈RNλ×Nλ, and nonnegative vectors aλ∈RNλ, such that -Eλ-1aλ>0, and construct the set of functions Fλ(xλ):Xλ↦R of the form (26). In addition, let us consider a domain Xρ=X1⊗⋯⊗Xλ⊗⋯⊗Xl∈Rl and construct a function ρ(x):Xρ↦R as(28)ρx=∑λFλxλ,with x=x1,…,xlT. First note that since -Eλ-1aλ>0 for every λ, 0∈Xρ. From Lemma 5 (in particular Remark 6), we can always find elements x∈Xρ such that ρ(x)=0. From Lemma 5 it also follows that each Fλ(xλ) has a continuous negative derivative. Hence, we can conclude that ρ(x) satisfies the conditions of the implicit function theorem (see, e.g., [30]).In this way, for everyλ=1,…,l, there exists a domain X¯λ=X1⊗⋯⊗Xλ-1⊗Xλ+1⊗⋯⊗Xl∈Rl-1, a vector x¯λ=(x1,…,xλ-1,xλ+1,…,xl)T, and a well-defined Lipschitz continuous function Zλ:X¯λ↦R, such that ρ(Zλ(x¯λ),x¯λ)=0. Moreover, since Fλ(xλ) is monotonous decreasing and its derivative is lower bounded, we have that(29)∂Zλ∂xj=-∂Fj/∂xj∂Fλ/∂xλ<0,for  everyλandj≠λ,which implies that each Zλ(x¯λ) has negative (and lower bounded) partial derivatives. Before presenting the result below, let us define for convenience the set Dρ+=Xρ∩R>0l (respectively, Dρ-=Xρ∩R<0l).Theorem 7.
Letz=z1,…,zl-1T be a strictly positive vector and Zl(0)>0 (respectively, <0). Then, there exists a function X:R>0l-1↦Dρ+ (respectively, Dρ-), well defined and Lipschitz continuous, of the form:(30)Xz=Θz·z1,where Θ:R>0l-1↦Xl∩R+ (respectively, Xl∩R-), that solves the following system.(31)ρX=0X1-z1Xl=0⋯=0Xl-1-zl-1Xl=0Proof.
Since0∈Xρ, there are vectors x∗ in Xρ∩R>0l that solve ρ(x∗)=0. To prove the above assertion, note that xl∗=Zl(0)(∈Xl∩R+) satisfies ρ(xl∗,0)=0. By assumption xl∗>0, because of (29), as each xλ increases above zero, xl decreases below xl∗. Thus, repeating this for every λ, by continuity one finds such positive vectors x∗. A similar argument can be employed for Zl(0)∈Xl∩R-. In fact, such vectors are also solutions of system (31) for z∗=z1∗,…,zl-1∗T, with zλ∗=xλ∗/xl∗.
Using a vector functionH(x,z):Xρ×R>0l-1→Rl to represent the left hand side of (31), we then have that H(x∗,z∗)=0 for some z∗∈R>0l-1 and x∗∈Dρ+. In addition, H(x,z) is Lipschitz continuously differentiable everywhere in Dρ+ and R>0l-1. Moreover,(32)∇xHx,z=∂F1∂x1⋯∂Fλ∂xλ⋯∂Fl∂xl10⋯0-z101⋯0-z2⋯⋯⋯⋯⋯00⋯1-zl-1is nonsingular. Hence the implicit function theorem [30] ensures the existence of a function X(z):R>0l-1→Dρ+ well defined and Lipschitz continuous, such that H(X(z),z)=0. Let Θ(z)≡Xl(z); then from (31) we have for each λ=1,…,l-1 that Xλ(z)=zλΘ(z), which leads to expression (30).Corollary 8.
IfZl(0)=0, then Θ(z) is the identically zero function.Proof.
ForZl(0)=0, we have that ρ(0)=0, so for every z>0, X=0 is the only solution of (31), and therefore Θ(z)=0.Remark 9.
Forl=1, (28) reduces to ρ(x)=F1(x1), and because of Lemma 5, there exists only one x1∈X1 such that F1(x1)=0.


### 2.4.2. The Feasibility Condition

Provided that the logarithm off(x;χ) in (22) lies in the image of ST, any fλ in (19) is compatible with (4). We then say that the family expressed in (22) is feasible.Lemma 10 (feasibility condition).
f ( x ; χ ) in (22) is feasible for any χ∈U and x∈X(χ) such that (gq)Tln⁡f(x;χ)=0, for every q=1,…,δ.Proof.
Since feasibility implies thatln⁡f(x;χ) is in the image of ST, it must be orthogonal to all elements of the basis gq∣q=1,…,δ, which proves the statement.In preparing to explore the consequences of feasibility, for eachχ∈U we define a function ρχ:X(χ)→R, of the form:(33)ρχx=∑λFλxλ;χ,with Fλ≡gλT(χ)lnfλ(xλ;χ), and gλ(χ)=Gλχ. The following lemma describes possible solutions for ρχ(x)=0.Lemma 11.
For everyχ∈U, there exists a well-defined and Lipschitz continuous function Θχ:R>0l-1→R such that(34)Xχz=Θχz·z1,with Xχ(z)∈X(χ), that solves ρχ(x)=0. Moreover, Θχ(z) is strictly positive (respectively, negative) if and only if lnf(0;χ)∉Im(ST). Otherwise, it is the identically zero function. In addition, f(Xχ(z);χ)=f(X-χ(z);-χ).Proof.
l n  f ( 0 ; χ ) ∉ Im ⁡ ( S T ) implies that ρχ(0)≠0, so by Assumption 4, we can find a Zl(0)>0 (respectively, <0) so that Theorem 7 applies, and expression (34) follows with Θχ(z) strictly positive (respectively, negative).
Iflnf(0;χ)∈Im⁡(ST), we necessarily have that ρχ(0)=0, so by Corollary 8 we conclude that Θχ(z) is identically zero.
Concerning the final part, note that for everyχ∈U we have that ρχ(x)=-ρ-χ(-x), so the left hand side of system (31) in the proof of Theorem 7 becomes H(-x,z), which leads to Θχ(z)=-Θ-χ(z). From (34), then we have that Xχ(z)=-X-χ(z). Substituting the above equivalences in (22) we get that f(Xχ(z);χ)=f(-X-χ(z);χ)=f(X-χ(z);-χ) which concludes the proof.Lemma 12.
For everyz∈R>0l-1, there exists at least one vector χ∗∈U for which f(Xχ∗(z);χ∗) is feasible.Proof.
By PropositionB.2, we have that for every z∈R>0l-1 there exists a vector γ∈Rδ (of the form (B.12)) which makes (B.8) satisfy (B.9). This is equivalent to saying that (B.8) is feasible (see Lemma 10). The proof then follows by noting that (B.8) is of the form f(Xχ∗(z);χ∗), with χ∗=γ-1γ.Corollary 13.
Forl=1 there is at least one vector χ∗∈U for which (24) is feasible.Proof.
Ifl=1, then for every χ∈U there exists one scalar xχ∈X(χ) that solves ρχ(x)=0 (Remark 9). Thus for l=1, Xχ(z) in (34) reduces to a constant that depends on χ. Taking this into account, the proof follows the same argument as in Lemma 12.


## 2.4.1. A Class of Strictly Monotonous Functions

Let us consider an open intervalX on the real line and define a vector function f(x):R↦Rn as(25)fx=E-1-a+xg,where E∈Rn×n is C-Metzler (thus invertible, see Appendix A), and a,g∈Rn, with g≠0 and a≥0 with at least one positive component. Moreover, let E-1a<0, such that the values of the scalar x for which the vector function f(x) remains strictly positive lie in the interval X=(L-,L+), which includes the zero. Note that since E-1g≠0, X⊂R. Moreover, at least one of the extremes L-,L+ is lower (respectively, upper) bounded, depending on whether the nonzero elements of E-1g have equal or distinct signs (see proofs of Proposition 4.1 and Theorem 5.1 in [28]).Associated withf(x), we have a scalar function F(x):X↦R of the form:(26)Fx=gTlnfx.As it is shown in [28], this function is monotonous decreasing. Such result on monotonicity is formally stated in the lemma below.Lemma 5 (Theorem 5.1 in [28]).
F ( x ) in (26) is monotonous decreasing with a continuous lower bounded derivative on the domain X. Moreover,(27)limx→L+⁡Fx=-∞andlimx→L-⁡Fx=+∞.Remark 6.
One main consequence of Lemma5 is that F(x) defines a bijective map between the set X (function domain) and its range (codomain) R.Forl>1 and λ=1,…,l, let us consider nonzero vectors gλ∈RNλ, C-Metzler matrices Eλ∈RNλ×Nλ, and nonnegative vectors aλ∈RNλ, such that -Eλ-1aλ>0, and construct the set of functions Fλ(xλ):Xλ↦R of the form (26). In addition, let us consider a domain Xρ=X1⊗⋯⊗Xλ⊗⋯⊗Xl∈Rl and construct a function ρ(x):Xρ↦R as(28)ρx=∑λFλxλ,with x=x1,…,xlT. First note that since -Eλ-1aλ>0 for every λ, 0∈Xρ. From Lemma 5 (in particular Remark 6), we can always find elements x∈Xρ such that ρ(x)=0. From Lemma 5 it also follows that each Fλ(xλ) has a continuous negative derivative. Hence, we can conclude that ρ(x) satisfies the conditions of the implicit function theorem (see, e.g., [30]).In this way, for everyλ=1,…,l, there exists a domain X¯λ=X1⊗⋯⊗Xλ-1⊗Xλ+1⊗⋯⊗Xl∈Rl-1, a vector x¯λ=(x1,…,xλ-1,xλ+1,…,xl)T, and a well-defined Lipschitz continuous function Zλ:X¯λ↦R, such that ρ(Zλ(x¯λ),x¯λ)=0. Moreover, since Fλ(xλ) is monotonous decreasing and its derivative is lower bounded, we have that(29)∂Zλ∂xj=-∂Fj/∂xj∂Fλ/∂xλ<0,for  everyλandj≠λ,which implies that each Zλ(x¯λ) has negative (and lower bounded) partial derivatives. Before presenting the result below, let us define for convenience the set Dρ+=Xρ∩R>0l (respectively, Dρ-=Xρ∩R<0l).Theorem 7.
Letz=z1,…,zl-1T be a strictly positive vector and Zl(0)>0 (respectively, <0). Then, there exists a function X:R>0l-1↦Dρ+ (respectively, Dρ-), well defined and Lipschitz continuous, of the form:(30)Xz=Θz·z1,where Θ:R>0l-1↦Xl∩R+ (respectively, Xl∩R-), that solves the following system.(31)ρX=0X1-z1Xl=0⋯=0Xl-1-zl-1Xl=0Proof.
Since0∈Xρ, there are vectors x∗ in Xρ∩R>0l that solve ρ(x∗)=0. To prove the above assertion, note that xl∗=Zl(0)(∈Xl∩R+) satisfies ρ(xl∗,0)=0. By assumption xl∗>0, because of (29), as each xλ increases above zero, xl decreases below xl∗. Thus, repeating this for every λ, by continuity one finds such positive vectors x∗. A similar argument can be employed for Zl(0)∈Xl∩R-. In fact, such vectors are also solutions of system (31) for z∗=z1∗,…,zl-1∗T, with zλ∗=xλ∗/xl∗.
Using a vector functionH(x,z):Xρ×R>0l-1→Rl to represent the left hand side of (31), we then have that H(x∗,z∗)=0 for some z∗∈R>0l-1 and x∗∈Dρ+. In addition, H(x,z) is Lipschitz continuously differentiable everywhere in Dρ+ and R>0l-1. Moreover,(32)∇xHx,z=∂F1∂x1⋯∂Fλ∂xλ⋯∂Fl∂xl10⋯0-z101⋯0-z2⋯⋯⋯⋯⋯00⋯1-zl-1is nonsingular. Hence the implicit function theorem [30] ensures the existence of a function X(z):R>0l-1→Dρ+ well defined and Lipschitz continuous, such that H(X(z),z)=0. Let Θ(z)≡Xl(z); then from (31) we have for each λ=1,…,l-1 that Xλ(z)=zλΘ(z), which leads to expression (30).Corollary 8.
IfZl(0)=0, then Θ(z) is the identically zero function.Proof.
ForZl(0)=0, we have that ρ(0)=0, so for every z>0, X=0 is the only solution of (31), and therefore Θ(z)=0.Remark 9.
Forl=1, (28) reduces to ρ(x)=F1(x1), and because of Lemma 5, there exists only one x1∈X1 such that F1(x1)=0.


## 2.4.2. The Feasibility Condition

Provided that the logarithm off(x;χ) in (22) lies in the image of ST, any fλ in (19) is compatible with (4). We then say that the family expressed in (22) is feasible.Lemma 10 (feasibility condition).
f ( x ; χ ) in (22) is feasible for any χ∈U and x∈X(χ) such that (gq)Tln⁡f(x;χ)=0, for every q=1,…,δ.Proof.
Since feasibility implies thatln⁡f(x;χ) is in the image of ST, it must be orthogonal to all elements of the basis gq∣q=1,…,δ, which proves the statement.In preparing to explore the consequences of feasibility, for eachχ∈U we define a function ρχ:X(χ)→R, of the form:(33)ρχx=∑λFλxλ;χ,with Fλ≡gλT(χ)lnfλ(xλ;χ), and gλ(χ)=Gλχ. The following lemma describes possible solutions for ρχ(x)=0.Lemma 11.
For everyχ∈U, there exists a well-defined and Lipschitz continuous function Θχ:R>0l-1→R such that(34)Xχz=Θχz·z1,with Xχ(z)∈X(χ), that solves ρχ(x)=0. Moreover, Θχ(z) is strictly positive (respectively, negative) if and only if lnf(0;χ)∉Im(ST). Otherwise, it is the identically zero function. In addition, f(Xχ(z);χ)=f(X-χ(z);-χ).Proof.
l n  f ( 0 ; χ ) ∉ Im ⁡ ( S T ) implies that ρχ(0)≠0, so by Assumption 4, we can find a Zl(0)>0 (respectively, <0) so that Theorem 7 applies, and expression (34) follows with Θχ(z) strictly positive (respectively, negative).
Iflnf(0;χ)∈Im⁡(ST), we necessarily have that ρχ(0)=0, so by Corollary 8 we conclude that Θχ(z) is identically zero.
Concerning the final part, note that for everyχ∈U we have that ρχ(x)=-ρ-χ(-x), so the left hand side of system (31) in the proof of Theorem 7 becomes H(-x,z), which leads to Θχ(z)=-Θ-χ(z). From (34), then we have that Xχ(z)=-X-χ(z). Substituting the above equivalences in (22) we get that f(Xχ(z);χ)=f(-X-χ(z);χ)=f(X-χ(z);-χ) which concludes the proof.Lemma 12.
For everyz∈R>0l-1, there exists at least one vector χ∗∈U for which f(Xχ∗(z);χ∗) is feasible.Proof.
By PropositionB.2, we have that for every z∈R>0l-1 there exists a vector γ∈Rδ (of the form (B.12)) which makes (B.8) satisfy (B.9). This is equivalent to saying that (B.8) is feasible (see Lemma 10). The proof then follows by noting that (B.8) is of the form f(Xχ∗(z);χ∗), with χ∗=γ-1γ.Corollary 13.
Forl=1 there is at least one vector χ∗∈U for which (24) is feasible.Proof.
Ifl=1, then for every χ∈U there exists one scalar xχ∈X(χ) that solves ρχ(x)=0 (Remark 9). Thus for l=1, Xχ(z) in (34) reduces to a constant that depends on χ. Taking this into account, the proof follows the same argument as in Lemma 12.


## 3. Results and Discussion

Here we make use of the feasibility condition discussed in the previous section to derive a general representation of the set of all (positive) equilibrium solutions associated with system (1). This will be employed afterwards to identify potential classes of bifurcations.


### 3.1. Characterizing the Equilibrium Set

Depending on the dimensions of subspacesIm⁡(L)∩Im⁡(S) and Null(L) (Definition 2), two scenarios will be considered. In the first place, we assume that both are empty and get the following result.Theorem 14.
Let the dimensionsη=μ=0 (Definition 2), and lnf(0;χ)∉Im⁡(ST). Then, for each z∈R>0l-1 there is at least one vector χ∗∈U, a real number τ∈R with sgn⁡(τ)=sgn⁡(Θχ∗(z)), and vectors ξ∈Rm, such that(35)STξ=lnfz;χ∗(36)LTξ=-lnz(37)yTξ=lnτΘχ∗z.Moreover, any c=exp(ξ) satisfying (35)-(37) is an equilibrium solution of the kinetic system (1).Proof.
In order to prove thatξ simultaneously solves (35)-(37), let us define the matrix M=M0∣y, with M0=S∣L, and rewrite these equations as(38)MTξ=Qz,τ;χ∗,withQz,τ;χ∗=lnfz;χ∗-lnzlnτΘχ∗z.Suppose firstly that y∉Im(L)⊕Im(S). Then, a basis for the kernel of M consists of a vector set wq∣q=1,…,δ, where each element is of the form wqT=gqT0lT, with 0l being the l dimensional vector zero. Note that each element is in the null space by construction since the corresponding gq is in the null space of S, and the δ elements are linearly independent.
In addition, because of dimensionsη=μ=0, no other vector linearly independent to the above defined vector set can be found satisfying that Mw=0. If such a vector w∈Rn, with wT=wSTwLTwyT, and wS,wL,wy≠0 exists, then SwS+LwL+ywy=0. But such a relation would require the overlapping of the corresponding images of S, L, and y, which is against the assumption that the dimensions of subspaces Im(L)∩Im(S) and Null(L) are zero (η=μ=0).
Thatξ solves (38) is equivalent to saying that Q(z,τ;χ∗)∈Im(MT) for each vector z∈R>0l-1 and τ∈R, with sgn(τ)=sgn(Θχ∗(z)), which in turn requires that wqTQ(z,τ;χ∗)=0, for every q=1,…,δ.
Becauselnf(0;χ)∉Im(ST), by Lemma 11 we have that, for every χ∈U, Θχ(z) in (34) has a positive (respectively, negative) sign. From Lemma 12, it follows that for each vector z∈R>0l-1 there is at least one vector χ∗∈U such that f(Xχ∗(z);χ∗) is feasible (with Xχ∗(z) being a vector function of the form (34)) and we get expression (35).
Thatξ solves (38) for each z∈R>0l-1 and any τ∈R is a consequence of the structure of each vector of the basis, so that(39)wqTQz,τ;χ∗=gqTlnfz;χ∗=0.Suppose, on the other hand, that y∈Im(L)⊕Im(S). Then y can be expressed as a linear combination of the columns of L and S and the dimension of the kernel of M becomes δ+1. The new element of the basis can be expressed as a vector wδ+1T=wSTwLT1. In order for Q(z,τ;χ∗) to satisfy (38) we have that wδ+1TQ(z,τ;χ∗)=0. This equation can be solved in τ resulting in a function:(40)τz;χ∗=Θχ∗zΨz,χ∗, with Ψ(z,χ∗)=exp⁡wSTlnf(z;χ∗)-wLTln⁡z. Consequently, for each z and χ∗ there exists one τ and vectors ξ∈Rm satisfying (38).
In order to prove that anyc=exp(ξ) satisfying (35)-(37) is in fact an equilibrium solution, we combine the above expressions with (3)-(4) and make use of (19) to get(41)ψrc≡τΘχ∗zz-11,(42)fλc≡Eλ-1-aλ+zλΘχ∗zGλχ∗,λ≠l(43)flc≡El-1-al+Θχ∗zGlχ∗,where z-1 in (42) should be understood as the element-wise inverse of vector z. Combining (41)-(43) with (2), we also obtain that ϕλ(c)≡τGλχ∗ for λ=1,…,l. Finally, substituting previous expressions in (1), we get(44)c˙=τ∑λSλGλχ∗=0,which proves that any c=exp(ξ) satisfying (35)-(37) is in fact an equilibrium solution of the kinetic system (1).Ifη or μ in Definition 2 is positive, then Theorem 14 does not hold for all z∈R>0l-1, but for a subset of it. In preparing for this result, let us define the matrix:(45)R=0n-l×μ∣-----∣VP∣-----∣-----01×μ∣01×η,where P∈R(l-1)×μ collects in columns the corresponding basis for Null(L) and V∈R(n-1)×η verifies that M0V=0. Note that the number of columns of V and the dimension of Im(L)∩Im(S) coincide.Theorem 15.
Letlnf(0;χ)∉Im(ST) and subspaces Null(L) and Im(L)∩Im(S) be not necessarily empty. Then, there exists a set:(46)Z=z∈R>0l-1∣RTQz,τ;χ∗=0such that for each z∈Z, there is at least one vector χ∗∈U, a real number τ∈R with sgn⁡τ=sgn⁡(Θχ∗(z)), and vectors ξ∈Rm, which satisfy (35)-(37). Moreover, as in Theorem 14, any c=exp(ξ) satisfying (35)-(37) is an equilibrium solution of the kinetic system (1).Proof.
Thatξ solves (35) for each z∈R>0l-1 follows from expression (39) in Theorem 14. Note that it also solves (37) for any τ∈R if y∉Im(L)⊕Im(S), or for a τ satisfying (40), otherwise.
However, because subspacesNull(L) and Im(L)∩Im(S) are not necessarily empty, expression (36) does not hold for any z∈R>0l-1 but for a subset related to a basis for the kernel of M0. In this way, the proof reduces to show that such subset is in fact the one defined in (46). By the rank-nullity theorem, we have that(47)dim⁡NullM0=dim⁡NullS+dim⁡NullL+dim⁡Im⁡L∩Im⁡S=δ+μ+η.This means that the basis contains some μ+η extra elements, in addition to the δ ones of the form wqT=gqT0lT, used in the proof of Theorem 14. Note that such elements correspond with the columns of matrix R (45). Since ξ solving (38) is equivalent to Q(z,τ;χ∗)∈Im(MT), the relations below must hold for each z: (48)wqTQz,τ;χ∗=0,withq=1,…,δRTQz,τ;χ∗=0which implies that z∈Z. The remaining part of the proof follows the same steps as in Theorem 14.Corollary 16.
Iflnf(0;χ)∈Im(ST), the set of equilibrium solutions is given by c=exp(ξ), with ξ satisfying (35) and f(z;χ∗)≡f(0;χ), for every χ∈U.Proof.
Iflnf(0;χ)∈Im(ST), then it clearly satisfies (35) for some ξ∈Rm. That c=exp(ξ) defines the set of equilibrium solutions can be shown by noting that f(0;χ) implies for each λ=1,…,l and χ∈U that fλ(0;χ)=fλ∗ (19). Substituting this term in (2) we get that ϕλ(c)=0, thus making the right hand side of (1) equal to zero.Remark 17.
Note that because of Lemma11, the set of equilibrium solutions described by (35)-(37) is identical to that obtained with -χ∗It should be noted that Theorem14 allows us to parameterize the set of equilibrium solutions of (1) in terms of z∈R>0l-1, including χ∗∈U as a function of z, and τ. If in addition y∈Im(L)⊕Im(S), then a functional relationship (40) is established between τ and z.This follows from the results in AppendixB that essentially shows that, for every z, there exists at least one γ satisfying the feasibility condition (B.9) (Proposition B.2). γ relates to χ∗ as in Lemma 12; i.e., χ∗=γ-1γ.Following the notation used in this paper, the complex balance condition in CRNT [10] is verified whenever lnf(0;χ)∈Im(ST), which implies that ϕ(c)=0. As a consequence of this condition, the set of equilibrium solutions of (1) can be obtained from (35) and expressed as a function of (m-s) state variables. Under these circumstances, when the kinetic system (1) describes a chemical reaction network dynamics, there is only one equilibrium per compatibility class.We end up this section with a result that links the dimension of the null space ofS in system (1) with a particular type of bifurcation (folding or limit point bifurcation) as the mechanism underlying multiple equilibrium solutions on certain compatibility classes (11).As mentioned in the introduction, there are two main bifurcation phenomena linked to multiple equilibria in dynamic systems: One is the limit point bifurcation or folding of the equilibrium manifold. The second is the branching phenomenon. Before introducing the main result of this paper, we clarify what is meant here by folding and branching of the equilibrium. In Figure2 we provide bifurcation diagrams corresponding to each of the bifurcation phenomena.Figure 2
Two types of bifurcations leading to multistationarity. (a) Folding: multistationarity arises due to a folding of the equilibrium manifold. As it is seen in the figure, the multistability region is enclosed by two folding (or limit point) bifurcations with respect to a mass conservation constant. Importantly, the three equilibria indicated in the figure by red squares correspond to different vectorsz. (b) Branching: multistationarity arises via equilibrium branching; in this case the three equilibria indicated by red squares correspond to the same vector z (there are three different values of ξ for a given z).
(a) (b)The appearance of multiple equilibria through equilibrium folding is illustrated in Figure2(a). At the folding or limit point bifurcation, the equilibrium manifold is tangent to a stoichiometric compatibility class. Note that this corresponds to the object referred to in the literature as fold, saddle node, turning point, or limit point bifurcation [8, 31, 32] with bifurcation parameter being a mass conservation constant.In the case of the equilibrium branching, illustrated in Figure2(b), different equilibrium branches emerge from a branching point, in which two equilibrium curves with different tangents intersect [8].Importantly, if multiple equilibria appear via branching, this implies that for a vectorz there are multiple values of ξ that solve simultaneously (35)-(37). In case of uniqueness of solution for system (35)-(37) no branching is possible.Proposition 18.
Let the dimension of null space ofS be δ=1. Then, for each z∈R>0l-1 there can be at most one equilibrium solution in the interior of any compatibility class Ω(c0).Proof.
Forδ=1, the set U contains only two elements χ=+1 and χ=-1, which according to Lemma 11 lead to the same set of equilibrium solutions given by (35)-(37) (see also Remark 17).
Letc∗,c∗∗ be two equilibrium solutions that belong to the interior or Ω(c0), for some c0; then we have that BT(c∗-c∗∗)=0, which implies that (c∗-c∗∗)∈Im(S). In addition, suppose that there exists one vector z, such that both ξ∗=lnc∗ and ξ∗∗=lnc∗∗ satisfy (35)-(37). If that occurs, then we have that(49)STξ∗-ξ∗∗=0and therefore (c∗-c∗∗)T(lnc∗-lnc∗∗)=0. But the above equality holds if and only if c∗=c∗∗ (see, e.g., [28]), so for each z there can be at most one equilibrium solution in Ω(c0).Note that this result implies that for deficiency-one networks, equilibrium branching as illustrated in Figure2 is not possible. If multistationarity appears, it is therefore through a folding of the equilibrium manifold.For deficiencies larger than one, the above result does not hold any more. In fact, as we illustrate next on one example, it is possible to find multiple (not aligned) vectorsχ∗∈U, satisfying the conditions of Theorem 14.


### 3.2. Some Illustrative Examples

Example 1 (positive equilibrium set for aδ=1 polynomial system).
Let us consider the system discussed in Section2 (Figure 1). The family of solutions (22) consists of the three vector functions:(50)f1x1;χ=αβ1β+x1·γα+ββγβ,(51)f2x2;χ=2+x2·0,(52)f3x3;χ=20·11-x3·21,where we have used χ=+1, since δ=1. Note that by Lemma 11, fλ(xλ;+1)=fλ(xλ;-1). Θχ(z) in expression (34) becomes(53)Θ+1z=20·β-α2·β+γ·α+β·z1.As stated in Lemma 11, Θ+1(z) is identically zero for parameters satisfying α=20·β, which implies that lnf(0;χ)∈Im(ST). Otherwise, for α≠20·β, its sign (either positive or negative depending on the selected parameters) is preserved for every positive z1, z2, which implies that lnf(0;χ)∉Im(ST).
For this exampledim(Null(L))=0, but both the images of L and S overlap, resulting in dim(Im(L)∩Im(S))=1. We then have that R (45) becomes(54)R=0100-1010T,
Theorem15 applies, and the set of equilibrium solutions is given by expressions (35)-(37). This is valid for every z∈Z, where the set Z takes the following form.(55)Z=z∈R>0l-1∣z2=1+γz1Θ+1zβ20-Θ+1zSince y∈Im(L)⊕Im(S), τ is a function of z of the form (40), with(56)Ψz,+1=βz11+γz1Θ+1z.The set of all possible equilibrium states, as a function of z, is given by the following expression:(57)lnc=0000-10-10000-101-10-211-11-1000·lnf1z1Θ+1z1;+1ln⁡2ln⁡20-Θ+1z1-ln⁡z1.By Proposition18, for each z∈Z, which in this case is equivalent to each z1∈R>0, one can expect one equilibrium at most. Note that for lnf(0;+1)∈Im(ST) there can be only one equilibrium in each compatibility class. Otherwise, we cannot preclude the existence of multiple steady states in the same compatibility class. For the present example, these two scenarios are illustrated in Figure 3. On the one hand, for α=100,β=100, γ=10, there is a unique equilibrium for every compatibility class. On the other hand, for α=1,β=100, γ=10, the network has three equilibria for values of the mass conservation constant within a certain range (including σ=500). Multistationarity arises via equilibrium folding, in agreement with the results presented in the previous section for deficiency-one networks.Figure 3
Mass conservation constantσ=BTc versus z for two sets of parameters (Example 1). The continuous line (α=1, β=100, γ=10) shows three equilibria (stars) for the compatibility class at σ=500 (dotted line). The dashed line (α=100, β=100, γ=10) shows unique solutions in every compatibility class.Example 2 (kinetic system (ofδ=2) exhibiting branching).
Let us consider the weakly reversible reaction network taken from [12] in Figure 4, with reaction constants k12=k34=α and k23=k41=β. The dynamics of this l=1 network can be accommodated in the structure given by (1) with ϕ(c)=ψ1(c)a+Ef(c), ψ1(c)=exp(yTlnc) (see Remark 3) and(58)E=-β00β-α00α-β,a=α00,b=00β,y=30.For this example, we have that(59)S=-2-3-1231andG=1001-2-3,where the columns of G define a basis of the null space of S. This matrix can be partitioned as in (B.1) with GA∈R2×2, invertible. The family of solutions (24) becomes(60)fx,χ=αβ1αβ+x·-1β0-1α-1α1β2β·χ.Note that by Corollary 13, for each set of parameters (α,β) there must be at least one χ∗∈U for which f(x,χ) is feasible. In order to compute those unit vectors, we make use of the arguments in Proposition B.2 to find a vector γ≡x·χ that makes (60) feasible. Using (B.9) and (B.12), we get (61)φϵ=ϵ2ϵ3andγϵ=-βϵ2+α-αϵ3+βϵ2.Substituting the above expressions in (B.5) and reordering terms we finally get(62)ϵ-ϵ2+2αβϵ3-2αβ=0.Each solution ϵ∗ of (62) corresponds with a given γ∗ and in turn with a given χ∗=γ∗-1γ∗ (Corollary 13). Solutions of (62) for different values of the parameter ratio α/β are presented in Figure 5. For each parameter set below the line α=0.17β, system (62) has 3 positive real solutions corresponding with 3 different χ∗, and leading to 3 branches. As shown in Figure 6, this results in multiple equilibria in the form of bistability for each compatibility class. Compatibility classes are given by (11) with B=11T.Figure 4
Example2: weakly reversible network with l=1 and N=4.Figure 5
Example2: positive solutions of system (62) as a function of the parameter ratio α/β.Figure 6
Example2: equilibria (solid lines) for α=0.16, β=1 in different compatibility classes (dashed lines).


## 3.1. Characterizing the Equilibrium Set

Depending on the dimensions of subspacesIm⁡(L)∩Im⁡(S) and Null(L) (Definition 2), two scenarios will be considered. In the first place, we assume that both are empty and get the following result.Theorem 14.
Let the dimensionsη=μ=0 (Definition 2), and lnf(0;χ)∉Im⁡(ST). Then, for each z∈R>0l-1 there is at least one vector χ∗∈U, a real number τ∈R with sgn⁡(τ)=sgn⁡(Θχ∗(z)), and vectors ξ∈Rm, such that(35)STξ=lnfz;χ∗(36)LTξ=-lnz(37)yTξ=lnτΘχ∗z.Moreover, any c=exp(ξ) satisfying (35)-(37) is an equilibrium solution of the kinetic system (1).Proof.
In order to prove thatξ simultaneously solves (35)-(37), let us define the matrix M=M0∣y, with M0=S∣L, and rewrite these equations as(38)MTξ=Qz,τ;χ∗,withQz,τ;χ∗=lnfz;χ∗-lnzlnτΘχ∗z.Suppose firstly that y∉Im(L)⊕Im(S). Then, a basis for the kernel of M consists of a vector set wq∣q=1,…,δ, where each element is of the form wqT=gqT0lT, with 0l being the l dimensional vector zero. Note that each element is in the null space by construction since the corresponding gq is in the null space of S, and the δ elements are linearly independent.
In addition, because of dimensionsη=μ=0, no other vector linearly independent to the above defined vector set can be found satisfying that Mw=0. If such a vector w∈Rn, with wT=wSTwLTwyT, and wS,wL,wy≠0 exists, then SwS+LwL+ywy=0. But such a relation would require the overlapping of the corresponding images of S, L, and y, which is against the assumption that the dimensions of subspaces Im(L)∩Im(S) and Null(L) are zero (η=μ=0).
Thatξ solves (38) is equivalent to saying that Q(z,τ;χ∗)∈Im(MT) for each vector z∈R>0l-1 and τ∈R, with sgn(τ)=sgn(Θχ∗(z)), which in turn requires that wqTQ(z,τ;χ∗)=0, for every q=1,…,δ.
Becauselnf(0;χ)∉Im(ST), by Lemma 11 we have that, for every χ∈U, Θχ(z) in (34) has a positive (respectively, negative) sign. From Lemma 12, it follows that for each vector z∈R>0l-1 there is at least one vector χ∗∈U such that f(Xχ∗(z);χ∗) is feasible (with Xχ∗(z) being a vector function of the form (34)) and we get expression (35).
Thatξ solves (38) for each z∈R>0l-1 and any τ∈R is a consequence of the structure of each vector of the basis, so that(39)wqTQz,τ;χ∗=gqTlnfz;χ∗=0.Suppose, on the other hand, that y∈Im(L)⊕Im(S). Then y can be expressed as a linear combination of the columns of L and S and the dimension of the kernel of M becomes δ+1. The new element of the basis can be expressed as a vector wδ+1T=wSTwLT1. In order for Q(z,τ;χ∗) to satisfy (38) we have that wδ+1TQ(z,τ;χ∗)=0. This equation can be solved in τ resulting in a function:(40)τz;χ∗=Θχ∗zΨz,χ∗, with Ψ(z,χ∗)=exp⁡wSTlnf(z;χ∗)-wLTln⁡z. Consequently, for each z and χ∗ there exists one τ and vectors ξ∈Rm satisfying (38).
In order to prove that anyc=exp(ξ) satisfying (35)-(37) is in fact an equilibrium solution, we combine the above expressions with (3)-(4) and make use of (19) to get(41)ψrc≡τΘχ∗zz-11,(42)fλc≡Eλ-1-aλ+zλΘχ∗zGλχ∗,λ≠l(43)flc≡El-1-al+Θχ∗zGlχ∗,where z-1 in (42) should be understood as the element-wise inverse of vector z. Combining (41)-(43) with (2), we also obtain that ϕλ(c)≡τGλχ∗ for λ=1,…,l. Finally, substituting previous expressions in (1), we get(44)c˙=τ∑λSλGλχ∗=0,which proves that any c=exp(ξ) satisfying (35)-(37) is in fact an equilibrium solution of the kinetic system (1).Ifη or μ in Definition 2 is positive, then Theorem 14 does not hold for all z∈R>0l-1, but for a subset of it. In preparing for this result, let us define the matrix:(45)R=0n-l×μ∣-----∣VP∣-----∣-----01×μ∣01×η,where P∈R(l-1)×μ collects in columns the corresponding basis for Null(L) and V∈R(n-1)×η verifies that M0V=0. Note that the number of columns of V and the dimension of Im(L)∩Im(S) coincide.Theorem 15.
Letlnf(0;χ)∉Im(ST) and subspaces Null(L) and Im(L)∩Im(S) be not necessarily empty. Then, there exists a set:(46)Z=z∈R>0l-1∣RTQz,τ;χ∗=0such that for each z∈Z, there is at least one vector χ∗∈U, a real number τ∈R with sgn⁡τ=sgn⁡(Θχ∗(z)), and vectors ξ∈Rm, which satisfy (35)-(37). Moreover, as in Theorem 14, any c=exp(ξ) satisfying (35)-(37) is an equilibrium solution of the kinetic system (1).Proof.
Thatξ solves (35) for each z∈R>0l-1 follows from expression (39) in Theorem 14. Note that it also solves (37) for any τ∈R if y∉Im(L)⊕Im(S), or for a τ satisfying (40), otherwise.
However, because subspacesNull(L) and Im(L)∩Im(S) are not necessarily empty, expression (36) does not hold for any z∈R>0l-1 but for a subset related to a basis for the kernel of M0. In this way, the proof reduces to show that such subset is in fact the one defined in (46). By the rank-nullity theorem, we have that(47)dim⁡NullM0=dim⁡NullS+dim⁡NullL+dim⁡Im⁡L∩Im⁡S=δ+μ+η.This means that the basis contains some μ+η extra elements, in addition to the δ ones of the form wqT=gqT0lT, used in the proof of Theorem 14. Note that such elements correspond with the columns of matrix R (45). Since ξ solving (38) is equivalent to Q(z,τ;χ∗)∈Im(MT), the relations below must hold for each z: (48)wqTQz,τ;χ∗=0,withq=1,…,δRTQz,τ;χ∗=0which implies that z∈Z. The remaining part of the proof follows the same steps as in Theorem 14.Corollary 16.
Iflnf(0;χ)∈Im(ST), the set of equilibrium solutions is given by c=exp(ξ), with ξ satisfying (35) and f(z;χ∗)≡f(0;χ), for every χ∈U.Proof.
Iflnf(0;χ)∈Im(ST), then it clearly satisfies (35) for some ξ∈Rm. That c=exp(ξ) defines the set of equilibrium solutions can be shown by noting that f(0;χ) implies for each λ=1,…,l and χ∈U that fλ(0;χ)=fλ∗ (19). Substituting this term in (2) we get that ϕλ(c)=0, thus making the right hand side of (1) equal to zero.Remark 17.
Note that because of Lemma11, the set of equilibrium solutions described by (35)-(37) is identical to that obtained with -χ∗It should be noted that Theorem14 allows us to parameterize the set of equilibrium solutions of (1) in terms of z∈R>0l-1, including χ∗∈U as a function of z, and τ. If in addition y∈Im(L)⊕Im(S), then a functional relationship (40) is established between τ and z.This follows from the results in AppendixB that essentially shows that, for every z, there exists at least one γ satisfying the feasibility condition (B.9) (Proposition B.2). γ relates to χ∗ as in Lemma 12; i.e., χ∗=γ-1γ.Following the notation used in this paper, the complex balance condition in CRNT [10] is verified whenever lnf(0;χ)∈Im(ST), which implies that ϕ(c)=0. As a consequence of this condition, the set of equilibrium solutions of (1) can be obtained from (35) and expressed as a function of (m-s) state variables. Under these circumstances, when the kinetic system (1) describes a chemical reaction network dynamics, there is only one equilibrium per compatibility class.We end up this section with a result that links the dimension of the null space ofS in system (1) with a particular type of bifurcation (folding or limit point bifurcation) as the mechanism underlying multiple equilibrium solutions on certain compatibility classes (11).As mentioned in the introduction, there are two main bifurcation phenomena linked to multiple equilibria in dynamic systems: One is the limit point bifurcation or folding of the equilibrium manifold. The second is the branching phenomenon. Before introducing the main result of this paper, we clarify what is meant here by folding and branching of the equilibrium. In Figure2 we provide bifurcation diagrams corresponding to each of the bifurcation phenomena.Figure 2
Two types of bifurcations leading to multistationarity. (a) Folding: multistationarity arises due to a folding of the equilibrium manifold. As it is seen in the figure, the multistability region is enclosed by two folding (or limit point) bifurcations with respect to a mass conservation constant. Importantly, the three equilibria indicated in the figure by red squares correspond to different vectorsz. (b) Branching: multistationarity arises via equilibrium branching; in this case the three equilibria indicated by red squares correspond to the same vector z (there are three different values of ξ for a given z).
(a) (b)The appearance of multiple equilibria through equilibrium folding is illustrated in Figure2(a). At the folding or limit point bifurcation, the equilibrium manifold is tangent to a stoichiometric compatibility class. Note that this corresponds to the object referred to in the literature as fold, saddle node, turning point, or limit point bifurcation [8, 31, 32] with bifurcation parameter being a mass conservation constant.In the case of the equilibrium branching, illustrated in Figure2(b), different equilibrium branches emerge from a branching point, in which two equilibrium curves with different tangents intersect [8].Importantly, if multiple equilibria appear via branching, this implies that for a vectorz there are multiple values of ξ that solve simultaneously (35)-(37). In case of uniqueness of solution for system (35)-(37) no branching is possible.Proposition 18.
Let the dimension of null space ofS be δ=1. Then, for each z∈R>0l-1 there can be at most one equilibrium solution in the interior of any compatibility class Ω(c0).Proof.
Forδ=1, the set U contains only two elements χ=+1 and χ=-1, which according to Lemma 11 lead to the same set of equilibrium solutions given by (35)-(37) (see also Remark 17).
Letc∗,c∗∗ be two equilibrium solutions that belong to the interior or Ω(c0), for some c0; then we have that BT(c∗-c∗∗)=0, which implies that (c∗-c∗∗)∈Im(S). In addition, suppose that there exists one vector z, such that both ξ∗=lnc∗ and ξ∗∗=lnc∗∗ satisfy (35)-(37). If that occurs, then we have that(49)STξ∗-ξ∗∗=0and therefore (c∗-c∗∗)T(lnc∗-lnc∗∗)=0. But the above equality holds if and only if c∗=c∗∗ (see, e.g., [28]), so for each z there can be at most one equilibrium solution in Ω(c0).Note that this result implies that for deficiency-one networks, equilibrium branching as illustrated in Figure2 is not possible. If multistationarity appears, it is therefore through a folding of the equilibrium manifold.For deficiencies larger than one, the above result does not hold any more. In fact, as we illustrate next on one example, it is possible to find multiple (not aligned) vectorsχ∗∈U, satisfying the conditions of Theorem 14.


## 3.2. Some Illustrative Examples

Example 1 (positive equilibrium set for aδ=1 polynomial system).
Let us consider the system discussed in Section2 (Figure 1). The family of solutions (22) consists of the three vector functions:(50)f1x1;χ=αβ1β+x1·γα+ββγβ,(51)f2x2;χ=2+x2·0,(52)f3x3;χ=20·11-x3·21,where we have used χ=+1, since δ=1. Note that by Lemma 11, fλ(xλ;+1)=fλ(xλ;-1). Θχ(z) in expression (34) becomes(53)Θ+1z=20·β-α2·β+γ·α+β·z1.As stated in Lemma 11, Θ+1(z) is identically zero for parameters satisfying α=20·β, which implies that lnf(0;χ)∈Im(ST). Otherwise, for α≠20·β, its sign (either positive or negative depending on the selected parameters) is preserved for every positive z1, z2, which implies that lnf(0;χ)∉Im(ST).
For this exampledim(Null(L))=0, but both the images of L and S overlap, resulting in dim(Im(L)∩Im(S))=1. We then have that R (45) becomes(54)R=0100-1010T,
Theorem15 applies, and the set of equilibrium solutions is given by expressions (35)-(37). This is valid for every z∈Z, where the set Z takes the following form.(55)Z=z∈R>0l-1∣z2=1+γz1Θ+1zβ20-Θ+1zSince y∈Im(L)⊕Im(S), τ is a function of z of the form (40), with(56)Ψz,+1=βz11+γz1Θ+1z.The set of all possible equilibrium states, as a function of z, is given by the following expression:(57)lnc=0000-10-10000-101-10-211-11-1000·lnf1z1Θ+1z1;+1ln⁡2ln⁡20-Θ+1z1-ln⁡z1.By Proposition18, for each z∈Z, which in this case is equivalent to each z1∈R>0, one can expect one equilibrium at most. Note that for lnf(0;+1)∈Im(ST) there can be only one equilibrium in each compatibility class. Otherwise, we cannot preclude the existence of multiple steady states in the same compatibility class. For the present example, these two scenarios are illustrated in Figure 3. On the one hand, for α=100,β=100, γ=10, there is a unique equilibrium for every compatibility class. On the other hand, for α=1,β=100, γ=10, the network has three equilibria for values of the mass conservation constant within a certain range (including σ=500). Multistationarity arises via equilibrium folding, in agreement with the results presented in the previous section for deficiency-one networks.Figure 3
Mass conservation constantσ=BTc versus z for two sets of parameters (Example 1). The continuous line (α=1, β=100, γ=10) shows three equilibria (stars) for the compatibility class at σ=500 (dotted line). The dashed line (α=100, β=100, γ=10) shows unique solutions in every compatibility class.Example 2 (kinetic system (ofδ=2) exhibiting branching).
Let us consider the weakly reversible reaction network taken from [12] in Figure 4, with reaction constants k12=k34=α and k23=k41=β. The dynamics of this l=1 network can be accommodated in the structure given by (1) with ϕ(c)=ψ1(c)a+Ef(c), ψ1(c)=exp(yTlnc) (see Remark 3) and(58)E=-β00β-α00α-β,a=α00,b=00β,y=30.For this example, we have that(59)S=-2-3-1231andG=1001-2-3,where the columns of G define a basis of the null space of S. This matrix can be partitioned as in (B.1) with GA∈R2×2, invertible. The family of solutions (24) becomes(60)fx,χ=αβ1αβ+x·-1β0-1α-1α1β2β·χ.Note that by Corollary 13, for each set of parameters (α,β) there must be at least one χ∗∈U for which f(x,χ) is feasible. In order to compute those unit vectors, we make use of the arguments in Proposition B.2 to find a vector γ≡x·χ that makes (60) feasible. Using (B.9) and (B.12), we get (61)φϵ=ϵ2ϵ3andγϵ=-βϵ2+α-αϵ3+βϵ2.Substituting the above expressions in (B.5) and reordering terms we finally get(62)ϵ-ϵ2+2αβϵ3-2αβ=0.Each solution ϵ∗ of (62) corresponds with a given γ∗ and in turn with a given χ∗=γ∗-1γ∗ (Corollary 13). Solutions of (62) for different values of the parameter ratio α/β are presented in Figure 5. For each parameter set below the line α=0.17β, system (62) has 3 positive real solutions corresponding with 3 different χ∗, and leading to 3 branches. As shown in Figure 6, this results in multiple equilibria in the form of bistability for each compatibility class. Compatibility classes are given by (11) with B=11T.Figure 4
Example2: weakly reversible network with l=1 and N=4.Figure 5
Example2: positive solutions of system (62) as a function of the parameter ratio α/β.Figure 6
Example2: equilibria (solid lines) for α=0.16, β=1 in different compatibility classes (dashed lines).


## 4. Conclusions

We derived an explicit representation of the set of all positive equilibria attained by a class of polynomial dynamic systems. This result is formally stated in Theorem14. Such polynomial class turns out to accommodate in its state-space form a large class of weakly reversible networks with mass-action kinetics, which allowed us to explore in a straightforward manner the different mechanisms leading to multiple equilibria.Remarkably, network deficiency, as defined in Chemical Reaction Network Theory, plays a central role in the class of potential bifurcations. As stated in Proposition18, multiple equilibria on deficiency-one networks can only arise via folding of the equilibrium manifold. Bifurcation phenomena leading to multiple branches, on the other hand, demand higher deficiencies.The methods discussed here will hopefully pave the way to explore more complex classes of bifurcations linked to limit cycles or even chaotic attractors, for instance. In particular, Theorem14, combined with appropriate stability criteria, suggests a number of research directions oriented to the study of the effect of network structure and parameters on the resulting dynamic behavior, which can be extended to the analysis of its qualitative behavior under uncertainty. In the context of process systems, the results have relevant implications in the robust control of complex chemical reactors and the design of biochemical circuits.



---

*Source: 3912627-2018-12-12.xml*