A Some results required to prove Theorem 5.1

# Uniqueness of feasible equilibria for mass action law (MAL) kinetic systems

## Abstract

This paper studies the relations among system parameters, uniqueness, and stability of equilibria, for kinetic systems given in the form of polynomial ODEs. Such models are commonly used to describe the dynamics of nonnegative systems, with a wide range of application fields such as chemistry, systems biology, process modeling or even transportation systems. Using a flux-based description of kinetic models, a canonical representation of the set of all possible feasible equilibria is developed.

The characterization is made in terms of strictly stable compartmental matrices to define the so-called family of solutions. Feasibility is imposed by a set of constraints, which are linear on a log-transformed space of complexes, and relate to the kernel of a matrix, the columns of which span the stoichiometric subspace. One particularly interesting representation of these constraints can be expressed in terms of a class of monotonous decreasing functions. This allows connections to be established with classical results in CRNT that relate to the existence and uniqueness of equilibria along positive stoichiometric compatibility classes.

In particular, monotonicity can be employed to identify regions in the set of possible reaction rate coefficients leading to complex balancing, and to conclude uniqueness of equilibria for a class of positive deficiency networks. The latter result might support constructing an alternative proof of the well-known deficiency one theorem. The developed notions and results are illustrated through examples.

Keywords: Chemical reaction networks, kinetic systems, mass action law, network deficiency, feasible equilibrium, complex balanced equilibrium

This manuscript was published as: A. A. Alonso and G. Szederkényi. Uniqueness of feasible equilibria for mass action law (MAL) kinetic systems. Journal of Process Control, 48: 41–71, 2016. DOI link: http://dx.doi.org/10.1016/j.jprocont.2016.10.002

## Nomenclature

 Notation Description Defining/introducing eqn. or (sub)section Rn n-dimensional real space Rn>0 (Rn<0) n-dimensional positive (resp. negative) orthant Rn≥0 (Rn≤0) n-dimensional non-negative (resp. non-positive) orthant x>0 (x<0) each element of the vector x is positive (resp. negative) x≥0 (x≤0) each element of the vector x is non-negative (resp. non-positive) 1n∈Rn The n-dimensional vector with each element being one εi the ith standard basis vector in Rn sec. 2, sec.3 D(x) diagonal matrix D(x)∈Rn×n with components of x in the diagonal eq. (28),(73) m number of species sec. 2 n number of complexes sec. 2 ρ number irreversible chemical reaction steps in the network sec. 2 Rij rate of the reaction from complex i to complex j eq. (1) kij rate coefficient of the reaction from complex i to complex j eq. (1) ℓ the number of linkage classes sec. 2 λ integer for indexing linkage classes sec. 2 Lλ the λth linkage class sec. 2 Nλ the number of complexes in linkage class no. λ sec. 2 (footnote 1) jλ index of the reference complex in linkage class no. λ subsec. 2.1 c vector of concentrations (state variables) sec. 2 Y m×n dimensional molecularity matrix eq. (4) ψ:Rm→Rn monomial function of the kinetic dynamics, ψi(c)=∏mj=1cYjij eq. (1) ϕi net reaction rate corresponding to complex i eq. (6) Δ subspace containing im(Ak) eq. (10) Ξ stoichiometric subspace eq. (11) s dimension of the stoichiometric subspace subsec. 2.2 δ deficiency of the reaction network eq. (13)

## 1 Introduction

Deterministic reaction networks obeying mass action law (MAL) kinetics form an important subclass of kinetic systems, which in spite of their apparent simplicity, are able to describe a rich variety of dynamical behavior, that includes multiple equilibria conditions, oscillations or even chaos [19, 9]. Such networks are typically employed to describe the dynamics of open or closed chemical reaction systems, but over the last years they proved useful in modeling other system classes as well.

Reaction networks belong to the class of nonnegative (or positive) systems, the main characteristics of which being that the non-negative orthant is invariant for the dynamics. The application field of nonnegative systems extends far beyond chemistry, and includes dynamical models whose state variables are naturally nonnegative, as it is the case of biological systems in their many scales (from cells to ecological systems), or systems that can be transformed to be nonnegative, such as certain process models (e.g. heat exchangers, distillation columns, convection networks), economic, transportation or stochastic models [20]. With an appropriate selection of coordinates, even many classical mechanical and electrical models can be described in the nonnegative framework.

The main specialty of reaction networks within nonnegative polynomial models is the lack of so-called cross-effects, what defines an additional constraint between the monomial coefficients and exponents [19]. Still, the class of reaction network models is quite wide, and many non-chemical models can be brought into a kinetic form using simple transformations [9, 48]. Widely used examples of kinetic systems are compartmental systems [31] and Lotka-Volterra models into which most smooth nonlinear ODEs can be embedded [33]. These facts, clearly underline the importance of reaction network models and motivate us to attempt to look at general dynamical models through the glasses of kinetic systems.

The study of the relationships between chemical reaction structure and dynamic behaviour is the purpose of Chemical Reaction Network Theory (CRNT), a program formally proposed and developed in [4, 5, 40]. One of the earliest results on the relation between the solutions of nonlinear dynamical systems (including kinetic systems) and their associated directed graphs is published in [52]. Important cycle-related conditions on the stability of kinetic systems were given in [10]. An extensive stability analysis of reaction networks using algebraic and graph-theoretical tools can be found in [11]. Thermodynamically motivated Lyapunov-function-based stability analysis of kinetic systems, considering certain frequently applied model-simplification steps, is proposed in [32].

The seminal works in [35, 21] (collected in their most comprehensive form in [23]) explored the dynamic properties of MAL complex chemical systems, and contributed to equip CRNT with a mathematical formalism that has prevailed to present. It is important to remark here, that several different network structures may correspond to the same kinetic differential equations [35, 18]. Therefore, important network properties such as deficiency, weak reversibility or complex balancing, may vary among the possible reaction structures belonging to the same ODE model (see, e.g. [50, 36]).

One fundamental problem in CRNT is to decide from the structure or parameters of the network, whether it can exhibit or not multiple equilibria. An important early result in this field is the rigorous proof of the existence and uniqueness of thermodynamic equilibrium in a mixture of chemically reacting ideal gases [53]. The motivation in [49] was the computational analysis of large thermodynamical models. The work contains fundamental results about the existence and uniqueness of compositions minimizing the free energy.

In answering the questions about the properties of equilibria, the concept of network deficiency (a number that relates to reaction network structure and stoichiometry) has become central to characterize the network behavior. Two essential results of CRNT are the well-known deficiency zero and deficiency one theorems [23, 24] which (besides other important results) establish conditions for networks to have exactly one equilibrium point in each positive stoichiometric compatibility class [23]. This suggests network robustness with respect to parameter variability, and underlines the importance of the kinetic system class in general nonlinear systems theory.

CRNT has received renewed interest over the last years, particularly in the area of systems biology, because of its potential to explore and to analyse complex behavior and functionality in biological systems (e.g. [13, 42, 12, 45]). Most efforts were dedicated to investigate the relationships between reaction network structure and dynamic behaviour. In this regard, special mention should be made of the so-called injectivity property, investigated as a condition that relates to the singularity (or not) of the determinant of the Jacobian associated to a given dynamic system [16]. Algebraic and graph theoretical methods have been devised to check injectivity, and therefore uniqueness of equilibria [16, 17]. In the same direction, extensions to cope with instabilities have been developed in [42]. From different perspectives, a number of necessary and sufficient conditions for a given network structure and stoichiometry to accommodate multiple equilibria have been also recently proposed in [12, 45].

A particularly interesting class of chemical networks are the reversible ones, either in the strict thermodynamic sense, in which every elementary reaction step is reversible, or in a weak reversibility context. Reversibility leads to a particular set of positive equilibria which is known as detailed balance if each reaction step is equilibrated by a reverse one, or complex balanced if the network is weakly reversible.

At this point, it must be remarked that equilibrium should be understood along the sequel in the sense given in dynamic systems, irrespectively of whether it corresponds to thermodynamic equilibrium or to a particular steady-state on a chemical reactor. Note, however, that in agrement with thermodynamics, instabilities in the dynamics of reaction systems (when taking place on a homogeneous medium in isothermal conditions) require the reaction domain to be open to mass exchange with the environment.

Because of microreversibility, most chemical systems, when closed to mass and energy exchanges with the environment, satisfy the principle of detailed balance equilibrium, resulting into stable equilibria [29]. As discussed in [30] and [28], irreversibility can be allowed within a reaction network, as limit cases of reversible steps under a thermodynamic consistency condition (known as the Wegscheider condition) which necessarily assumes microreversibility.

The notion of complex balancing (also known as cyclic balancing or semi-detailed balancing), on the other hand, generalizes the detailed balance condition to any weakly reversible network. The structure of complex balanced systems has been explored in [15] and shown to be a toric variety with unique and stable equilibrium points (see also [51]). Extensions to cope with more general classes of kinetic systems have been investigated in [43, 46].

It is important to mention here the recent fundamental results on the proof of the Global Attractor Conjecture which says that any equilibrium point of a complex balanced mass action system is globally stable. A proof for the single linkage class case was given in [3], while a possible general proof based on differential inclusions was described in [14].

CRNT, as it stands nowadays within the field of applied mathematics, offers an extraordinary potential in system’s theory for analysis and design of complex dynamic systems of polynomial type, what in turn may cover a wide spectrum of chemical and biological systems. Unfortunately, many of its results remain at a large extent unexploited, when not unnoticed, in the fields of process systems and engineering.

Among the reasons that hamper application might be certain advanced mathematical tools and the intensive use of graph theory that are often not well-enough known to engineers. Some practical questions that demand attention relate to the link between dynamic behavior of a given mechanism and parameter sets (reaction rate coefficients), or to the design of a chemical/biochemical network with some pre-specified behaviour (e.g bistable, oscillatory, etc).

In this contribution we present some conditions that ensure feasibility of equilibrium solutions for weakly reversible mass action law (MAL) systems. They are linked to the notion of âfamily of solutionsâ, a concept originally derived in [44, 45] to study multiplicity phenomena as a function of network parameters.

In deriving what it will be referred in the sequel as feasibility conditions, we exploit a flux-based form of the model equation. Within such structure, the time evolution of the species concentration vector is expressed as the product of a matrix denoted by , whose columns span the stoichiometric subspace of the reaction system, and a vector function that is related to concentrations through a class of stable Metzler matrices [6].

As we will show, feasibility relates to the orthogonality between a log-transformed vector function of reaction complexes and the kernel of matrix . Based on this observation, feasibility conditions will be expressed in terms of certain functions that can be employed to identify admissible equilibria within the positive orthant of the concentration space. It will be shown that such functions are monotonous in their respective argument and take the zero within their domain, what will allow us to establish links with existence and uniqueness of equilibria along positive stoichiometric compatibility classes for MAL kinetic systems. In this context, connections between monotonicity and two classical results in CRNT theory that relate to complex balanced equilibrium [34, 35], and to a class of positive deficiency networks [22, 23, 24], will be discussed.

Finally, it must be remarked that the potential interest of the notion of complex balancing in the context of process control is to characterize stable operation regimes in open systems, where the principle of detailed balance does not necessarily hold. This may allow, for instance, the selection or manipulation of exchange fluxes so to preserve stability of the resulting (open to the environment) reaction system, via appropriate process optimization and/or feed-back control (see e.g. [41]). Future directions may also involve the detection or design of networks having multiple equilibria.

The paper is organized as follows: Section 2 introduces a formal description of chemical reaction networks. The graph structure underlying a reaction network, and its algebraic counterpart, will be described in Section 3. Section 4 presents a flux-based form canonical representation of the equilibrium set, that includes some feasibility conditions. Relationships between network structure and monotonicity of feasibility conditions will be established in Section 5. Connections between monotonicity of feasibility functions and some classical results on uniqueness and stability of equilibria will be discussed in Section 6.

## 2 Preliminaries: Reaction Network Structure and Dynamics

Let be the number of chemical species which react by irreversible chemical reaction steps, and the corresponding vector of species concentrations, defined as mole number per unit of volume. Each reaction step transforms some set of chemicals, usually referred to as reactants, into a set of reaction products. In CRNT, reactants and reaction products receive the name of reaction complexes. Complexes and reaction steps describe a graph where complexes correspond to nodes and reaction steps to directed edges.

Formally, a graph involving complexes linked by irreversible reaction steps can be constructed by associating to each complex a set with integer elements, and a vector . The elements of the set are the indices of the complexes that are directly reachable (i.e. by one reaction step) from . From now on, we will refer to each complex by the corresponding index . Vector has as entries the (positive) stoichiometric coefficients of the molecular species that participate in complex .

The graph structure is then built by linking every complex to . This process results in a number of connected components known in CRNT as linkage classes. For each linkage class , we define the set which contains as elements the indexes of the complexes that belong to that linkage class1.

Complexes are connected within a linkage class by sequences of irreversible reaction steps that define directed paths. Two complexes are strongly linked if they can be mutually reached from each other by directed paths (trivially, every complex is strongly linked to itself). A maximal set of pairwise strongly linked complexes defines a strong terminal linkage class if no other complex can be reached from its nodes. In this work we will consider only networks in which every linkage class contains just one strong terminal linkage class.

A linkage class is said to be weakly reversible if any pair of its complexes is strongly linked. Weakly reversible networks are those composed by weakly reversible linkage classes. A particular type of weakly reversible linkage class is a reversible linkage class if each reaction step is itself reversible, so that for every and , we have that . The rate , at which a set of reactants in complex is transformed into a set of products in complex , will be assumed to be mass action, so that:

 Rij(c)=kijψi(c),  with  ψi(c)=m∏j=1cyjij≡cyi, (1)

where is the stoichiometric vector corresponding to complex . The reaction systems we consider in this work will take place under isothermal conditions, what makes any reaction rate parameter constant. Whenever c is a strictly positive vector, the following alternative representation for may be more convenient:

 lnψi(c)=yTilnc, (2)

where the natural logarithm operator acts on any vector element-wise. Let be the vector containing as entries the monomials described in (1), then the previous expression can be written in matrix form as:

 lnψ(c)=YTlnc, (3)

where is the so-called molecularity matrix which collects as columns the stoichiometric vectors associated to the complexes of the network.

### 2.1 The dynamics of reaction networks

Following the classical work by Feinberg [21], the time evolution of species concentrations on a well-mixed reaction medium at constant temperature can be described by a set of ordinary differential equations that we write as:

 ˙c=Y⋅Ak(ψ(c))=Y⋅∑λAλk(ψ(c)), (4)

where is a linear operator defined as:

 Aλk(ψ)≡∑i∈Lλψi∑j∈Iikij⋅(εj−εi), (5)

with denoting the th standard unit vector employed to represent axes on a cartesian coordinate system. Let us define the net reaction rate flux around a complex , as the signed sum of in- and out-flowing fluxes, i.e. as a function of the form:

 ϕi(ψ)=∑{j|i∈Ij}Rji(ψ)−∑j∈IiRij(ψ), (6)

where the first summation at the right hand side extends to all source complexes in the network from which there exists a reaction step to product complex , and is represented by .

We can express in (5) in terms of fluxes (6), by selecting any reference complex , and adding and subtracting from the right hand side of (5) so that:

 Aλk(ψ)=∑i∈Lλψi∑j∈Iikij⋅(εj−εjλ)−∑i∈Lλ⎛⎝∑j∈Iikijψi⎞⎠⋅(εi−εjλ).

After switching subindexes, re-ordering the summations for the first term at the right hand side and making use of (6), we get the following equivalent expressions:

 Aλk(ψ) = ∑i∈Lλ⎛⎜⎝∑{j|i∈Ij}kjiψj⎞⎟⎠⋅(εi−εjλ)−∑i∈Lλ⎛⎝∑j∈Iikijψi⎞⎠⋅(εi−εjλ) (7) = ∑i∈Lλϕi(ψ)⋅(εi−εjλ).

For convenience, the reference complex will be chosen from the corresponding strong terminal linkage class. Since vectors are orthogonal, by using (7), we have that for every . Let , then we also have that and therefore:

 ∑i∈Lλϕi(ψ)=⎛⎝∑i∈Lλεi⎞⎠TAλk(ψ)=0. (8)

Note that fluxes in (6) (as well as the linear operator in (5)) are implicitly dependent on the reaction rate coefficients associated to the reaction steps in the linkage class. By inspection of (7), it can be concluded that the image of lies on the subspace defined as follows:

 Δ=Δ1+⋯+Δλ+⋯+Δℓ, (9)

where

 Δλ=span{εi−εjλ | i∈Lλ}% forλ=1,⋯,ℓ. (10)

and the sum of vector spaces and is defined as:

 V1+V2={v1+v2 | v1∈V1, v2∈V2}.

Since vectors in are linearly independent, they form a basis for the subspace , thus . In addition, since the subspaces are orthogonal:

 dim(Δ)=∑λ(Nλ−1)=n−ℓ.

This implies that if and only if for all . Consequently, if a positive concentration vector c exists compatible with a zero flux condition for every complex in the network, that vector should be an equilibrium for system (4). Such equilibrium condition, known as complex balanced [34], is formally defined as follows:

###### Definition 2.1

(Complex Balanced Equilibrium) Any vector such that (Eqn (6)) for every is called a complex balanced equilibrium solution.

A subclass of complex balanced equilibrium, particularly meaningful from a thermodynamic point of view as it relates to microreversibility ([37, 29]), is the detailed balance equilibrium which we define next:

###### Definition 2.2

(Detailed Balance Equilibrium) If the network is reversible (i.e. for every and , we have that ) any vector such that (where is of the form (1)) is called as a detailed balance equilibrium solution.

### 2.2 The stoichiometric subspace

Similarly to the subspace , we define the stoichiometric subspace as:

 Ξ=Ξ1+⋯+Ξλ+⋯+Ξℓ,

where:

 Ξλ=span{yi−yjλ | i∈Lλ}forλ=1,⋯,ℓ. (11)

In what follows it will be more convenient to collect the elements from each of the sets and their union, column-wise in matrices and , respectively, so that:

 S=[S1⋯Sλ⋯Sℓ]. (12)

Let , which eventually coincides with the rank of , then it follows from the rank-nullity theorem that the dimension of the kernel (null space) of will be:

 δ=n−ℓ−s. (13)

This number is known in CRNT as the deficiency of the network. In a similar way, we can define the deficiency of each linkage class as the dimension of the kernel of so that , where . Since , it is not difficult to conclude that linkage class and network deficiencies relate as:

 δ≥∑λδλ. (14)

Let be a basis for the kernel of , and express each vector in terms of sub-vectors (one per linkage class), so that:

 (gr)T=[(gr1)T⋯(grλ)T⋯(grℓ)T], for r=1,⋯,δ. (15)

Using the above description, equation can be re-written as:

 ∑λSλgrλ=0  for  r=1,…,δ. (16)

We will be particularly interested in solutions of system (4) on the convex region resulting from the intersection of the non-negative (respectively positive) orthant in the concentration space and a certain linear variety associated to the stoichiometric subspace , the result known in CRNT as a stoichiometric (respectively, positive stoichiometric) compatibility class. Given a reference concentration vector , the stoichiometric compatibility class can be formally defined as:

 Ω(c0)={c∈Rm  |  c≥0,PT(c−c0)=0}, (17)

where is a full rank matrix whose columns span the orthogonal complement . The corresponding positive stoichiometric compatibility class can expressed as . In passing, let us define the function as . Such function is constant along trajectories (4), since by combining (7) and (4) we have that:

 ˙c=∑λ∑i∈Lλϕi(ψ(c))(yi−yjλ), (18)

and the columns of are orthogonal to , hence . In other words, is an invariant of motion for system (4). From this observation it is not difficult to conclude that any trajectory that starts in a compatibility class will remain there.

### 2.3 Some examples of chemical reaction networks

A reversible chemical reaction network

Let us consider a reaction network involving molecular species we label as , each of them constituted by a combination of three types of functional groups (or atoms) we denote as , and . The (reversible) chemical reaction steps that take place are:

 A2B+C⇆AC+ABAB+2C⇆AC2BAC2B⇆AC+CB (19)

Molecular species and functional groups are related as follows: , , , , , . The network consists of complexes:

 {C1,C2,C3,C4,C5}≡{M1+M4,M2+M3,M3+2M4,M5,M2+M6}.

Making use of the formal description previously discussed, the sets that indicate which complexes are reached from complex become, for this example:

 I1={2}I2={1}I3={4}I4={3,5}I5={4}. (20)

The corresponding stoichiometric vectors associated to each complex are written as columns in the molecularity matrix :

 Y=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣100000100101100102000001000001⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (21)

The graph representation is depicted in Figure 1, and comprises two linkage classes and . The net reaction fluxes (6) around each complex are:

 ϕ1(ψ)=k21ψ2−k12ψ1ϕ2(ψ)=k12ψ1−k21ψ2ϕ3(ψ)=k43ψ4−k34ψ3ϕ4(ψ)=k34ψ3+k54ψ5−(k43+k45)ψ4ϕ5(ψ)=k45ψ4−k54ψ5 (22)

Note that if and are chosen as reference complexes, by relation (8), the fluxes associated to the reference become:

 ϕ1(ψ)=−ϕ2(ψ)andϕ3(ψ)=−(ϕ4(ψ)+ϕ5(ψ)).

The image of coincides with subspace (Eqn (9)) with and of the form:

 Δ1=span{(ε2−ε1)}Δ2=span{(ε4−ε3),(ε5−ε3)}. (23)

Matrices , employed in Section 2.2 to define (column-wise) the corresponding stoichiometric subspaces (11), are of the form:

 S1=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−111−100⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,  S2=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0001−1−1−2−21001⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

The dimension of the stoichiometric subspace, which coincides with the rank of matrix , is . Hence, network deficiency is what means that no vector other than the zero vector exists such that .

The numbers of atoms (or functional groups) - remain constant, provided that reactions take place on a closed domain (i.e. no mass exchanges with the environment occur). Let be total concentrations for , and on the closed and homogeneous domain. Mole-number balances result in the following set of linear relations:

 [A]0=2[A2B]+[AC]+[AB]+[AC2B][B]0=[A2B]+[AB]+[AC2B]+[CB][C]0=[AC]+[C]+2[AC2B]+[CB]

where brackets indicate chemical species concentrations. Previous relations can be written in matrix form as follows:

 PT(c−c0)=0  with  P=⎡⎢⎣211010101011010121⎤⎥⎦T,

where c is the vector of chemical species concentrations and is the (constant) concentration of functional groups/atoms. It must be noted that the above expression is employed in (17) to characterize the set of compatibility classes. As discussed in Section 2.2, because (full rank), the columns of define a basis for the orthogonal complement of the stoichiometric subspace.

An irreversible network

Let us consider the following set of irreversible reactions:

 S+I→2II⇆R→S (24)

This network can be interpreted as an extension of the SIR epidemic model [38] which describes the effect of a disease on a large population. Individuals on the population are classified either as those susceptible to the disease (), infected () or those recovered from the disease (). In this extension, individuals under recovering may evolve either to those susceptible to the disease or directly infected again. In the CRNT formalism, the network comprises three species, with concentrations , and , and complexes, numbered as:

 {C1,C2,C3,C4,C5}≡{2I,S+I,S,R,I}.

Graph structure is depicted in Figure 2. Molecularity matrix for this network reads:

 Y=⎡⎢⎣011002100100010⎤⎥⎦. (25)

Choosing as reference complexes and , the matrices become:

 S1=⎡⎢⎣1−10⎤⎥⎦,  S2=⎡⎢⎣−1−10110⎤⎥⎦.

Because both matrices are full rank, the dimension of and is and , respectively. Thus and . However, matrix is rank deficient since vector is parallel to the vector in the second column of . Consequently , and , verifying inequality (14). A basis for the kernel of is given by the vector .

For this network, the basis that spans the orthogonal complement of the stoichiometric subspace corresponds to . Each compatibility class is given by (see (17)), the region of non-negative concentrations () satisfying:

 [S]+[I]+[R]=PTc0,

with being a constant vector. A representation of a compatibility class for the reaction network considered is presented in Figure 3, with .

## 3 Linkage classes, Graphs and Compartmental Matrices

The nature of equilibrium solutions in chemical reaction networks (e.g. positivity, instability of equilibrium points, or the possibility of multiple equilibria) is at a large extent determined by the graph structure of each linkage class and the properties of some matrices associated to it, that belong to the class of compartmental matrices [26].

In this section, we describe such matrices and discuss their properties, with emphasis on invertibility and non-negativity of their inverses. The main results, summarized in Lemma 3.1, will be extensively employed in the sequel. For the sake of completeness, we introduce a derivation from scratch, while establishing connections with known facts in the field of positive linear systems and non-negative matrices.

To that purpose, we introduce a graph related to the graph description of a linkage class, and a matrix associated to it. The properties of this matrix will be studied by constructing an auxiliary linear dynamic system and examining the corresponding equilibrium.

A directed graph is constructed by a set of vertices containing a distinguished vertex and a set of edges. The first set of vertices , with indexes will be referred to as nodes, while the remaining vertex will represent the ‘environment’. To any directed edge for and in (i.e. ) there corresponds a scalar weight . In addition, for every , such that or , we associate scalar weights and , respectively. Such weights will be collected as entries in (non-negative) vectors , with (respectively, ) if there is no directed edge (respectively, ). As in the description of linkage classes (Section 2) we say that two nodes are strongly linked if they can be reached from each other by directed paths. The maximal set of strongly linked nodes (not passing through the environment) in from which there are no outgoing edges to other nodes, defines a strong terminal set, we will refer to as .

Since in this work we are interested in linkage classes with just one strong terminal linkage class, the graphs we consider will only contain one strong terminal set. The set of non-terminal nodes is defined as . Some examples of directed graphs are illustrated in Figure 4.

We associate to the graph a state vector , where each component corresponds to a node. For each node , we define an internal net flux and a net exchange with the environment . Each flux is of the form:

 ϕi(z)=∑{j|i∈Ij}Vjizj−zi∑j∈IiVij, (26)

where as in Section 2, the indexes of the nodes that are directly reached from node are grouped in a set , and refers to all nodes with edges directed to . In addition, the net exchange with the environment is expressed as . From (26), and similarly to expression (8), it is straightforward to see that:

 ∑i∈Lϕi(z)=0. (27)

We consider that for each node , the state will evolve in time as a function of the corresponding internal and external net fluxes, so that . Combining this expression with (27), the dynamics of the system can be described as:

 ˙z=n∑i=1ϕi(z)(εi−ε1)−Bz+a, (28)

where , a diagonal matrix with the components of b in the diagonal. Equation (28) can be re-written in the alternative form:

 ˙z=Wz+a (29)

where matrix , and is the matrix that contains as off-diagonal components the coefficients in (26), with if there is no directed edge . The diagonal elements for are of the form:

 Vii=−n∑j=1j≠iVij.

Both matrices and are compartmental [26], 2 and belong to the class of Metzler matrices [6]. It is known from e.g. [7], that the eigenvalues of a compartmental matrix are either zero or they have negative real parts. Such conclusion can be also reached from the structure of compartmental matrices by applying the Gershgorin disc theorem [27].

###### Definition 3.1

We say that the matrix associated to the system (29) is C-Metzler if its entries are of the form:

 Wij≥0for  i≠jWii=−(bi+∑j≠iWji)with  bi≥0  and  Wii<0, (30)

with at least one positive associated to the strong terminal set (i.e. ).

###### Proposition 3.1

Consider system (29) with and nonnegative initial conditions . Then for every . 3

Proof: In order to prove the statement all we need is to show that the flow associated to the differential system on the boundary of the positive orthant is either aligned to the boundary or oriented to the interior of the orthant. Before we compute the flow, let us define the set which characterizes the -th facet of the positive orthant. The inner product between the flow induced by (28) (equivalently (29)) on any element , and the unit vector orthogonal to takes the form:

 εTk˙z=ϕk(z)−bkzk+ak=∑{j|k∈Ij}Vjkzj+ak≥0 (31)

where the equivalence at the right hand side holds since in . Thus, at the boundary of , the flow associated to the differential system will be either aligned to the boundary () or oriented to the interior of the positive orthant (i.e. ). Repeating the argument for all values of completes the proof.

Next, we will study the equilibrium of system (29) that results from some constant non-negative vectors a and b. To that purpose, let be the set containing those non-terminal nodes that are directly linked to any node in the strong terminal set (the different partitions are illustrated in Figure 5). Let us introduce functions , and , where:

 ωp=∑i∈Lpεi,  ωq=∑i∈Lqεi. (32)

The time derivative of along system (28) is of the form:

 ˙σq = ωTqa−∑i∈Lqbizi−∑i∈Lℓ ∑j∈IiVijzi (33) ˙σp = ωTpa−∑i∈Lpbizi+∑i∈Lℓ ∑j∈IiVijzi (34)
###### Proposition 3.2

Consider system (29) with and C-Metzler (Definition 3.1). Then, for any , is the only equilibrium solution of (29).

Proof: Since , by Proposition 3.1 we have that will remain in the positive orthant, and will be non-negative, and .

For every node it is clear that is the only possible equilibrium. Otherwise, from (33), it follows that for all times, what would make to become negative. In addition, note that for to be an equilibrium point for node , requires the equilibrium states for all nodes linked to (thus ) to be zero as well, otherwise . Repeating the argument upstream to the nodes linked to , we conclude that zero is the only possible equilibrium for every node in .

Since is C-Metzler (i.e. there exists at least one index for which ), zero is also the only possible equilibrium for every node in the terminal set. Suppose, on the contrary, that there exists a positive equilibrium for some . In the limit, expression (34) would become:

 limt→∞˙σp=−∑i∈Lpbiz∗i<0,

but cannot become negative, thus . Finally, since is an equilibrium point for such that , the equilibrium states for all nodes linked to must be zero as well, otherwise . repeating the argument for every we have that .

###### Remark 3.1

Using the terminology of [26], defines a compartmental system with one trap, where the notion of trap is equivalent to the definition of strong terminal set used in this paper. Therefore, zero is an eigenvalue of with multiplicity . is also a compartmental matrix, and due to the choice of , the compartmental system corresponding to contains no traps. Therefore, is of full rank and thus, assuming , the only equilibrium of (29) is .

###### Proposition 3.3

Consider system (29), with matrix being C-Metzler (Definition 3.1). Equilibrium will be non-negative and globally asymptotically stable. Moreover, let some entry of vector positive (i.e. ). Then, for all nodes reached from node by directed paths, we have that .

Proof: First note that because matrix in (29) is compartmental, its eigenvalues are either zero or they have negative real parts. Thus, to show that the equilibrium is globally asymptotically stable it only remains to prove that no zero eigenvalue exists. Actually, this is the case, since from Proposition 3.2, the only equilibrium solution for which is . Since all eigenvalues have negative real part, is invertible and the equilibrium is globally asymptotically stable.

In proving the second part of the statement, we have that since , the equilibrium in node must be:

 z∗k=∑{j|k∈Ij}Vjkz∗j+akbk+∑j∈IkVkj>0

and so is the case for all reached from , so that , what completes the proof.

###### Remark 3.2

The full rank property of and the uniqueness of the equilibrium for are equivalent, since this latter means that the dimension of the kernel of is zero. From this, it follows that cannot have zero eigenvalues, and from the compartmental property of we obtain that all the real parts of its eigenvalues are negative. This means that is a stability matrix. Since a can be considered as a bounded input, is a unique asymptotically stable equilibrium point of (29). Since is a full-rank M-matrix, it is inverse-nonnegative [47], i.e. all entries of are non-positive. This implies that is a nonnegative vector.

###### Lemma 3.1

Any C-Metzler matrix is non-singular and its inverse non-positive. Let its associated graph be numbered so that the first p nodes are in (thus, the remaining are in ). Then can be partitioned as:

 N=[NpNpq∅Nq], (35)

with the zero matrix, and strictly negative matrices, and non-positive. Moreover, for each node , every node not reached from by directed paths will correspond with an entry .

Proof: That (being C-Metzler) is non-singular and therefore invertible has been shown in the proof of Proposition 3.3 (see also Remark 3.2). In addition, note that the equilibrium solution can be written as , and choose . For each , the corresponding equilibrium will be , where represents the -th column of matrix .

According to Proposition 3.3, the equilibrium state for all nodes reached from by directed paths must be positive, and thus the corresponding entries must be negative. This holds for all . Thus, using the order given in the statement of the Lemma, we can conclude that is strictly negative. In addition, every can be reached from what makes strictly negative as well. The zero matrix appears since none of the nodes can be reached from nodes . Finally, zero entries in