Clustering-based H2 -State Aggregation of Positive Networks and
Its Application to Reduction of Chemical Master Equations
Takayuki Ishizaki†1 , Kenji Kashima2 , Antoine Girard3 , Jun-ichi Imura1 , Luonan Chen4 and Kazuyuki Aihara5

Abstract— In this paper, based on a notion of network
clustering, we propose a state aggregation method for positive
systems evolving over directed networks, which we call positive
networks. In the proposed method, we construct a set ofclusters (i.e., disjoint sets of state variables) according to akind of local uncontrollability of systems. This method preservesinterconnection topology among clusters as well as stability andsome particular properties, such as system positivity and steadystate characteristic (steady-state distribution). In addition, we
derive an H2 -error bound of the state discrepancy caused by theaggregation. The efficiency of the proposed method is shown
through the reduction of a chemical master equation representing the time evolution of the Michaelis-Menten chemical
reaction system.

I. I NTRODUCTION
Dynamical systems evolving on large-scale networks,
whose behaviors are determined by the interaction of a number of subsystems, have been widely studied over the past
decades. Examples of the network systems include social
networks, biology networks, spread of infection; see [1] foran overview. For such systems, it is crucial to address a
network structure preserving model reduction problem. This
structured reduction enables us to analyze coarse propertiesof the given large-scale networks such as averaged behavior.
As one of possible approaches, singular perturbation-basedstate aggregation techniques have been studied in, e.g., [2],[3]. However, this kind of approach does not explicitly take
into account the effect of external inputs. Furthermore, manykinds of structure-preserving model reduction methods havebeen considered in literature. For example, the papers [4],[5] address a model reduction problem while preserving
some underlying structure of systems such as the Lagrangian
structure and the second-order structure. However, these
problems are not formulated based on the premise of the
network structure. In addition, even though the paper [6]has proposed a network structure preserving reduction for† Research Fellow of the Japan Society for the Promotion of Science
1 Department of Mechanical and Environmental Informatics, Graduate
School of Information Science and Engineering, Tokyo Institute of Technology; 2-12-1, Meguro, Tokyo:{ishizaki@cyb., imura@}mei.titech.ac.jp
2 Department of Systems Innovation, Graduate School of Engineering
Science, Osaka University; 1-3, Machikaneyama, Toyonaka, Osaka:kashima@sys.es.osaka-u.ac.jp
3 Laboratory Jean Kuntzmann, University of Grenoble; B.P. 53, 38041
Grenoble, France: Antoine.Girard@imag.fr
4 Affiliatioin with position: Chinese Academy of Sciences, The Universityof Tokyo: lnchen@sibs.ac.cn
5 Institute of Industrial Science, University of Tokyo; 4-6-1 Komaba,
Meguro ward, Tokyo: aihara@sat.t.u-tokyo.ac.jp

network systems, a priori knowledge on partitioning (clustering) of the subsystems is required. Moreover, the relation
between partitioning and the reduction error is not discussed.
In contrast with these existing approaches, this paper
proposes a state aggregation method for positive systems
evolving on directed networks, which we call positive networks, based on a notion of network clustering. In thismethod, we construct a set of clusters (i.e., disjoint subsets
of state variables) according to local uncontrollability ofsystems. The aggregation of the constructed clusters undersuitable weights provides a reduced model that preservesinterconnection topology among clusters as well as stability
and some particular properties, such as system positivity andsteady-state characteristic. In addition, we derive an H2 error bound of the state discrepancy caused by the aggregation. The method proposed in this paper coincides witha generalization of our state aggregation approach proposedin [8], where systems evolving on undirected networks are
considered and the resultant aggregation error is evaluatedin terms of the H∞ -norm.
Furthermore, this paper also provides an application to
the reduction of chemical master equations (CMEs), which
describe the time evolution of chemical reactions by a set of
linear ordinary differential equations. The CME expressionof chemical reactions belongs to a class of continuous-timeMarkovian processes, where each state variable represents
the realization probability of a molecule distribution. The
state transition diagram can be regarded as a positive network on a multidimensional lattice (see Section IV for thedetails). The proposed method provides a reduced model that
preserves specific properties of the CMEs such as system
positivity and steady-state distribution.
This paper is organized as follows: In Section II, we
describe a positive system evolving on directed networks,
called positive networks, and formulate a clustering-basedstate aggregation problem for this class of systems. In Section III, using a projected controllability gramian of positive
networks, we devise a state aggregation method where weconstruct clusters redundant in a suitable sense. In addition,
we provide an algorithm to construct the redundant clusters.
Section IV shows the efficiency of the proposed method byapplying it to the reduction of chemical master equations,which belong to a class of the positive networks. Finally,
Section V concludes this paper.
NOTATION The following notation is to be used: R: theset of real numbers, R+ (R0+ ): the set of positive (nonnegative) real numbers, N: the set of natural numbers, In :the unit matrix of the size n × n, eni : the ith column vectorof In , span(M ): the range space spanned by the column
vectors√ of a matrix M, ∥M ∥F : the Frobenius norm of M ,i.e., tr[M M ∗ ], ∥M ∥: the induced 2-norm of a matrix M ,i.e., the maximum singular value of M , diag(v): the diagonal
matrix whose diagonal entries are the entries of a vectorv, Diag(M1 , . . . , Mn ), the block diagonal matrix whose
diagonal blocks are composed of matrices M1 , . . . , Mn .
For a set of natural numbers I, the cardinality of I isdenoted by |I|. Furthermore, enI ∈ Rn×|I| denotes the matrix
whose column vectors are composed of eni for each i ∈ I,
i.e., enI = [eni1 , . . . , enim ] ∈ Rn×m for I = {i1 , . . . , im }.
A matrix M (respectively, a transfer matrix G) is said
to be marginally stable if its eigenvalues (poles) are allin the closed left-half plane, and all eigenvalues (poles)
with zero real value are simple roots. In particular, M(respectively, G) is said to be stable if they are in the
open left-half plane. Furthermore, M is said to be reducible
if it can be placed into block upper-triangular form by
simultaneous row and column permutations. Conversely,
M is said to be irreducible if it is not reducible. Finally, the H2 and H∞ -norm of∫ a stable transfer matrix∞
1/2
G are denoted by ∥G(s)∥H2 := ( −∞ ∥G(jω)∥2F dωand
2π )∥G(s)∥H∞ := supω∈R ∥G(jω)∥, respectively.
II. P ROBLEM F ORMULATIONIn this paper, we deal with the following positive linear
systems evolving on directed networks:
Definition 1: Define
Mn := {M = {mi,j } ∈ Rn×n : mi,j ≥ 0, ∀i ̸= j,irreducible, marginally stable}. (1)A linear systemẋ = Ax + bu

(2)

is said to be positive network (A, b) if A ∈ Mn and b ∈ Rn0+ .
In this paper, for simplicity of discussion, we only deal
with single-input systems, even though the generalization tomulti-input systems is available. Furthermore, the assumptionof the irreducibility in (1) can be relaxed under a suitablesituation (see Section III-D below for the details).
This class of systems includes, e.g., spatially-discrete
reaction-diffusion systems, electric circuit systems and continuous time Markovian processes. The state trajectory ofthese systems does not escape from the non-negative orthant,i.e., Rn0+ , under any non-negative input signals and initialconditions. Systems having such a non-negative property,
called positive systems, often appear in ecology, industrial
engineering and socio-economics [9], [10]. For A = {ai,j }and b = {bi }, Fig. 1 depicts the interconnection topology ofstates (network structure) of positive networks.
In matrix theory, the largest real part eigenvalue of theMetzler matrices is called the Frobenius eigenvalue, and theassociated eigenvector is called the Frobenius eigenvector,
whose entries are all non-negative [9], [10]. In what follows,
λF (M ) denotes the Frobenius eigenvalue of M ∈ Mn .
nFurthermore, vl (M ) ∈ R1×n
0+ and vr (M ) ∈ R0+ denote theFig. 1.

Depiction of Positive Networks (ai,j ≥ 0, i ̸= j).

left and right Frobenius eigenvectors such that ∥vl (M )∥ =∥vr (M )∥ = 1. First, we state the following lemma that sums
up the properties of A ∈ Mn :Lemma 1: Any A ∈ Mn is equipped with the following
properties:• the Frobenius eigenvalue λF (A) is simple and real
• the entries of both left and right Frobenius eigenvectors
are all positive, i.e., vlT (A), vr (A) ∈ Rn+• A has no eigenvalue on the imaginary axis except theorigin.
This lemma ensures that if A ∈ Mn has the zeroeigenvalue λF (A) = 0, its left and right eigenspaces arenecessarily one-dimensional. Furthermore, defineM†n := {M ∈ Mn : vr (M ) = vlT (M )}.
The following diagonal transformation is to be useful forsystem analyses:
Lemma 2: Given A ∈ Mn , define]1[(3)
D := diag−1 (vr (A)) diag(vlT (A)) 2 .
Then, DAD−1 ∈ M†n holds.
This lemma guarantees that any A ∈ Mn is diagonallysimilar to DAD−1 ∈ M†n . Note that DAD−1 has a same
network topology (Boolean structure) as that of A. The matchof the left and right Frobenius eigenvectors of A ∈ M†nworks important roles in the following arguments. Based on
this lemma, we, as necessary, assume A ∈ M†n without lossof generality. Next, in order to formulate a clustering-basedstate aggregation problem, we give the following definition:
Definition 2: Let (A, b) be a positive network. The familyof an index set {I[l] }l∈L for L := {1, . . . , L} is called acluster set (its element is referred to as a cluster) if eachelement
I[l] is a disjoint subset of {1, . . . , n} and it satisfies∪l∈L I[l] = {1, . . . , n}. Furthermore, an aggregation matrix
compatible with {I[l] }l∈L is defined by
P := Diag(p[1] , . . . , p[L] )Π ∈ RL×n

(4)

with p[l] ∈ R1×|I[l] | such that ∥p[l] ∥ = 1, and the permutation
matrix Π := [enI[1] , . . . , enI[L] ]T ∈ Rn×n for enI[l] ∈ Rn×|I[l] | .
Then, the aggregated model of (A, b) associated with P isgiven by (PAPT , Pb).
The aggregation matrix P clearly satisfies PPT = IL . This
clustering-based state aggregation represents the aggregation

CoarseAggregated Model

Posi!ve Network
Fine
Fig. 2.

Depiction of Clustering-based State Aggregation.

of the original state x[l] ∈ R|I[l] | into the aggregated statesx̂[l] ∈ R under the aggregation weight of p[l] ∈ R1×|I[l] | ;see Fig. 2. Let us consider recovering x[l] ∈ R|I[l] | from|I[l] |
T
the aggregated states by pT.
[l] x̂[l] = (p[l] p[l] )x[l] ∈ R
According to this embedding, the transfer functions of thepositive network and the aggregated model are defined byg(s) = (sIn − A)−1 b
(5)
T
T −1
ĝ(s) = P (sIL − PAP ) Pb,
respectively. In this notation, the state aggregation problem
to be considered is formulated as follows:
Problem 1: Let (A, b) be a positive network. Given a
constant ϵ ∈ R+ , find an aggregation matrix P in (4) suchthat the aggregated model (PAPT , Pb) is a positive networkand satisfies ∥g(s) − ĝ(s)∥H2 ≤ ϵ.
III. M AIN R ESULTS
A. Exact Aggregation
In what follows, we only consider the case of λF (A) = 0
because the results derived here is almost straightforwardly
applied to stable positive networks. In this situation, we are
required to take care of the stability of the error system g − ĝin Problem 1. To ensure the stability of the error system, we
provide the following form of the aggregation weight:
Lemma 3: Given a positive network (A, b), assume A ∈M†n . For any cluster set {I[l] }l∈L , ifp[l] :=

vl (A)enI[l]∥vl (A)enI[l] ∥1×|I |∈ R+ [l] ,(enI[l] )T g(s) = pT[l] ĝ[l] (s),(6)

Proof: Consider a matrix q[l] ∈ R(|I[l] |−1)×|I[l] | such
|I[l] |×|I[l] |T Tthat [pTis unitary. Define Q by replac[l] , q[l] ] ∈ R
ing p[l] in (4) with q[l] for each l ∈ L, where we allow empty
q[l] if |I[l] | = 1. Considering the coordinate transformationof the error system by the unitary matrix [PT , QT ]T ∈ Rn×n ,we haveg(s) − ĝ(s) = Ξ(s)QT Qg(s)(7)
where Ξ(s) := PT (sIL − PAPT )−1 PA + In . To prove theclaim, it suffices to show that Ξ and Qg are both stable,
instead marginally stable.
First, note that the definition of p[l] in (6) guarantees
span(vlT (A)) ⊂ span(PT ). This provides
(8)

(9)

where p[l] is defined in (6).
This redundancy represents the linear dependence of thetransfer functions corresponding to I[l] , namely the uncontrollability of the cluster state x[l] = (enI[l] )T x. In orderto characterize (9) in an algebraic manner, we introduce a
projected controllability gramian of (A, b) as shown in thefollowing lemma:
Lemma 4: Given a positive network (A, b), let W ∈R(n−1)×nthat [vlT (A), W T ] ∈ Rn×n is unitary. Define∫ ∞suchT
W AW T τ
ΦW := 0 e
W b(eW AW τ W b)T dτ . Then, the positive semi-definite matrix
Φ := W T ΦW W ∈ Rn×n

then the error system g − ĝ associated with P in (4) is stable.

vl (A)PT P = vl (A), PT Pvr (A) = vr (A),where the property of A ∈ M†n ensures the second equality.
Using (8), we have Qvr (A) = QPT Pvr (A) = 0, which implies that the eigenspace of (A, Q) associated with λF (A) =0 is unobservable. Thus, all the poles of Qg are in the open
left-half plane.
Next, we prove the stability of Ξ. Since the positivity
1×|I |of p[l] ∈ R+ [l] holds for all l ∈ L, the off-diagonal
entries of PAPT are all non-negative. Furthermore, using(8), we have λF (PAPT ) = 0 and vl (PAPT ) = vl (A)PT ,which ensure the marginal stability of PAPT ∈ ML . In
addition, (8) provides vl (PAPT )PA = vl (A)A = 0, which
implies that the eigenspace of (PAPT , PA) associated withλF (PAPT ) = 0 is uncontrollable. Hence, all the poles of Ξ
are as well in the open left-half plane.
Lemma 3 provides a suitable aggregation weight p[l] toguarantee the stability of the error system. This is basedon the fact that the eigenspace of (A, b) associated with
λF (A) = 0 is exactly preserved into that of (PAPT , Pb).
Next, let us consider a simple situation under which a cluster
is exactly redundant in the following sense:
Definition 3: Let (A, b) be a positive network. A cluster
I[l] is said to be redundant if there exists a scaler rational
function ĝ[l] such that(10)

is independent of the choice of W .
The positive semi-definite matrix Φ in (10) represents thecontrollability gramian of the stable state-space of (A, b).
The positive semi-definite matrix Φ algebraically characterizes the redundancy of clusters as follows:
Theorem 1: Given a positive network (A, b), assume A ∈M†n . Define Φ ∈ Rn×n in (10). A cluster I[l] is redundant if
and only if there exists a row vector ϕ̂[l] ∈ R1×n such that
(enI[l] )T Φ 21 = pT[l] ϕ̂[l](11)where p[l] is defined in (6) and Φ 12 denotes a Cholesky
factor such that Φ = Φ 12 ΦT1 . Furthermore, if all clusters are2
redundant, then the aggregated model (PAPT , Pb) associatedwith P in (4) is a positive network and satisfies
g(s) = ĝ(s).

(12)

Proof: [proof of (11) ⇔ (9)] Note that (11) holds
if and only if there exists q[l] ∈ R(|I[l] |−1)×|I[l] | such
T T
that [pTis unitary and q[l] (enI[l] )T Φ 12 = 0 holds.
[l] , q[l] ]Furthermore, it is also equivalent to

Theorem 2: Given a positive network (A, b), assume A ∈M†n . Define Φ ∈ Rn×n in (10). If all cluster I[l] are θredundant, then the aggregated model (PAPT , Pb) associatedwith P in (4) is a positive network and satisfies

q[l] (enI[l] )T Φ 12 {q[l] (enI[l] )T Φ 21 }T = 0.

∥g(s) − ĝ(s)∥H2 ≤Thus, we have∫ ∞Tq[l] (enI[l] )T eAτ W T W bbT W T W eA τ enI[l] qT[l] dτ = 00⇔ q[l] (enI[l] )T eAt W T W b = 0, ∀t ∈ [0, ∞)⇔ q[l] (enI[l] )T (sIn − A)−1 W T W b = 0.

∑L

l=1 |I[l] |(|I[l] | − 1) and Ξ is defined in (7).

Proof: First, ∥g − ĝ∥H2 ≤ ∥Ξ∥H∞ ∥QT Qg∥H2 follows
from (7). Using (13) and the fact that Φ is the controllability
gramian of (A, W T W b), we have1

q[l] (enI[l] )T (sIn − A)−1 W T W b = q[l] (enI[l] )T g(s).

q[l] (enI[l] )T vr (A) = q[l] (enI[l] )T PT Pvr (A) = 0(14)

follows from (8). These imply that the eigenspaces of(A, W T W b) and (A, q[l] (enI[l] )T ) associated with λF (A) = 0
are uncontrollable and unobservable, respectively. Hence,(13) follows from the elimination of the uncontrollable andunobservable state-space.
[proof of (12)] The redundancy of all clusters implies that
Qg = 0 in (7). Hence, (12) follows.
Theorem 1 implies that the redundancy of a cluster I[l]is characterized by the linear dependence among the row
vectors of Φ 12 . More specifically, the behavior of xi isalgebraically translated into the ith row vector of Φ 12 . Note,however, that the redundancy (9) or equivalently (11) is
restrictive in general. This is because, it represents a kindof local uncontrollability such that the controllable subspaceof x[l] = (enI[l] )T x ∈ R|I[l] | is one-dimensional.
B. Redundancy Relaxation
In what follows, aiming at significant order reduction,
we relax (9) through its equivalent characterization of (11).
To this end, we introduce the following relaxation of theredundancy:
Definition 4: Let (A, b) be a positive network. Define Φ ∈Rn×n in (10). A cluster I[l] is said to be θ-redundant if there
exists a row vector ϕ̂[l] ∈ R1×n such that
1

≤ |I[l] | 2 θ, θ ∈ R+

= (tr[QT QΦQT Q]) 2 = ∥QΦ 12 ∥F .

(13)

Note that (9) is equivalent to the existence of q[l] such that(13) is equal to 0. Here, vl (A)W T W b = 0 follows from[vlT (A), W T ] being unitary and

F

(16)

∥QT Qg(s)∥H2 = ∥QT Q(sIn − A)−1 W T W b∥H2

In what follows, let us prove(enI[l] )T Φ 21 − pT[l] ϕ̂[l]where α :=√α ∥Ξ(s)∥H∞ θ(15)

where p[l] is defined in (6) and Φ 21 denotes a Cholesky factor
such that Φ = Φ 12 ΦT1 .
2

In this definition, the constant θ represents the degree
1of the redundancy (the normalization |I[l] | 2 is introduced
for a technical reason). Clearly, (15) includes (11), and isequivalent to (11) if θ = 0. Next, we propose to construct acluster set {I[l] }l∈L such that all clusters are θ-redundant.
This notion yields an aggregated model as shown in thefollowing theorem:

T T
T
is
Define ∆[l] := (enI[l] )T Φ 12 − pT[l] ϕ̂[l] . Since [p[l] , q[l] ]nT 1unitary, we have q[l] (eI[l] ) Φ 2 = q[l] ∆[l] , where we allow
empty q[l] if |I[l] | = 1. The θ-redundancy implies ∥∆[l] ∥F ≤1
|I[l] | 2 θ. Thus, we have∥QΦ 12 ∥2F ≤L∑∥q[l] ∥2F ∥∆[l] ∥2F ≤

l=1

L
∑

|I[l] |(|I[l] | − 1)θ2 .

l=1

Hence, (16) follows.
This theorem shows a linear dependence between the
degree of the redundancy θ in (15) and the aggregation error
in terms of H2 -norm. This implies that we can use θ asa design parameter to regulate the coarseness of resultant
aggregated models.
C. Cluster Construction
We give an algorithm for cluster determination. When weconstruct a cluster, we first choose an index i0 , then findindices satisfying
ϕi − vi vi−1
ϕi0 F ≤ θ, ∀i ∈ I[l]0(17)where ϕi ∈ R1×n denotes the ith row vector of Φ 12 andvi ∈ R denotes the ith entry of vl (A). The condition (17)is a sufficient condition of (15), where we take ϕ̂[l] =∥vl (A)enI[l] ∥vi−1
ϕi0 . Based on (17), an algorithm to make0
clusters is constructed as follows: Suppose that, at the lth
step, θ-redundant clusters {I[1] , . . . , I[l] } ̸= {1, . . . , n} areconstructed. At the (l + 1)th step, a subsequent cluster I[l+1]is constructed by the procedure of(i) choose i0 ∈ {1, . . . , n}\{I[1] , . . . , I[l] }, not belonging
to any cluster, and set I[l+1] = {i0 }(ii) for each i ∈ {1, . . . , n}\{I[1] , . . . , I[l] }, calculate thenorm in (17)(iii) if the norm is less than θ, update I[l+1] ← {I[l+1] , i}.
Repeating this procedure yields a cluster set {I[l] }l∈L satisfying the θ-redundancy in (15). In this construction algorithm,
degrees of freedom for the choice of i0 in step (i) and the
choice of ϕ̂[l] remain.

D. Generalization to a Class of Reducible Systems
In this subsection, for marginally stable positive networks, we consider relaxing the irreducibility assumed in
Definition 1. Note that for reducible A, we cannot apply
the diagonal transformation in Lemma 2 because the left
and right Frobenius eigenvectors of A have possibly zero
entries. This spoils the property in (8), fundamentally usedin Section III-A. Conversely, if once (8) is guaranteed, theassumption of A ∈ M†n is not required in Theorem 2.
Without loss of generality, let us assume that reducible A
is in the form of


A1,1
K
∑
 ..

n×n
.
.
A= .
, n=
nk (18)
∈R
.
AK,1

···

k=1

AK,K

where the diagonal blocks Ak,k ∈ Mnk are irreducible forall k ∈ {1, . . . , K}. In addition, we assume that†
• AK,K ∈ MnK , without loss of generality• λF (A) = 0 is simple• the entries of the left eigenvector associated with λF (A)are all positive, i.e., vl (A) ∈ R1×n
+ .
Though the second and third items are not necessarily trivial
for reducible A, chemical reaction systems described by thechemical master equation, in fact, satisfy these assumptions(see Section IV for details). Under these assumptions, it isreadily verified that λF (A) = λF (AK,K ) = 0. In addition,vl (A) and vr (A) are in the form of vl (A) = [v+ , vl (AK,K )]1×(n−nK )and vr (A) = [0, vl (AK,K )]T where v+ ∈ R+
denotes some positive row vector.
We construct θ-redundant clusters imposing a particularconstraint so as to guarantee (8). Denote the index set corresponding to the kth block by Nk , i.e., Ak,k = (enNk )T AenNk
holds for each k ∈ {1, . . . , K}. Let us construct clusters
separately within N1:K−1 := {N1 , . . . , NK−1 } and NK .
More specifically, we make {I[l] }l∈L such thatl0∪

I[l] = N1:K−1 ,

l=1

L
∪

I[l] = NK(19)

l=l0 +1

for some l0 ∈ L. Such a cluster set assures that theaggregation matrix is in theform of P = Diag(P1:K−1 , PK )∑l0 × K−1
nk
k=1
and PK ∈ R(L−l0 )×nK .
where P1:K−1 ∈ R
TSince vl (AK,K )PK PK = vl (AK,K ) holds, this ensures (8).
Thus, P associated with a θ-redundant cluster set {I[l] }l∈Lsatisfying (19) provides an aggregated model satisfying (16).
IV. A PPLICATION TO R EDUCTION OF C HEMICALM ASTER E QUATIONIn this section, we deal with the following chemical reaction system, called the Michaelis-Menten system composedof the 4 kinds of molecules S1 , . . . , S4 :

S3

c1

 S3
c2
c3

→ S4 + S2 ,where the kth entry denotes the number of Sk for each k ∈{1, . . . , 4}.
Hereafter, denote each realizable distribution by ξi ∈ N4
for i ∈ {1, . . . , n}, and the realization probability of ξiat the time t by xi (t) ∈ [0, 1]. Under this denotation, thetime evolution of (20) is represented by the chemical masterequation (CME)ẋ = Ax, x(0) = [1, 0, . . . , 0]T(21)

where x := [x1 , . . . xn ]T ∈ [0, 1]n and A ∈ Rn×n is
marginally stable and Metzler (see [11] for details). Notethat the initial condition x(0) represents that the realization
probability x1 (t) of ξ1 := [N0 , N0 , 0, 0]T is 1 at t = 0.
Furthermore, the sum of all states constantly remains to
be 1 for all t ∈ [0, ∞) because each xi represents therealization probability of the distribution ξi . This implies
that A in (21) satisfies λF (A) = 0 and vl (A) = [1, . . . , 1].
Moreover, due to the irreversibility of the the second reactionin (20), A is reducible and satisfies vr (A) = limt→∞ x(t) =[0, . . . , 0, 1]T . Fig. 3 depicts the transition diagram of the
distributions of (20) in the CME expression, where the pairof (xi ; ξi ), i.e., the realization probability and the realizable
distribution, is assigned at each node i. In this figure, thehorizontal and vertical transitions correspond to the first andthe second reaction in (20), respectively.
Remark 1: In general, the transition diagrams of chemical
reactions composed of κ kinds of molecules are described bya subgraph of the multidimensional lattice Nκ . Furthermore,
A is marginally stable and Metzler, and is uniquely determined from chemical reactions to be considered. In addition,
it satisfies that λF (A) = 0 is simple, vl (A) = [1, . . . , 1]
and vr (A) = limt→∞ x(t). It should be emphasized thatsince the system dimension tends to be considerably large(O(N0κ )), the analytical/numerical treatments of CMEs are
not necessarily easy [11]. Thus, toward practical analyses, weare required to construct a mesoscopic model that describes
essential behavior of the microscopic expression in (21).
B. H2 -aggregation of Chemical Master Equation

A. Chemical Master Equation

S1 + S2

where c1 , . . . , c3 ∈ R0+ denote the reaction rate constants.
Let us consider the situation where N0 molecules of each
S1 and S2 are present at the initial time. In this situation, allrealizable distributions of the numbers of the molecules are
enumerated by the n = n(N0 ) := (N02 + 3N0 + 2)/2 kindsof vectors, i.e.

 
 



N0
N0 − 1N0 − 20
 N0   N0 − 1   N0 − 2 



, 
, 
 , . . . ,  N0  ∈ N4 ,
 0  
 

 0 
1
2
0
0
0
N0

(20)

In what follows, we show the efficiency of the proposed
method through the reduction of the CME in (21). For (20),
let the reaction rate constants be c1 = 1, c2 = 0.1 andc3 = 3. Furthermore, let the initial number of molecules be
N0 = 140, which yields the n = 10011th dimensional CME.
From Fig. 4, which shows the resultant order of aggregatedmodels versus the value of θ, we can see that the order ofProbability
Fig. 3. State Transition Diagram of Chemical Master Equation (ξ1 =[N0 , N0 , 0, 0]T , . . . , ξn(N0 ) = [0, N0 , 0, N0 ]T ).

Time
Fig. 5. State Trajectory of Original System (the solid line, 10011th order)
and Aggregated Model (the broken line, 1077th order).

Reduced Ordersystems described by the chemical master equation, we haveobtained a mesoscopic model that well approximates theoriginal system as well as preserves specific properties as
the chemical master equation.
ACKNOWLEDGMENT

Value ofFig. 4.

Order of Resulting Models versus Values of θ.

aggregated models decreases as θ increases. This means thatthe θ-redundancy appropriately captures the coarseness ofaggregated models. In particular, when we take θ = 5×10−5 ,the order of the aggregated model is L = 1077 and theresultant relative error is ∥g − g∥H2 /∥g∥H2 = 0.0239.
Furthermore, Fig. 5 shows the state trajectory of theoriginal system x ∈ Rn (the solid line, n = 10011th order)
and that of the aggregated model PT x̂ ∈ Rn (the broken
line, L = 1077th order). We can see from this figure that the
behavior of the CME in (21) is well approximated by that of
the aggregated model. In addition, the aggregated model not
only approximates the behavior of (21), but also preserves
specific properties of the CME, namely it is equipped withthe following properties:• the positivity of the system (non-negativity of the state
trajectory) is preserved• the sum of all states constantly remains to be 1
• the steady-state distribution is exactly preserved• each state of the aggregated model represents averaged
states of the original CME
where the second to fourth items are assured by (8).
V. C ONCLUSION
In this paper, we have proposed a model reduction method
for positive systems evolving on directed networks, called
positive networks. In this method, we construct a set ofclusters based on a kind of local controllability of the statespace. Aggregating the cluster set under suitable weighting,
we obtain a reduced model that preserves interconnection
topology among the clusters as well as some specific properties, such as stability, steady-state characteristic and system
positivity. Furthermore, we have derived an H2 -error boundof the state discrepancy caused by the aggregation. In addition, applying the proposed method to chemical reaction

This research is partially supported by the Aihara Innovative Mathematical Modelling Project, the Japan Society
for the Promotion of Science (JSPS) through the “Funding
Program for World-Leading Innovative R&D on Science andTechnology (FIRST Program),” initiated by the Council for
Science and Technology Policy (CSTP). This research is also
supported by Grant-in-Aid for JSPS Fellows. The authors
would like to thank D. Cheong from Tokyo Institute ofTechnology for fruitful discussion.
R EFERENCES[1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang,“Complex networks: Structure and dynamics,” Physics Reports, vol.
424-4-5, pp. 175–308, 2006.
[2] N. Aoki, “Control of large-scale dynamic systems by aggregation,”
IEEE Transactions on Automatic Control, vol. 13-3, pp. 246–253,
1968.
[3] J. Chow and P. Kokotovic, “Time scale modeling of sparse dynamic
networks,” IEEE Transactions on Automatic Control, vol. 30-8, pp.
714–722, 1985.
[4] S. Lall, P. Krysl, and J. E. Marsden, “Structure-preserving model
reduction for mechanical systems,” Physica D: Nonlinear Phenomena,
vol. 184-1-4, pp. 304–318, 2003.
[5] R. C. Li and Z. Bai, “Structure-preserving model reduction,” Applied
Parallel Computing, vol. 3732, pp. 323–332, 2006.
[6] H. Sandberg and M. Murray, “Model reduction of interconnected linear
systems,” Optimal Control Applications and Methods, vol. 30-3, pp.
225–245, 2009.
[7] E. Yeung, J. Gonalves, H. Sandberg, and S. Warnick, “Network
structure preserving model reduction with weak a priori structural
information,” in Joint 48th IEEE Conference on Decision and Control
and 28th Chinese Control Conference, 2009, pp. 3256–3263.
[8] T. Ishizaki, K. Kashima, J. Imura, and K. Aihara, “Reaction-diffusion
clustering of single-input dynamical networks,” in 50th IEEE Conference on Decision and Control, and European Control Conference,
2011, pp. 7837–7842.
[9] L. Farina and S. Rinaldi, Positive Linear Systems: Theory and Applications, P. Hilton and H. Hochstadt, Eds.
Wiley-Interscience
Publication, 2000.
[10] L. Benvenuti and L. Farina, “A tutorial on the positive realization
problem,” IEEE Transactions on Automatic Control, vol. 49-5, pp.
651–664, 2004.
[11] B. Munsky and M. Khammash, “The finite state projection approachfor the analysis of stochastic noise in gene networks,” IEEE Transactions on Automatic Control, vol. 53, Special Issue, pp. 201–214,
2008.