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.
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:where the state vector is nonnegative and usually denotes species concentrations. The summation extends over , and for each we have and of the form:where are a class of invertible Metzler matrices [29], are nonnegative vectors, and both contain system parameters.
Assumption 1. In this contribution, we specialize on pairs such that the resulting vector (i.e., strictly positive). (As it is extensively discussed in [28],any weakly reversible chemical reaction network satisfies this assumption.)
Each function in expression (2) is a component of the vector function , which together with takes the form:where , is a (positive) coefficient vector, and . We refer to (1) as a dynamic polynomial system, which can be written in compact form:where (with ) reads and is of the form with subvectors as in (2).
Definition 2. We consider two subspaces associated with matrices in (3) and in (4)-(5) with the following dimensions:
Remark 3. When , expression (3) reduces to
Although not restricted to them, for and appropriate matrices and , 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.
2.2. Invariance and Equilibrium
All possible equilibrium solutions of system (5) necessarily correspond to either being a zero vector or lying in the null space of . Let ; then by the rank-nullity theorem, the dimension of the null space of is . Such number is known in CRNT as the deficiency of the chemical reaction network, so we will refer to the null space of as the deficiency subspace. Let be a basis of the null space of , and express each vector in terms of subvectors such thatFrom the above description, the equality , for every , can be expressed as
Note that by using matrices of the form:expression (9) can be rewritten as . On occasions, we may use the matrix .
Assumption 4. We assume in what follows that system (1) is such that is full rank.
Starting from an arbitrary reference concentration vector , we define a setwhere the columns of are a basis of . Note that is an invariant of motion for system (1), since , so that any trajectory that starts in 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 is known as the stoichiometric compatibility class at .
2.3. Representing CRN Dynamics
A CRN involving 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 , , , , and 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 and 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.Note that the third column is the zero vector (the complex corresponds to the environment). Let be the rate at which a set of reactants in complex is transformed into a set of products in complex , and assume that it obeys the mass-action law so thatUnder constant volume and flux and taking into account the fact that only and are exchanged with the environment, balance equations for the five chemical species become where 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 Figure 1 correspond with columns , , and and are identified as , , and , respectively. They are employed to construct the stoichiometric subspace as discussed in [28], thus leading to (1), where the matrices are of the form:Vectors and matrices in (2) are as follows:where , represent the input product and substrate concentrations, respectively, and the residence time of the system. Matrices , , and are C-Metzler (see Lemma A.3 in the appendix) so their inverses are nonpositive:As in Assumption 1, the network is reversible (thus weakly reversible) and satisfies that , for .
We obtain expression in (3) by constructing and from the reference complexes as follows:where and . 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 . Since the rank of in (5) is and , , and so that , the dimension of the null space is . A basis for the null space is given by the vector , so that , , and .
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 , a vector function aswhere and , with being any element of the set:For any , there exists an interval where remains positive. Note that, since by Assumption 1, the intervalis never empty and contains the zero. We define an interval and consider the vector function aswith . Note that for every we have that (i.e., strictly positive). For convenience, let us define the following matrices:to rewrite the vector function (22) as
2.4.1. A Class of Strictly Monotonous Functions
Let us consider an open interval on the real line and define a vector function aswhere is C-Metzler (thus invertible, see Appendix A), and , with and with at least one positive component. Moreover, let , such that the values of the scalar for which the vector function remains strictly positive lie in the interval , which includes the zero. Note that since , . Moreover, at least one of the extremes is lower (respectively, upper) bounded, depending on whether the nonzero elements of have equal or distinct signs (see proofs of Proposition 4.1 and Theorem 5.1 in [28]).
Associated with , we have a scalar function of the form: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]). in (26) is monotonous decreasing with a continuous lower bounded derivative on the domain . Moreover,
Remark 6. One main consequence of Lemma 5 is that defines a bijective map between the set (function domain) and its range (codomain) .
For and , let us consider nonzero vectors , C-Metzler matrices , and nonnegative vectors , such that , and construct the set of functions of the form (26). In addition, let us consider a domain and construct a function aswith . First note that since for every , . From Lemma 5 (in particular Remark 6), we can always find elements such that . From Lemma 5 it also follows that each has a continuous negative derivative. Hence, we can conclude that satisfies the conditions of the implicit function theorem (see, e.g., [30]).
In this way, for every , there exists a domain , a vector , and a well-defined Lipschitz continuous function , such that . Moreover, since is monotonous decreasing and its derivative is lower bounded, we have thatwhich implies that each has negative (and lower bounded) partial derivatives. Before presenting the result below, let us define for convenience the set (respectively, ).
Theorem 7. Let be a strictly positive vector and (respectively, ). Then, there exists a function (respectively, ), well defined and Lipschitz continuous, of the form:where (respectively, ), that solves the following system.
Proof. Since , there are vectors in that solve . To prove the above assertion, note that satisfies . By assumption , because of (29), as each increases above zero, decreases below . Thus, repeating this for every , by continuity one finds such positive vectors . A similar argument can be employed for . In fact, such vectors are also solutions of system (31) for , with .
Using a vector function to represent the left hand side of (31), we then have that for some and . In addition, is Lipschitz continuously differentiable everywhere in and . Moreover,is nonsingular. Hence the implicit function theorem [30] ensures the existence of a function well defined and Lipschitz continuous, such that . Let ; then from (31) we have for each that , which leads to expression (30).
Corollary 8. If , then is the identically zero function.
Proof. For , we have that , so for every , is the only solution of (31), and therefore .
Remark 9. For , (28) reduces to , and because of Lemma 5, there exists only one such that .
2.4.2. The Feasibility Condition
Provided that the logarithm of in (22) lies in the image of , any in (19) is compatible with (4). We then say that the family expressed in (22) is feasible.
Lemma 10 (feasibility condition). in (22) is feasible for any and such that , for every .
Proof. Since feasibility implies that is in the image of , it must be orthogonal to all elements of the basis , which proves the statement.
In preparing to explore the consequences of feasibility, for each we define a function , of the form:with , and . The following lemma describes possible solutions for .
Lemma 11. For every , there exists a well-defined and Lipschitz continuous function such thatwith , that solves . Moreover, is strictly positive (respectively, negative) if and only if . Otherwise, it is the identically zero function. In addition, .
Proof. implies that , so by Assumption 4, we can find a (respectively, ) so that Theorem 7 applies, and expression (34) follows with strictly positive (respectively, negative).
If , we necessarily have that , so by Corollary 8 we conclude that is identically zero.
Concerning the final part, note that for every we have that , so the left hand side of system (31) in the proof of Theorem 7 becomes , which leads to . From (34), then we have that . Substituting the above equivalences in (22) we get that which concludes the proof.
Lemma 12. For every , there exists at least one vector for which is feasible.
Proof. By Proposition B.2, we have that for every there exists a vector (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 , with .
Corollary 13. For there is at least one vector for which (24) is feasible.
Proof. If , then for every there exists one scalar that solves (Remark 9). Thus for , 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 subspaces and (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 (Definition 2), and . Then, for each there is at least one vector , a real number with , and vectors , such thatMoreover, any 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 , with , and rewrite these equations asSuppose firstly that . Then, a basis for the kernel of consists of a vector set , where each element is of the form , with being the dimensional vector zero. Note that each element is in the null space by construction since the corresponding is in the null space of , and the elements are linearly independent.
In addition, because of dimensions , no other vector linearly independent to the above defined vector set can be found satisfying that . If such a vector , with , and exists, then . But such a relation would require the overlapping of the corresponding images of , , and , which is against the assumption that the dimensions of subspaces and are zero ().
That solves (38) is equivalent to saying that for each vector and , with , which in turn requires that , for every .
Because , by Lemma 11 we have that, for every , in (34) has a positive (respectively, negative) sign. From Lemma 12, it follows that for each vector there is at least one vector such that is feasible (with being a vector function of the form (34)) and we get expression (35).
That solves (38) for each and any is a consequence of the structure of each vector of the basis, so thatSuppose, on the other hand, that . Then can be expressed as a linear combination of the columns of and and the dimension of the kernel of becomes . The new element of the basis can be expressed as a vector . In order for to satisfy (38) we have that . This equation can be solved in resulting in a function: with . Consequently, for each and there exists one and vectors satisfying (38).
In order to prove that any satisfying (35)-(37) is in fact an equilibrium solution, we combine the above expressions with (3)-(4) and make use of (19) to getwhere in (42) should be understood as the element-wise inverse of vector . Combining (41)-(43) with (2), we also obtain that for . Finally, substituting previous expressions in (1), we getwhich proves that any 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 , but for a subset of it. In preparing for this result, let us define the matrix:where collects in columns the corresponding basis for and verifies that . Note that the number of columns of and the dimension of coincide.
Theorem 15. Let and subspaces and be not necessarily empty. Then, there exists a set:such that for each , there is at least one vector , a real number with , and vectors , which satisfy (35)-(37). Moreover, as in Theorem 14, any satisfying (35)-(37) is an equilibrium solution of the kinetic system (1).
Proof. That solves (35) for each follows from expression (39) in Theorem 14. Note that it also solves (37) for any if , or for a satisfying (40), otherwise.
However, because subspaces and are not necessarily empty, expression (36) does not hold for any but for a subset related to a basis for the kernel of . 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 thatThis means that the basis contains some extra elements, in addition to the ones of the form , used in the proof of Theorem 14. Note that such elements correspond with the columns of matrix (45). Since solving (38) is equivalent to , the relations below must hold for each : which implies that . The remaining part of the proof follows the same steps as in Theorem 14.
Corollary 16. If , the set of equilibrium solutions is given by , with satisfying (35) and , for every .
Proof. If , then it clearly satisfies (35) for some . That defines the set of equilibrium solutions can be shown by noting that implies for each and that (19). Substituting this term in (2) we get that , thus making the right hand side of (1) equal to zero.
Remark 17. Note that because of Lemma 11, the set of equilibrium solutions described by (35)-(37) is identical to that obtained with
It should be noted that Theorem 14 allows us to parameterize the set of equilibrium solutions of (1) in terms of , including as a function of , and . If in addition , then a functional relationship (40) is established between and .
This follows from the results in Appendix B that essentially shows that, for every , there exists at least one satisfying the feasibility condition (B.9) (Proposition B.2). relates to as in Lemma 12; i.e., .
Following the notation used in this paper, the complex balance condition in CRNT [10] is verified whenever , which implies that . As a consequence of this condition, the set of equilibrium solutions of (1) can be obtained from (35) and expressed as a function of 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 of 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 Figure 2 we provide bifurcation diagrams corresponding to each of the bifurcation phenomena.
(a)
(b)
The appearance of multiple equilibria through equilibrium folding is illustrated in Figure 2(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 Figure 2(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 vector 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 of be . Then, for each there can be at most one equilibrium solution in the interior of any compatibility class .
Proof. For , the set contains only two elements and , which according to Lemma 11 lead to the same set of equilibrium solutions given by (35)-(37) (see also Remark 17).
Let be two equilibrium solutions that belong to the interior or , for some ; then we have that , which implies that . In addition, suppose that there exists one vector , such that both and satisfy (35)-(37). If that occurs, then we have thatand therefore . But the above equality holds if and only if (see, e.g., [28]), so for each there can be at most one equilibrium solution in .
Note that this result implies that for deficiency-one networks, equilibrium branching as illustrated in Figure 2 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 , satisfying the conditions of Theorem 14.
3.2. Some Illustrative Examples
Example 1 (positive equilibrium set for a polynomial system). Let us consider the system discussed in Section 2 (Figure 1). The family of solutions (22) consists of the three vector functions:where we have used , since . Note that by Lemma 11, . in expression (34) becomesAs stated in Lemma 11, is identically zero for parameters satisfying , which implies that . Otherwise, for , its sign (either positive or negative depending on the selected parameters) is preserved for every positive , , which implies that .
For this example , but both the images of and overlap, resulting in . We then have that (45) becomesTheorem 15 applies, and the set of equilibrium solutions is given by expressions (35)-(37). This is valid for every , where the set takes the following form.Since , is a function of of the form (40), withThe set of all possible equilibrium states, as a function of , is given by the following expression:
By Proposition 18, for each , which in this case is equivalent to each , one can expect one equilibrium at most. Note that for 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 , , there is a unique equilibrium for every compatibility class. On the other hand, for , , the network has three equilibria for values of the mass conservation constant within a certain range (including ). Multistationarity arises via equilibrium folding, in agreement with the results presented in the previous section for deficiency-one networks.
Example 2 (kinetic system (of ) exhibiting branching). Let us consider the weakly reversible reaction network taken from [12] in Figure 4, with reaction constants and . The dynamics of this network can be accommodated in the structure given by (1) with , (see Remark 3) andFor this example, we have thatwhere the columns of define a basis of the null space of . This matrix can be partitioned as in (B.1) with , invertible. The family of solutions (24) becomesNote that by Corollary 13, for each set of parameters there must be at least one for which is feasible. In order to compute those unit vectors, we make use of the arguments in Proposition B.2 to find a vector that makes (60) feasible. Using (B.9) and (B.12), we get Substituting the above expressions in (B.5) and reordering terms we finally getEach solution of (62) corresponds with a given and in turn with a given (Corollary 13). Solutions of (62) for different values of the parameter ratio are presented in Figure 5. For each parameter set below the line , system (62) has positive real solutions corresponding with different , and leading to 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 .
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 Theorem 14. 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 Proposition 18, 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, Theorem 14, 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.
Appendix
A. Some Classes of Compartmental Matrices
Definition A.1 (see [33]). A matrix is compartmental if (1), for , ;(2), for .
Definition A.2 (see [28]). An invertible compartmental matrix is C-Metzler if its entries are of the form:
Lemma A.3 (see [28]). The inverse of a C-Metzler matrix is nonpositive.
B. Some Results of Existence
In preparing for the results below, we partition the full rank matrix discussed in Section 2, possibly after some reordering, aswhere is invertible and (note that ). In addition, by using (34) from Lemma 11, we can rewrite (24) asApplying the same reordering as to get (B.1), we can partition the matrices at the right hand of (B.2) aswhere is invertible.
Lemma B.1. Let the systemwith invertible, , and , have a solution for a nonzero vector . Suppose that is a nonsingular matrix where and , are inverses of diagonal matrices with components of (respectively, ) in the diagonal. Then a solution exists for the system with .
Proof. Let us construct the vector function : We have that (i) , by assumption; (ii) is Lipschitz continuous and differentiable for every and ; (iii) its Jacobian with respect to is of the form: in the above can be inverted using blockwise inversion whenever (i) is a square invertible matrix and (ii) is a nonsingular matrix. Since (i) holds because is nonsingular and (ii) is consequence of the assumption that is nonsingular, then the matrix is nonsingular for every , and .
Thus, from the implicit function theorem (e.g., [30]) there exist some continuous and differentiable maps and , such that . In particular for , we have a set of vectors that solve (B.4)-(B.5). This completes the proof.
Proposition B.2. Let us consider the vector:where , , and the pair of vectors relate as , with , and . For every there exists at least one such that
Proof. From Lemma 11, for each and , the right hand side of (B.2) makes (33). Equivalently, vector (B.8) with , solves . Choosing one and a , we compute one that is substituted in (B.8) leading to a set of positive vectors satisfyingActually, suppose that ; then the pair is a solution of the system: with . Applying Lemma B.1, it follows that a solution also exists for system (B.4)-(B.5) with . Let be one such solution. The corresponding is computed from the structure of (B.8) asBy construction, for this , the pair satisfies (B.9).
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work has been partially financed with grants PIE201230E042, grant Salvador de Madariaga (PR2011-0363), and AGL2015-67504-C3-2-R. Manuel Pájaro acknowledges support from Spanish MINECO Fellowship No. BES-2013-063112 and PIE201870E041.