

HUNGARIAN JOURNAL
OF INDUSTRY AND CHEMISTRY

VESZPRÉM
Vol. 42(2) pp. 103–107 (2014)

ON THE PARAMETRIC UNCERTAINTY OF WEAKLY REVERSIBLE
REALIZATIONS OF KINETIC SYSTEMS

GYÖRGY LIPTÁK B1 , GÁBOR SZEDERKÉNYI1,2 AND KATALIN M. HANGOS1,3

1Process Control Research Group, MTA SZTAKI, Kende u. 13-17, Budapest, 1111, HUNGARY
2Faculty of Information Technology, Péter Pázmány Catholic University, Práter u. 50/a, Budapest, 1083, HUNGARY

3Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10,
Veszprém, 8200, HUNGARY
BE-mail: lipgyorgy@gmail.com

The existence of weakly reversible realizations within a given convex domain is investigated. It is shown that the domain of
weakly reversible realizations is convex in the parameter space. A LP-based method of testing if every element of a convex
domain admits weakly reversible realizations is proposed. A linear programming method is also presented to compute
a stabilizing kinetic feedback controller for polynomial systems with parametric uncertainty. The proposed methods are
illustrated using simple examples.

Keywords: parametric uncertainty; computational methods; optimization; kinetic systems

Introduction

The notion of parametric robustness is well-known and
central in linear and nonlinear systems and control theory
[1]. It is used for ensuring a desirable property, such as
stability, in a given domain in the parameter space around
a nominal realization having the desired property.

The aim of the paper is to extend the notions and tools
of parametric robustness for a class of positive polyno-
mial systems, namely a class of kinetic systems. Only the
very first steps are reported here that offer a computation-
ally efficient method for checking one of the many impor-
tant properties of kinetic systems, their weak reversibility.

Basic Notions and Methods

The basic notions and tools related to reaction kinetic sys-
tems and their realizations are briefly summarized in this
section.

Kinetic Systems, their Dynamics and Structure

Deterministic kinetic systems with mass action kinetics
or simply chemical reaction networks (CRNs) form a
wide class of non-negative polynomial systems, that are
able to produce all the important qualitative phenomena
(e.g. stable/unstable equilibria, oscillations, limit cycles,
multiplicity of equilibrium points and even chaotic be-
haviour) present in the dynamics of nonlinear processes
[2].

The general form of dynamic models studied in this
paper is the following

ẋ = M ·ψ(x), (1)

where x ∈ Rn is the state variable and M ∈ Rn×m. The
monomial vector function ψ : Rn → Rm is defined as

ψj(x) =

n∏
i=1

x
Yij
i , j = 1, . . . ,m (2)

where Y ∈ Nn×m0 . The system Eq.(1) is kinetic if and
only if the matrix M has a factorization

M = Y ·Ak. (3)

The Kirchhoff-matrix Ak has non-positive diagonal and
non-negative off-diagonal elements and zero column
sums. The matrix pair (Y,Ak) is called the realization
of the system Eq.(1).

The chemically originated notions: The chemically
originated notions of kinetic systems are as follows: the
species of the system are denoted by X1, . . . ,Xn, and
the concentrations of the species are the state variables
of Eq.(1), i.e. xi = [Xi] ≥ 0 for i = 1, . . . ,n. The
structure of kinetic systems is given in terms of its com-
plexes Ci, i = 1, . . . ,m that are non-negative linear
combinations of the species i.e. Ci =

∑n
j=1[Y ]jiXj for

i = 1, . . . ,m, and therefore Y is also called the complex
composition matrix.

The reaction graph: The weighted directed graph (or
reaction graph) of kinetic systems is G = (V,E), where


104

V = {C1,C2, . . . ,Cm} and E denote the set of ver-
tices and directed edges, respectively. The directed edge
(Ci,Cj) (also denoted by Ci → Cj ) belongs to the re-
action graph if and only if [Ak]j,i > 0. In this case,
the weight assigned to the directed edge is Ci → Cj is
[Ak]j,i.

Stoichiometric subspace: Stoichiometric subspace S
is given by the span of the reaction vectors

S = {[Y ]·i − [Y ]·j | [Ak]ij > 0}. (4)

The stoichiometric compatibility classes of a kinetic sys-
tem are the affine translations of the stoichiometric sub-
space: (x0 + S)∩Rn≥0.

Structural Properties and Dynamical Behaviour

It is possible to utilize certain structural properties of ki-
netic systems that enable us to effectively analyze the sta-
bility of the system.

Deficiency: There are several equivalent ways to de-
fine deficiency. We will use the following definition

δ = dim(ker(Y )∩ Im(BG)), (5)

where BG is the incidence matrix of the reaction graph.
It is easy to see that deficiency is zero if ker(Y ) = {0}
or equivalently rank(Y ) = m.

Weak reversibility: A CRN is called weakly reversible
if whenever there exists a directed path from Ci to Cj in
its reaction graph, then there exists a directed path from
Cj to Ci. In graph theoretic terms, this means that all
components of the reaction graph are strongly connected
components.

Deficiency zero theorem: A weakly reversible kinetic
system with zero deficiency has precisely one equilibrium
point in each positive stoichiometric compatibility class
that is locally asymptotically stable (conjecture: globally
asymptotically stable).

Computing Weakly Reversible Realizations Formulated
as an Optimization Problem

In this section, first a method for computing weakly re-
versible realization based on Ref.[3] is briefly presented.
We assume that we have a kinetic polynomial system of
the form Eq.(1).

We use the fact known from the literature that a real-
ization of a CRN is weakly reversible if and only if there
exists a vector with strictly positive elements in the kernel
of Ak, i.e. there exists b ∈ Rn+ such that Ak · b = 0 [4].
Since b is unknown, too, this condition in this form is not
linear. Therefore, we introduce a scaled matrix Ãk

Ãk = Ak ·diag(b) (6)

where diag(b) is a diagonal matrix with elements of b.
It is clear from Eq.(6) that Ãk is also a Kirchhoff matrix

and that 1 ∈ Rm (the m-dimensional vector containing
only ones) lies in kernel of Ãk. Moreover, it is easy to see
that Ãk defines a weakly reversible network if and only
if Ak corresponds to a weakly reversible network. Then,
the weak reversibility and the Kirchhoff property of Ãk
can be expressed using the following linear constraints

m∑
i=1

[
Ãk

]
ij

= 0, j = 1, . . . ,m

m∑
i=1

[
Ãk

]
ji

= 0, j = 1, . . . ,m

[
Ãk

]
ij
≥ 0, i,j = 1, . . . ,m, i 6= j[

Ãk

]
ii
≤ 0, i = 1, . . . ,m. (7)

Moreover, the equation Eq.(3) is transformed by diag(b)
(we can do this, because diag(b) is invertible):

M ·diag(b) = Y ·Ak ·diag(b)︸ ︷︷ ︸
Ãk

(8)

Finally, by choosing an arbitrary linear objective function
of the decision variables Ãk and b, weakly reversible real-
izations of the studied kinetic system can be computed (if
any exist) in a LP framework using the linear constraints
Eq.(7) and (8).

Weakly Reversible CRN Realizations

In this section, first the convexity of the weakly reversible
Kirchhoff matrix will be shown. After that the practical
benefits of this property will be demonstrated in the field
of system analysis and robust feedback design.

Convexity of the Weak Reversibility in the Parameter
Space

Theorem 1. Let A(1)k and A
(2)
k be m × m weakly re-

versible Kirchhoff matrices. Then the convex combination
of the two matrices remains weakly reversible.

Proof. The idea behind the proof is based on Ref.[5].
A Kirchhoff matrix is weakly reversible if and only if
there is a strictly positive vector in its kernel. Therefore
strictly positive vectors p1, p2 exist such as A

(1)
k ·p1 = 0

and A(2)k · p2 = 0. Let us define the following scaled
Kirchhoff matrix: Â(1)k = A

(1)
k · diag(p1) and Â

(2)
k =

A
(2)
k · diag(p2). These scaled matrices have identical

structures to the original ones. Moreover, Â(1)k ·1
(m) = 0

and Â(2)k ·1
(m) = 0 where the vector 1(m) denotes the m

dimensional column vector composed of ones. For that

(λÂ
(1)
k + (1−λ)Â

(2)
k ) ·1

(m) = 0. (9)

for any λ ∈ [0,1]. Therefore the convex combination
of the original two realizations has to be weakly re-
versible.


105

(a) (b) (c)

Figure 1: A weakly reversible reaction graphs of the
three realizations (Y,A(1)k ), (Y,A

(2)
k ) and (Y,A

(3)
k )

Weak Reversibility of CRN Realizations with Parametric
Uncertainty

We assume that a CRN with parametric uncertainty is
given as

ẋ = M ·ψ(x), (10)
where x ∈ Rn is the state variable, ψ ∈ Rn → Rm
contains the monomials and the matrix M ∈ Rn×m is an
element of the following set

M =
{ l∑
i=1

αiMi | (∀i : αi ≥ 0) ∧
l∑

i=1

αi = 1
}
.

(11)
The goal is to find a method for checking the weak re-
versibility of the system Eq.(10) for all matrices M ∈M.

When all vertices Mi have a weakly reversible real-
ization (Y,A(i)k ) then any element of the set M has a re-
alization (Y,Ak) such that Ak is the convex combination
of the Kirchhoff matrices A(i)k . The obtained realization
Ak will be weakly reversible due to Theorem 1. There-
fore, it is enough to compute a weakly reversible realiza-
tion for each matrix Mi by using the previously presented
LP-based method.

A Simple Example

Let us consider the following polynomial system

[
ẋ1
ẋ2

]
= M ·


 x1x2
x1x2


 , (12)

where M is an arbitrary convex combination of the fol-
lowing three matrices

M1 =

[
0 1 −1
1 −1 0

]
,

M2 =

[
−1 1 0
1 −1 0

]
,and

M3 =

[
0 1 −1
0 0 0

]
.

Figure 2: A weakly reversible realization of the convex
combination M = 0.2M1 + 0.4M2 + 0.4M3

In order to show weak reversibility for all possible con-
vex combinations, we have to find a weakly reversible re-
alization for each matrix M1, M2 and M3. The resulting
weakly reversible reaction graphs are depicted in Fig.1,
while Fig.2 illustrates an inner point realization which is
weakly reversible too.

Computing Kinetic Feedback for a Polynomial
System with Parametric Uncertainty

Besides the possible application of the above described
LP-based method for robust stability analysis, it can also
be used for stabilizing feedback controller design. For
this purpose, a generalized version of our preliminary
work on kinetic feedback computation for polynomial
systems to achieve weak reversibility and minimal defi-
ciency [6] is used here.

The Feedback Design Problem

We assume that the equation of the open-loop polynomial
system with linear constant parameter input structure is
given as

ẋ = M ·ψ(x) + Bu, (13)

where x ∈ Rn is the state vector, u ∈ Rp is the input and
ψ ∈ Rn → Rm contains the monomials of the open-loop
system. The input matrix is B ∈ Rn×p, the correspond-
ing complex composition matrix is Y with rank m, and
M ∈ Rn×m is an element of the following set

M =

{
l∑

i=1

αiMi | (∀i : αi ≥ 0) ∧
l∑

i=1

αi = 1

}
.

(14)
Moreover, a positive vector x ∈ Rn>0 being the desired
equilibrium point is given as a design parameter. Note
that the above polynomial system is not necessarily ki-
netic, i.e. not necessarily positive, and may not have a
positive equilibrium point at all.


106

The aim of the feedback is to set a region in the state
space R ⊆ Rn≥0 where x is (at least) a locally asymptot-
ically stable equilibrium point of the closed-loop system
for all M ∈M.

For this purpose we are looking for a feedback in the
form

u = Kψ(x) (15)

which transforms the open-loop system into a weakly re-
versible kinetic system with zero deficiency for all M ∈
M with the given equilibrium point x.

Feedback Computation

Similarly to the realization computation, the matrix K
will be determined by solving an LP problem. The con-
vexity result shows that it is enough to compute one
weakly reversible realization (Y,A(r)k ) in each vertex Mr
to ensure weak reversibility for all possible closed-loop
systems. All realizations will have zero deficiency, be-
cause of the rank condition rank(Y ) = m [7].

First we note, that the realization (Y,A(r)k ) that corre-
sponds to the closed-loop system is

Mr + B ·K = Y ·A
(r)
k . (16)

where the matrix A(r)k should be Kirchhoff

m∑
i=1

[Ãk]ij = 0, j = 1, . . . ,m

[Ãk]ij ≥ 0, i,j = 1, . . . ,m, i 6= j

[Ãk]ii ≤ 0, i = 1, . . . ,m. (17)

In order to obtain a weakly reversible closed-loop system
with an equilibrium point x, the matrix A(r)k should be
weakly reversible and has to have the vector ψ(x) in its
right kernel, i.e.

A
(r)
k ·ψ(x) = 0. (18)

Finally, by choosing an arbitrary linear objective func-
tion of the decision variables A(1)k , . . . ,A

(l)
k and K, the

feedback gain K can be computed (if it exists) in a LP
framework using the linear constraints Eqs.(16-18).

With the resulting feedback gain K, the point x will
be an equilibrium point of all possible closed-loop sys-
tems, and x will be locally asymptotically stable in the
region S = (x+S)∩Rn≥0, where S is the stoichiometric
subspace of the closed-loop system.

Example

Let the open-loop system be given as

ẋ = M


 x1x2x2x3

x1


 +


 01

0


u (19)

Figure 3: Weakly reversible realization of the
closed-loop system, where

M = 0.6M1 + 0.2M2 + 0.2M3

where M is an arbitrary convex combination of the fol-
lowing three matrices:

M1 =


 −1 1 02 1 2

1 −1 0


 ,

M2 =


 0 0 01 1 3

0 0 0


 ,and

M3 =


 0 1 −12 0 3

0 −1 1


 .

The desired equilibrium point x = [1 1 1]T .
We are looking for a feedback law with gain K which

transforms the matrices Mi into weakly reversible kinetic
systems with the given equilibrium point.

By solving the feedback design LP optimization prob-
lem using the linear constraints Eqs.(16-18), the com-
puted feedback is in the following form:

u =
[
2 1 2

]
ψ(x). (20)

Fig.3 depicts a weakly reversible realization of the
closed-loop system. The obtained closed-loop system in
an inner point of the convex set M has the following sto-
ichiometric subspace:

S = span




 11

0


−


 01

1


 ,

 11

0


−


 10

0




 .
(21)

Therefore, the equilibrium point x will be asymptotically
stable with the region S = (x+S)∩Rn≥0. Note, that one
should choose the initial value of the state variables from
S.

Fig.4 shows the time dependent behaviour of the
closed-loop solutions started from different initial points
in S.


107

Figure 4: Time-domain simulation of the closed-loop
system

Conclusion

It is shown in this paper that the domain of weakly re-
versible realizations is convex in the parameter space.
This property is utilized for developing methods in sys-
tem analysis and robust control design. An LP-based op-
timization method is proposed for testing if every element
of a convex domain given by its extremal matrices admits
a weakly reversible realization. An LP-based feedback
design method is also proposed that guarantees stability
with a desired equilibrium point. The proposed methods
are illustrated with simple examples.

REFERENCES

[1] SILJAK D.: Parameter space methods for robust con-
trol design: a guided tour, Aut. Contr., IEEE Trans-
actions on, 1989, 34(7), 674–688

[2] ANGELI D.: A tutorial on chemical network dynam-
ics, Eur. J. Contr., 2009, 15, 398–406

[3] SZEDERKÉNYI G., HANGOS K. M., TUZA Z.: Find-
ing weakly reversible realizations of chemical re-
action networks using optimization, MATCH Com-
mun. Math. Comp. Chem., 2012, 67, 193–212

[4] FEINBERG M., HORN F.: Chemical mechanism
structure and the coincidence of the stoichiometric
and kinetic subspaces, Arch. Rational Mechanics and
Analysis, 1977, 66(1), 83–97

[5] RUDAN J., SZEDERKÉNYI G., HANGOS K.M.,
PÉNI T.: Polynomial time algorithms to determine
weakly reversible realizations of chemical reaction
networks, J. Math. Chem., 2014, 52(5), 1386–1404

[6] LIPTÁK G., SZEDERKÉNYI G., HANGOS K.M.: Ki-
netic feedback computation for polynomial systems
to achieve weak reversibility and minimal deficiency,
13th Eur. Control Conf. ECC, Strasbourg, France,
2014, 2691–2696

[7] ANGELI D., LENHEE, P.D., SONTAG E.D.: On
the structural monotonicity of chemical reaction
networks, Proc. 45th IEEE Conf. Decision and
Control, 2006, 7–12