First-order chemical reaction networks I: theoretical considerations

O R I G I NA L PA P E R

First-order chemical reaction networks I: theoretical considerations

Roland Tóbiás¹ · László L. Stacho² · Gyula Tasi¹

Received: 21 February 2016 / Accepted: 4 June 2016

© Springer International Publishing Switzerland 2016

Abstract Our former study Tóbiás and Tasi (J Math Chem 54:85,2016) is continued, where a simple algebraic solution was given to the kinetic problem of triangle, quad- rangle and pentangle reactions. In the present work, after defining chemical reaction networks and their connectedness, first-order chemical reaction networks (FCRNs) are studied on the basis of the results achieved by Chellaboina et al. (Control Syst 29:60, 2009). First, it is proved that anFCRNis disconnected iff its coefficient matrix is block diagonalizable. Furthermore, mass incompatibility is used to interpret the reducibility of subconservative networks. For conservativeFCRNs, the so-called marker network is introduced, which is linearly conjugate to the original one, to describe the zero eigenvalue associated to the coefficient matrix of anFCRN. Instead of using graph- theoretical concepts, simple algebraic tools are applied to present and solve these problems. As an illustration, an industrially important ten-component (formal)FCRN is presented which has algebraically exact solution.

Keywords First-order reaction network·Algebraic model·Network decomposition· Mass incompatibility·Marker network·Multiplicity of the zero eigenvalue

1 Introduction

The graph-theoretical formalism of the chemical reaction networks (CRNs) was elab- orated in the 1970s by Horn, Jackson and Feinberg [16,18,25] which has received

B

tasi@chem.u-szeged.hu

1 Department of Applied and Environmental Chemistry, University of Szeged, Rerrich B. tér 1, Szeged 6720, Hungary

2 Bolyai Institute, University of Szeged, Aradi Vértanúk tere 1, Szeged 6720, Hungary

permanent applications [11,12,17,41]. Within the framework of this theory, numer- ous new kinetic concepts (e.g. linkage class, stoichiometric compatibility class, weak reversibility, complex and detailed balanced network) were introduced and results of high interest (multistability [11,12], deficiency [17], global attractor [10] and persis- tency theorems [2,36], conjugacy ofCRNs [13,26,27]) were achieved.

CRNscan be defined not only in the language of graph theory. In the study pub- lished by Chellaboina et al. [7], matrix-vector notation is applied and the dynamical equations are set up in terms of vector-matrix exponentiation. This approach resulted in significant progress in the characterization of the solutions connected to the mass- balance relations. As far as the algebraic treatment of the first-order linear systems of differential equations (FLSODE) is concerned, the books of Pontryagin [38] and Kailath [28] should be referred, where relevant information can be found about their structural properties, especially about their stability.

The qualitative theory of first-order reaction networks (FCRNs) was limited to the problems (stability [29], physical realizability [9], observability, controllability, identifiability [8], decomposability [34]) of compartmental systems [19,22,23], in chemical terms: isomerization reaction networks,IRNs. During the last two decades, the properties of arbitraryFCRNs have also been investigated in detail by Bernstein et al. [4,5,7].

In this study, first we summarize the significant results of Chellaboina et al. [7]

related toFCRNs and then we answer some further questions (e.g. mass incompati- bility, multiplicity of the zero eigenvalue of the coefficient matrix). Instead of using graph-theoretical concepts, simple algebraic tools are applied to present and solve these problems. As an example, an industrially important ten-component (formal) FCRNis presented which has an algebraically exact solution.

2 Preliminaries

2.1 Algebraic model ofCRNs

Consider a homogeneous reaction system with K chemical components and R ele- mentary reactions at constant temperature, pressure and volume according to the next scheme [7,43]:

K j=1

di jAj → K

j=1

gi jAj

i ∈Z⁺_R

(1)

wheredi j andgi j are the left and right stoichiometric coefficients of the species Aj

in the ith reaction. Equation (1) can be written with the matrices D = di j

and G=

gi j

as follows:

DA→GA (2)

whereA= Aj

denotes the vector of the chemical components. In particular, regard- ing elementary reactions,di jandgi j are nonnegative integers. If each of the reactions is of first-order, thenD=

δ_ηij

whereδ_ηijstands for the Kronecker delta and there

exists a uniqueηi ∈Z⁺_K for everyi ∈Z⁺_R. It is convenient to introduce the so-called stoichiometric matrixS=

νi j

:S=G−D.

The time dependence of the concentration vectorC= cj

is described with the following system of ordinary differential equations (mass-balance relations):

C˙ =S^Tρ (3)

whereC˙ =

˙ cj

is the derivative ofCwith respect to time,ρ= {ri}is the vector of the reaction rates,S^T is the transpose ofS. Reaction rates are approached via thelaw of kinetic mass-action[43]:

ri =ki

K l=1

c^d_l^il (4)

whereki >0 is theith elementary rate coefficient which depends on the temperature and, occasionally, on the pressure as well.

Assuming mass-action kinetics, the quadruple D,G,A,k is called the chem- ical reaction network (CRN) if k = {ki}. For arbitrary permutation matrices PR ∈ {0,1}^R×R and PK ∈ {0,1}^K×K,D,G,A,k is equivalent to the network PRDPK,PRGPK,P^T_KA,PRk.

To characterize the behavior ofCRNs concerning mass conservation, three impor- tant concepts [14] are necessary to mention. A networkD,G,A,kis

– conservativeif there existsM∈(0,∞)^Ksuch thatSM=0Rwhere0R = {0}^R; – subconservativeif there is anM∈(0,∞)^KthatSM∈(−∞,0]^R;

– superconservativeifSM∈[0,∞)^Rwith an appropriateM∈(0,∞)^K. 2.2 Decomposition of CRNs

Consider the partitions I =

Ik:k∈Z⁺_NC

and J =

Jl :l∈Z⁺_NC of Z⁺_R andZ⁺_K where NC is the number of the partition cells. The couple(I,J)decom- poses D,G,A,k if di j = gi j = 0 for i ∈ Ik and j ∈ Jl with k = l. Suppose that the couples (I,J) and

I,J

decompose D,G,A,k where I =

I_l:l ∈Z⁺_N C

and J =

J_l:l ∈Z⁺_N C

. Then I,J

with I = Ik∩I_l:k∈Z⁺_NC∧l∈Z⁺_N

C ∧Ik∩I_l=∅

andJ=

Jk∩J_l:k∈Z⁺_NC∧l ∈ Z⁺_N

C∧Jk∩J_l=∅

does the same. Clearly, there is a finest decomposing couple of partitions(I^∗,J^∗)with a maximal number of cellsN_C^∗ whereI^∗=

I_l^∗:l ∈Z⁺_N∗ C

andJ^∗=

J_l^∗:l∈Z⁺_N∗ C

. Notice, there are permutation matricesPRandPK block diagonalizingDandG:

PRDPK =diag_i∈Z+ N∗C

(Di) (5)

PRGPK =diag_i∈Z+

∗(Gi) (6)

whereDi,Gi ∈ R^Ri×Ki,diag_i∈Z+ N∗C

(Di) =diag

D1,D2, . . . ,DN_C^∗

,Ri and Ki is the cardinality of I_i^∗ andJ_i^∗(i ∈ ZN_C^∗). The matrices PR andPK partitionate the vectorsAandk, too:

P^T_RA=

A^T₁ A₂^T . . . A^T_N∗ C

T

(7) PRk=

k₁^T k^T₂ . . . k^T_N∗ C

T

(8) wherekiandAiare of the typeRi×1 andKi×1. The networkDi,Gi,Ai,kiis called amaximally independent subnetworkofD,G,A,k, and the set

Di,Gi,Ai,ki : i ∈ Z⁺_N∗

C

is a network decompositionof D,G,A,k. A network isconnectedif N_C^∗ =1, otherwise, it isdisconnected.

The matricesDandGcan be written in the form of Eqs. (5) and (6) iff the matrix H=G+Dis transformed by appropriate matricesPRandPK as follows:

PRHPK =diag_i∈Z+ N∗C

(Hi) . (9)

whereHi =Gi+Di. IfN_C^∗ =1, thenHis termed block diagonalizable. To determine PRandPKrequired by Eq. (9), the study of Schuster and Schuster [40] is worth citing.

2.3 Dynamics ofFCRNs

Henceforward, first-order reaction networks (FCRNs) are examined exclusively. The networkD,G,A,kisfirst-orderifD=

δηij

. AnFCRNis referred to asan iso- merization reaction network(IRN) ifG=

δκij

whereκi ∈Z⁺_K(i∈Z⁺_R). Regarding FCRNs, we get

C˙ =FC (10)

whereF = fj m

is the time-independent Jacobian ofC˙ with respect to the vector C.F is called the coefficient matrix of FCRNwhich can also be expressed in the subsequent way [7]:

F=S^Tdiag(k)D (11) where diag(k)=diag(k1,k2, . . . ,kR). The entry fi j can be written as follows:

fi j = R l=1

R m=1

νliδlmklδηmj = R m=1

νmiδηmjkm. (12)

Sinceνmi =gmi−δ_ηmi, the diagonal and off-diagonal entries need to be calculated separately:

fi j =

⎧⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎩ R m=1 ηm=i

(gmi−1)km, ifi = j R

m=1 ηm=j

gmikm, otherwise.

(13)

Equation (13) clearly shows that the relation fi j ≥0 implies ifi = j, thereforeFis a Metzlerian [7].

Equation (10) is a FLSODEwith initial conditionC0 = C(0) = c0j

whose solution is

C(t)= ∞ k=0

t^k k!

d^(k)C(t) dt^k

t=0

= ∞ k=0

t^k

k!F^kC0=exp(Ft)C0 (14) where the matrix exp(Ft)is used [30,31]. After linear algebraic transformations [32, 43], we get

C(t)=FC0V⁻¹E(t) (15) where

FC0 =

C0,FC0,F²C0, . . . ,F^K−1C0

(16) E(t)=

e^λ1t,t e^λ1t,t²e^λ1t, . . . ,t^μ1−1e^λ1t, . . . ,e^λLt,t e^λLt,t²e^λLt, . . . ,t^μL−1e^λLt T

(17) V=

E(0),E˙(0),E¨(0) , . . . ,E^(K−1)(0)

. (18)

FC0,E(t)andVin Eqs. (16)–(18) are the Krylov matrix, thetime evolution vector and the Vandermonde matrix, respectively. The parameterLis the number of distinct eigenvalues ofF,λkis thekth eigenvalue,μkis its multiplicity(k∈Z⁺_L). The columns of Vcorrespond to the various derivatives of the vector E(t) taken at t = 0. If the eigenvalues of the matrix F are calculated numerically, Eq. (15) is called the semianalytical solutionof Eq. (10).

Notice, that several algorithms (e.g. classical integration [24], transfer function [1,33], matrix and convolution [37] methods) can be used for solving Eq. (10).

2.4 Nonnegativity and semistable nonnegative equilibria

Since the entries ofC(t)represent the concentrations belonging to the species of an FCRN, thereforeC(t) ∈ [0,∞)^K (nonnegativity) andC_∞ ∈ [0,∞)^K (semistable nonnegative equilibrium) are indubitable for allt ≥ 0 whereC_∞ = limt→∞C(t).

Furthermore, based on [7], the nonnegativity and stability of Eq. (14) are investigated.

(The analysis of nonnegativity concerning to arbitraryCRNs can be also found in the work by Volpert [44].)

Nonnegativity of the functionC(t) =exp(Ft)C0implies from the well-known theorem that exp(Ft)∈ [0,∞)^K×K is valid for allt ≥ 0 iffFis a Metzler matrix [20]. As this is fulfilled by Eq. (13), so exp(Ft)C0∈[0,∞)^KifC0∈[0,∞)^K.

To examine the stability behavior of Eq. (14), some notions need to be recalled.

The vectorC˜ is an equilibrium point ifFC˜ =0K.C˜ is Lyapunov stable if there exists aθ > 0 for allt ≥ 0 andυ >0 such thatC− ˜C < υifC˜ −C0 < θ where denotes the Eucledian norm. A Lyapunov-stabile pointC˜ is regarded semistable if there isυ >0 such thatC˜ −C0< υinvolves Lyapunov stability of the vectorC_∞ for everyC0. It is proved [3] thatC˜is semistable iff each eigenvalue ofFhas a negative real part or is zero and the number of linearly independent eigenvectors related to the zero eigenvalue is equal to its (algebraic) multiplicity. A sufficient condition can also be derived [20]: ifF^TM∈(−∞,0]^K forM∈(0,∞)^K, thenC˜ =0K is a semistable equilibrium point of Eq. (10). This constraint is fulfilled by subconservative systems because

F^TM=

S^Tdiag(k)D T

M=D^Tdiag(k)SM (19) where SM ∈ (−∞,0]^R, consequentlyF^TM ∈ (−∞,0]^K. Notice, that the vector Mcan be replaced by the vector of molar masses; therefore, if these quantities are available, there is no need to search for such a vector.

3 Results

3.1 Relation of block diagonalizability of F and H

In this section it is revealed that the following two predicates are equivalent forFCRNs:

A)His block diagonalizable,B)Fis block diagonalizable.

Proof ofA)⇒B)Suppose thatPRandPKtransformHinto block diagonal form in terms of Eq. (9). Then bothDandGare block diagonal and diag(k)is partitioned into blocks by the permutation matrixPR:

PRdiag(k)P^T_R=diag_i∈Z+ N∗C

(diag(ki)) . (20)

Adapting these considerations, F=diag_i∈Z+

N∗ C

(Si)^Tdiag_i∈Z+ N∗

C

(diag(ki))diag_i∈Z+ N∗

C

(Di) (21)

where Si = Gi −Di. In this expression, the matrixi = S_i^Tdiag(ki)Di can be utilized which yieldsF=diag_i∈Z+

N∗ C

(i), thereforeA)⇒B).

Proof ofB)⇒A). From Eq. (13), it follows that fi j vanishes only if each term of the sum is identically zero as a result of gmi ≥ 0 andkm > 0. This means that gmi =0 comes from fi j = 0 andηm = j, i.e. Aj is not converted into Ai; fi j = fj i = 0 implicates that Ai and Aj do not take part in the same reaction. LetF = diag_i∈Z+

N∗ C

(i) (i ∈ R^Ki×Ki;i ∈ Z⁺_N∗

C)be considered. In this respect, Acan be partitioned by a permutation matrixPK =

δi j

in the form of Eq. (7). SinceFis a block diagonal matrix, the species belonging toAi andAj do not participate in the same conversion(i = j;i,j ∈ Z⁺_N∗

C). There needs to be such a permutation matrix PRwhich block diagonalizesD:

PRD=diag_i∈Z+ N∗

C

(Di) (22)

where Di ∈ {0,1}^Ri×Ki and Ri is the number of reactions taking place among the species ofAi(i ∈Z⁺_N∗

C). This permutationPRhas an impact on the shape of the matrix G:

PRG=

1^T 2^T . . . ^T_N∗ C

T

(23) wherei ∈[0,∞)^Ri×K. Leti be partitioned into the subsequent form:

i =

Gi1 Gi2 . . . G_{i N}∗ C

(24) whereGi j = [0,∞)^Ri×Kj(i,j ∈ Z⁺_N∗

C). With the help of Eqs. (22) and (24),PRG can be written as follows:

PRG=

⎛

⎜⎜

⎜⎝

G11 G12 . . . G1N_C^∗

G21 G22 . . . G2N_C^∗

... ... ... ...

GN_C^∗1GN_C^∗2. . .GN_C^∗N_C^∗

⎞

⎟⎟

⎟⎠. (25)

However,Gi j =0Ki×Kj(0Ki×Kj = {0}^Ki×Kj;i = j)is caused by the block diag- onality ofF, in one word,PRG = diag_i∈Z+

N∗ C

(Gii). In this instance,PRHis block diagonal, i.e.B)⇒A). Our further statements are drawn up for connectedFCRNs, namely, the matrixFis not deemed to be block diagonalizable.

3.2 Mass incompatibility in subconservative systems

Recall from Sect. 2.4 that subconservative systems exhibit semistable behavior.

Supposing a nonnegative integer-type matrix G, subconservativity implies certain restrictions on the kinetics of the networkD,G,A,k. It is enough to mention, the species Ai cannot be converted into Aj if Mi < Mj (mass incompatibility). In the sequel, we investigate the consequences of the conditionSM ∈(−∞,0]^R. Themth entry of the vectorSMis written as follows:

SMm = GMm− DMm = K

j=1

gm jMj−M_ηm. (26)

On the basis of the relation betweenMj andM_ηm,GMmis divided into three terms:

SMm =

j∈Z⁺K

M_j<M_ηm

gm jMj +M_ηm

j∈Z⁺K

M_j=M_ηm

gm j+

j∈Z⁺K

M_j>M_ηm

gm jMj−M_ηm. (27)

SinceSMm ≤0, the entrySMm can be estimated from below by neglecting the terms related toMj <M_ηm andMj =M_ηm:

0≥ SMm≥

j∈Z⁺K

M_j>M_ηm

gm jMj −M_ηm ≥M_m^⊗

j∈Z⁺K

M_j>M_ηm

gm j−M_ηm (28)

whereM_m^⊗= min

j∈ZK+

Mj>M_ηm

Mj. Rearranging Eq. (28), we obtain

j∈Z⁺K

Mj>M_ηm

gm j ≤ M_ηm Mm^⊗

<1 (29)

whereM_ηm/M_m^⊗<1. Since the entries ofGare nonnegative integers, we havegm j =0 forMj >M_ηm, establishing mass incompatibility. Based on all these, we get

SMm =

j∈Z⁺K

Mj<M_ηm

gm jMj+M_ηm

j∈Z⁺K

Mj=M_ηm

gm j−M_ηm. (30)

EstimatingSμm from below and simplifying byM_ηm, we conclude

1≥

j∈Z⁺K

M_j=M_ηm

gm j. (31)

This can be interpreted as at most one species Aj is produced in each reaction with the massMj =M_ηm, involvinggm j =1. Assuming thatMl =M_ηm andgml=1 for somel∈Z⁺_K, due to Eq. (30), we acquire

0≥

j∈Z⁺K

M_j<M_ηm

gm jMj +M_ηm −M_ηm =

j∈Z⁺K

M_j<M_ηm

gm jMj. (32)

In other words,gm j =0 for every j ∈Z⁺_K in case ofMj >M_ηm.

Mass incompatibility affects not onlyGbut alsoF. It is easily conceded that fi j =0 ifMi >Mj(i,j ∈Z⁺_K): sinceAjis not transformed intoAi, hence everygm jis equal to zero in Eq. (13), i.e. fi j =0. Consequently, permuting the chemical components in accordance with decreasing order of their masses, a lower block triangular matrix is received:

F=

⎛

⎜⎜

⎜⎝ F11

F21 F22

... ... ...

FN_W FN_W2. . .FN_WN_W

⎞

⎟⎟

⎟⎠ (33)

where NW is the number of species with different masses and the blocks explicitly not marked are zeros. This form is really advantageous because the eigenvalues ofF are identical to the ones of the matricesFii with significantly smaller size(i∈Z⁺_NW).

The role of block triangularity is expounded in more detail in Sect.3.4.

3.3 Concept of the marker network

Most chemical processes are studied in a closed system therefore paying particular attention to conservativeFCRNs is advisable. In this section it is explored how to trans- form a conservative network into a (formal) network containing only isomerization reaction steps.

Equation (19) clearly shows that F^TM = 0K is a consequence of the criterion SM=0R. Introducing the “all-onces vector”1K ∈ {1}^K, we can write

F^TM=F^Tdiag(M)1K =0K. (34) Equation (34) is invariant under the multiplication with an arbitrary matrix, therefore

diag⁻¹(M)F^Tdiag(M)1K =

diag(M)Fdiag⁻¹(M)T

1K =F^T1K =0K

(35) where diag⁻¹(M)=diag

M₁⁻¹,M₂⁻¹, . . . ,M_K⁻¹

andF=diag(M)Fdiag⁻¹(M).

It follows that

f_{i j} = K k=1

K l=1

δi kMifklδl jM_l⁻¹= Mi

Mj

fi j

i∈Z⁺_K

. (36)

Since fi j ≥0, therefore f_{i j} ≥ 0(i = j), and the diagonal entries f_ii can be deter- mined in the following way [see Eq. (35)]:

f_ii = − K

k=1 k=i

f_ki. (37)

Equation (35) suggests, there is a (formal)IRN δ_ηlj

, δ_κlj

, A_j

, k_l

with the coefficient matrixFwhere the mass of the species A_j is equal to one(η_l, κ_l ∈ Z⁺_K;l ∈Z⁺_R). Let us attempt to define such a network. At first, the following constraints are considered:η_k=η_l ⇒κ_k =κ_l;κ_k =κ_l⇒η_k =η_l;η_k =κ_k(k=l;k,l ∈Z⁺_R) excluding „self reactions” and „parallel reactions” from A_η

l into A_κ

l. Let Rbe the number of nonzero off-diagonal entries in the matrixF, andk_l = f_κ

lηl = 0. This network contains A_i → A_jwith the rate coefficient f_{j i} iff f_{j i} =0(i = j).

In this context,

D,G,A,k is referred to as themarker networkofD,G,A,k ifD =

δ_ηlj

,G = δ_κlj

,A = A_j

andk = k_l

. Next we prove thatF = G−DT

diag k

D. Denoting the matrix

G−DT

diag k

DbyF˜ = f˜i j

, off-diagonal entries f˜i j can be given

f˜i j =

R

m=1 ηm=j

δ_κmif_κ mηm =

R

m=1 ηm=j

δ_κmif_κ

mj, (38)

in agreement with Eq. (13). Obviously f_κ

mj = 0 is fulfilled just in case ofκm =i.

Only a unique indexm∈Z⁺_R has this property, thus f˜i j = f_{i j}. Using Eq. (38):

f˜ii = R m=1 ηm=i

δ_κmi−1 f_κ

mηm = −

R

m=1 ηm=i

f_κ

mi (39)

where δ_κmi = 0 is applied at the hand of η_m = κ_m . Regarding the caseη_m = i, summation needs to be performed on the row indices of nonzero off-diagonal entries related to theith column ofF:

f˜ii = − K

k=1 k=i; ˜f_ki=0

f_ki = − K

k=1 k=i

f_ki = f_ii. (40)

It can be shown thatC =diag(M)CwhereCis the concentration vector associa- ted to the vectorA. For this purpose, we need to construct the subsequent system of equations:

C˙=FC. (41)

Multiplying Eq. (10) by diag(M), Eq. (41) is obtained, since diag(M)C˙ =diag(M)F

diag⁻¹(M)diag(M)

C=Fdiag(M)C. (42) Moreover, the relationC =diag(M)Cestablishes a linear conjugacy relationship [26] between the vectorsCandC, i.e.

D,G,A,k is a linearly conjugate network ofD,G,A,k.

Based on all these, it is evident that the network

D,G,A,k indicates the dynamic behavior ofD,G,A,k, thus in certain cases (for instance, in character- izing the eigenvalues of the matrixF) examining the properties of the marker network may be sufficient.

3.4 Multiplicity of the zero eigenvalue in conservativeFCRNs

As presented for subconservativeFCRNs in Sect.3.2, if the molar masses of the species are not equal, an appropriate permutation matrixPK transforms the matrixFinto a lower block triangular form. In conservative systems, structure of Eq. (33) can be refined, furthermore, it can be derived under what conditions a diagonal block has a zero eigenvalue.

Knowing thatFis a Metzlerian with the effect ofF^T1K =0K, it is corresponded to the Laplacian of a simple weighted digraph, accordingly, the singularity of the diagonal blocks can be handled on the basis of graph theoretical principles [35]. In this section, recalling the works by Hearon and Taussky [22,23,42], we follow a path along which the problem of the zero eigenvalue can be studied purely with algebraic tools. Before, some linear algebraic notions need to be introduced.

The matrixQ= qi j

(Q∈C^K×K)isirreducible[6] ifK =1 or for each partition U,U

ofZ⁺_K(U,U = ∅)there existu ∈U andu ∈ Usuch thatquu = 0.Qis reducibleif it is not irreducible. Notice, that if a matrixQ∈[0,∞)^K×Kis irreducible, then max_i^K₌₁|ωi|is a simple eigenvalue whereωi denotes theith eigenvalue ofQ(i ∈ Z⁺_K)(Perron–Frobenius theorem [6]). In case ofK >1,Qis reducible iff it can be written in the next form:

PKQP^T_K =

⎛

⎜⎜

⎜⎝ Q11

Q21 Q22

... ... ...

Qn_T1Qn_T2. . .Qn_Tn_T

⎞

⎟⎟

⎟⎠ (43)

wherePK is a suitable permutation matrix andnT >1. IfnT ≥1 and the blocksQii

are irreducible, the matrixPKQP^T_K is called theFrobenius normal form[21] of the matrixQwhich can be determined by several methods (e.g. Tarjan algorithm [39]).

Later we also use a theorem proved by Taussky [42]: ifQ∈R^K×Kis an irreducible Metzlerian and−qii ≥K

j=1,j=iqj i(i ∈Z⁺_K), then det(Q)=0 is equivalent to the equation Q^T1K = 0K. (Although Taussky stated this theorem for the matrix−Q, the change of sign above is trivial, which is justified by the further considerations.)

Based on the Pearron–Frobenius theorem, it is also shown by Hearon [23] that the zero eigenvalue is simple for every irreducible MetzlerianQ∈R^K×Kwith the property of Q^T1K =0K.

After this brief algebraic overview, the characterization ofFis continued. Since Fis a Metzlerian with the relationF^T1K =0K, the conditions of Taussky’s theorem (aside from the irreducibility) are fulfilled. We state that

PKFP^T_K1K =0K (44)

PKFP^T_K

i j ≥0

i= j;i,j ∈Z⁺_K

(45) wherePKis an arbitrary permutation matrix. To prove Eq. (44), consider the following equation:

F^T1K =F^TP^T_KPK1K =PKF^TP^T_K1K =

PKFP^T_K T

1K =0K (46) where the relationsF^T1K = 0K andPK1K =1K are adapted. The nonnegativity of the off-diagonal entries can be easily verified by performing the multiplications in Eq. (45)(i = j;i,j ∈Z⁺_K):

PKFP^T_K

i j = K k=1

K l=1

PKi k f_kl

P^T_K

l j = K k=1

K l=1

δ_π(i)kf_klδ_π(j)l = f_{π( i)π(j)}≥0 (47) where the permutation π : Z⁺_K ↔ Z⁺_K connected to the matrixPK is introduced and the inequality f_{π( i)π(j)} ≥ 0 is applied becauseπ (i)=π (j)in case ofi = j.

Choosing a suitable permutation matrixPK, we can get the Frobenius normal form of matrixF:

PKFP^T_K =

⎛

⎜⎜

⎜⎝ F₁₁ F₂₁ F₂₂

... ... ...

F_NT1F_NT2. . .F_NTNT

⎞

⎟⎟

⎟⎠ (48)

whereF_ii ∈ R^Zi×Zi,Zi >1(i ∈ Z⁺_NT)andNT ≥1 is the number of the diagonal blocks in Eq. (48) which is equal to the number of the strongly connected components in the reaction graph associated to the network

D,G,A,k [17]. The Frobenius normal form of the matrixFis also generated by the matrixPK(Fii ∈R^Zi×Zi;i ∈Z⁺_NT):

PKFP^T_K =

⎛

⎜⎜

⎜⎝ F11

F21 F22

... ... ...

FN_T1FN_T2. . .FN_TN_T

⎞

⎟⎟

⎟⎠. (49)

As the matrixPKFP^T_K is Metzlerian, soF_{j i} ∈[0,∞)^Zj×Zi(i< j;i,j∈Z⁺_NT)and Taussky’s theorem holds for the irreducible matricesF_ii(det

F_ii

=0⇔F^T_ii1Zi = 0Zi). Let us suppose that det

F_ii

= 0 for somei ∈ Z⁺_NT. In this case, Eq. (46) involves the subsequent connection:

N_T

j=i

F^Tj i1Zj =F^Tii1Zi +

N_T

j=i+1

F^Tj i1Zj =

N_T

j=i+1

F^Tj i1Zj =0Zi. (50)

which involves the relationsF_{j i} =0Z_j×Z_i(i < j;i,j ∈Z⁺_NT).F_iiis calledterminal if det

F_ii

=0, otherwise it isnonterminal. (These notions are analogous to the ones of the terminal and nonterminal strongly connected components from [35].) If the diagonal blockF_ii is terminal, then zero is a simple eigenvalue, namely,F_ii obeys Hearon’s theorem. Put it another way,F_ii is singular with a simple zero eigenvalue iff it is a terminal block. Similar results are achieved by Foster [19], but in a slightly more complicated way.

The singularity of the blocksFii can be extended to the blocksFii which case is generally not treated in the literature. If the matrixF_ii is terminal, then det(Fii)=0 andFj i = 0Zj×Zi(i < j;i,j ∈ Z⁺_NT), involving that the multiplicity of the zero eigenvalue in the matrixFis equal to the number of the terminal blocks in its Frobenius normal form.

It is well-known from the Abel–Ruffini theorem that the roots of the characteristic polynomial associated to the matrixFii(i ∈ Z⁺_NT)can be symbolically expressed in general case iff its degree is not greater than four. Consequently, the eigenvalues of an arbitrary matrixFcan be given by closed formulas if Zi ≤ 4 or Zi ≤ 5 and det(Fii)=0(i ∈Z⁺_NT).

If either of the conditions above is in force for some blockFii, then it is advan- tageous to find its eigenvalues with the (linear, quadratic, cubic and quartic) root formulas. This is due to the fact that the convergence of the numerical eigenvalue algorithms is not guaranteed in certain cases. (These methods may exhibit oscillating or chaotic behavior.)

4 Simulation

As a real chemical application, an industrially important ten-component network, the mechanism for reductive dechlorination and elimination of chlorinated ethene sys- tem [15] is modeled (Fig.1) whose species have algebraically closed concentration functions. In this network, denoted byN, the reactions can be treated as formal isomer- ization processes whose rate coefficients are experimentally determined in aqueous solution on zinc catalysts with a specific surface area of 1.0 m²cm⁻³[15]. The values