- Email: [email protected]
Steady-state analysis of probabilistic Boolean networks
Available online at www.sciencedirect.com
Journal of the Franklin Institute 356 (2019) 2994–3009 www.elsevier.com/locate/jfranklin
Steady-state analysis of probabilistic Boolean networks Jinfeng Pan a, Jun-e Feng a,∗, Min Meng b b School
a School of Mathematics, Shandong University, Jinan 250100, China of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798, Singapore
Received 20 April 2018; received in revised form 18 January 2019; accepted 27 January 2019 Available online 1 February 2019
Abstract This paper investigates steady-state distributions of probabilistic Boolean networks via cascading aggregation. Under this approach, the problem is converted to computing least square solutions to several corresponding equations. Two necessary and sufficient conditions for the existence of the steady-state distributions for probabilistic Boolean networks are given firstly. Secondly, an algorithm for finding the steady-state distributions of probabilistic probabilistic Boolean networks is given. Finally, a numerical example is given to show the effectiveness of the proposed method. © 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
1. Introduction Developing models to analyse gene regulatory networks is one of the aims in systems biology. Numbers of mathematical models have been proposed to study gene regulatory networks, such as Bayesian networks [1], Boolean networks (BNs) [2], Petri nets [3], probabilistic Boolean networks (PBNs) [4], linear regression model [5]. Among them, BNs and PBNs, which were originally introduced by Kauffman [2], are some of the most attractive models. One BN is defined by a set of nodes (genes) {x1 , x2 , . . . , xn }, xi ∈ {0, 1}, and a set of Boolean functions { f1 , f2 , . . . , fn }, with the update rules as xi (t + 1) = fi (x j (t ), j ∈ Ni ). Hence, for a ∗
Corresponding author. E-mail addresses: [email protected] (J. Pan), [email protected] (J.-e. Feng), mengminmath@ gmail.com (M. Meng). https://doi.org/10.1016/j.jfranklin.2019.01.039 0016-0032/© 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2995
BN with n nodes, there are 2n states. On the basis of update rules, the state transition graph of a BN can be uniquely determined, and in this graph, the out-degree of any state is 1. In addition, due to the uncertainty of gene regulation and errors of experimental noise, it is appropriate that one update rule of a node in a BN is randomly chosen at each time among the candidates of Boolean functions, which results in the appearance of PBNs [6]. As an extension of BNs, PBNs have been widely studied in recent years. A PBN can be viewed as a collection of BNs. For a PBN with n nodes {x1 , x2 , . . . , xn }, the update rule f(i) of xi is chosen from { fi1 , fi2 , . . . , fiq(i) } in probability. Then for the whole network, there are ni=1 q(i) BNs. In this paper, the vector-valued function fl ( fl(l ) , fl(2) , . . . , fl(n) ) is used to determine a context of the PBN. For each gene, at each time point, a BN is randomly chosen to determine the network’s transition. A PBN can be regarded as a homogeneous Markov chain. Therefore, some results of Markov chains can also be used to analyse PBNs. Recently, a new matrix product, semi-tensor product of matrices, is proposed by Cheng and his co-workers [7]. Based on it, BNs can be converted to algebraic forms, under which many interesting results were obtained [8–11]. Nowadays, the semi-tensor product is applied to solve problems for PBNs as well, such as the controllability [12–14], stability [15,16], observability [17], synchronization [18,19], optimal problems [20], steady-state distributions [21,22]. When considering the long-run behavior of a PBN, finding its steady-state distributions is an important aspect of analysing biological systems. This steady-state distribution can provide a first-order statistical information for PBNs. On this basis, the network and the impact of different genes can be understood. In addition, one can figure out how to control some genes in the network so that the whole system can evolve to the target state or the expected steadystate probability distribution. Then we can develop and research therapeutic gene intervention or gene control policy [6]. If a PBN has steady-state distribution, then in the long run, the probability that the chain is in a given state is stationary, independent of the initial state. But if the PBN has no steady-state distribution, such question need not be answered. The steady-state distributions can be obtained based on transition matrices. Up to now, several references were devoted to finding the steady-state distributions of PBNs [23–25]. Context-reduction mapping and an approximation method to calculate the steady-state distributions of context-sensitive PBNs were proposed in [23]. Monte Carlo methods [24] were used to calculate the steadystate distributions of PBNs. In [25], an approximation method for computing the steady-state distributions of PBNs was presented by neglecting some BNs with very small probabilities during the construction of the transition probabilistic matrices. In [26], the long run behavior of large-scale PBN was analyzed based on its stochastic conjunctive normal form network. In [27,28], an algorithm was proposed to find the attractors of large-scale Boolean networks by network aggregation, which reduces the computation significantly. Transition probabilistic matrix-based algorithm to find the steady-state distributions of PBNs is complicated even if the number of states is not large. Hence, in this paper, the problem to find the steady-state distributions of PBNs with special structures is also considered by using the cascading aggregation. Using semi-tensor product, the probabilistic transition matrix of each block can be obtained. The main contributions of this paper are concluded as follows: For PBNs, two equivalent conditions for the existence of the steady-state distribution are proposed on the basis of the probabilistic transition matrices. The network aggregation is firstly used to find steady-state distributions of PBNs. Based on cascading aggregation, the problem can be solved by calculating the least square solutions to several corresponding equations. By the method presented in this paper, the computational complexity is reduced.
2996
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
The rest of this paper is organized as follows. Section 2 introduces some fundamental definitions and notations. The definition and properties of Markov chains, and two necessary and sufficient conditions for the existence of the steady-state distributions for PBNs, as well as a method to compute steady-state distributions for PBNs, are derived in Section 3. In Section 4, a numerical example is given to illustrate the effectiveness of the proposed method. Section 5 makes a brief conclusion. 2. Preliminaries and problem formulation In this section, we present some properties about semi-tensor product of matrices that will be used in this paper, and the definition of network aggregation of PBNs. 2.1. Probabilistic Boolean networks A PBN can be regarded as an extension of BN with a probabilistic setting. Consider the following PBN:
⎧ x1 (t + 1) = f (1) ([x j (t )], j ∈ N1 ), ⎪ ⎪ ⎪ ⎨x2 (t + 1) = f (2) ([x j (t )], j ∈ N2 ), .. ⎪ . ⎪ ⎪ ⎩ xn (t + 1) = f (n) ([x j (t )], j ∈ Nn ),
(1)
where [x1 , x2 , . . . , xn ]T ∈ D n is the state vector, D = {0, 1}. Ni is the index set of xi . Assume that there are q(i) candidates of f(i) as fl(i) ([x j (t )], j ∈ Nl(i) ), l = 1, 2, . . . , q(i),
(2)
where Nl(i) is the index set corresponding to fl(i) , and the probability of choosing fl(i) as the update function is cl(i) , i.e., cl(i) = Prob( f (i) = fl(i) ), l = 1, 2, . . . , q(i).
(3)
The following relation q(i)
cl(i) = 1
(4)
l=1
must be satisfied. Ni , i = 1, 2, . . . , n, can be obtained as Ni :=
q(i)
Nl(i) .
(5)
l=1
Due to the stochasticity in genetic networks, the next state of gene i is determined by q(i) (i) (i) Boolean functions f1(i) , f2(i) , . . . , fq(i) , with probabilities c1(i) , c2(i) , . . . , cq(i) , respectively. For a n PBN of n genes, there are i=1 q(i) possible Boolean networks, each of which is a possible realization of the genetic network and can be considered as a context of the PBN. For the (1) (2) (n) jth context, assume the network function is given by f j = ( f j(1) , f j(2) , . . . , f j(n) ), then the n (i) probability that the jth context is selected is c j = i=1 c j(i) .
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2997
Example 2.1. Consider a PBN with three nodes. The update rules are given as follows: f = 1
f1(1) = x1 (t ) ∨ x2 (t ), c1(1) = 0.8, f2(1) = x1 (t ), c2(1) = 0.2,
f 2 = f1(2) = x1 (t ), c1(2) = 1, (3) f = x3 (t ), c1(3) = 0.5, f 3 = 1(3) f2 = ¬x3 (t ), c2(3) = 0.5. Then we have q(1) = 2, q(2) = 1, q(3) = 2, N1(1) = {1, 2}, N2(1) = {1}, N1 = N1(1) ∪ N2(1) = {1, 2}, N2 = N1(2) = {1}, N3 = N1(3) ∪ N2(3) = {3}, and Prob( f 1 = f1(1) , f 2 = f1(2) , f 3 = f1(3) ) = 0.4. 2.2. Algebraic forms of Boolean functions For ease of presentation, we first give some notations, many of which are standard in literature [7]. Rm × n is the set of m × n real matrices. For a matrix M ∈ Rm × n , the (i, j)th entry is denoted by Mij . M > 0 means Mij > 0 for all i, j. coli (M) (rowi (M )) is the ith column (row) of matrix M. col(M) (row(M )) is the set
of columns (rows) of M. For a given vector x = (x1 , x2 , . . . , xn ), the l1 norm of x is xl1 = ni=1 |xi |. δni is the ith column of the identity matrix In . n = {δni |i = 1, 2, . . . , n}. Assume that L ∈ Lm×n , L = [δmi1 δmi2 . . . δmin ], and its shorthand form is written as L = δm [i1 i2 . . . in ]. Given a set S = {c1 , c2 , . . . , cl }, 1 ≤ ci ≤ 2n , ci = cj , δ2n S is used to represent the set {δ2c1n , δ2c2n , . . . , δ2cln }. Given a matrix L ∈ Rm × n , assume
l P = { p1 , p2 , . . . , pl }, pi ∈ {1, 2, . . . , m}, we denote L(P, s) to represent k=1 L pk s , s ∈ {1, 2, . . . , n}. Assume A ∈ Rm × r , B ∈ Rn × r , the Khatri–Rao product of A and B is defined as follows: A ∗ B = [col1 (A ) col1 (B), col2 (A ) col2 (B), . . . , colr (A ) colr (B)]. To use matrix expression of logic, we identify T = 1 ∼ δ21 , F = 0 ∼ δ22 , thus D ∼ 2 , where ∼ means the logical equivalence. Subsequently, every Boolean function f : D n → D, can be equivalently expressed as a mapping f : n2 → 2 . Lemma 2.1. ([7]) Consider a Boolean function f (x1 , x2 , . . . , xn ) with xi ∈ 2 . Denote x = x1 x2 · · · xn , then there exists a unique matrix M f ∈ L2×2n , called the structural matrix of f, such that f (x1 , x2 , . . . , xn ) = M f ni=1 xi = M f x.
(6)
2.3. Aggregation of PBNs In this section, we first introduce the network aggregation, which plays a vital role in studying the steady-state distributions of PBNs. Denote χ = {x1 , x2 , . . . , xn }, then χ can be partitioned into μ blocks χ = χ1 ∪ χ2 ∪ · · · ∪ χμ , where χ i = ∅, χ i ⊂ χ , and χi ∩ χ j = ∅ for i
= j. Let χi = {xi1 , xi2 , . . . , xini }, where ni is the number of state nodes in the ith block, and μ i=1 ni = n.
2998
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Fig. 1. Example of aggregation of a network with 9 nodes into 3 Boolean control networks.
If j ∈ Ni , we say that there is a directed edge from xj to xi . All the starting nodes of the incoming and outgoing edges of block χ i constitute sets Ui = {ui1 , . . . , uimi } and Yi = {yi1 , . . . , yipi }, respectively. Set Ui (Yi ) may be empty, when there are no incoming (outgoing) edges of block χ i . Let C be the set of starting nodes whose edges cut by the aggregation, then we have C = ∪μi=1 Ui = ∪μi=1 Yi . Therefore the subnetwork i with nodes in χ i and inputs in Ui is a probabilistic Boolean control network, which can be described as follows: i :
xi j (t + 1) = f (i j) (xi1 (t ), xi2 (t ), . . . , xini (t ), ui1 (t ), ui2 (t ), . . . , uimi (t )),
(7)
for i = 1, 2, . . . , μ and j = 1, 2, . . . , ni .cl(i j) = Prob( f (i j) = fl(i j) ), l = 1, 2, . . . , q(i j). Definition 2.2. [27] The aggregation is said to be cascading if (with possible order) the groups χ i = ∪ij=1 χ j , i = 1, 2, . . . , μ,
(8)
have zero in-degree. Definition 2.3. [27] The aggregation is said to be an acyclic aggregation, if its aggregation graph contains no cycle. Lemma 2.4. [29] An aggregation is a cascading aggregation if and only if it is an acyclic aggregation. The method to find the acyclic aggregation of a network has been given in literature [27]. In what follows, an example is given to show the aggregation of PBNs.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2999
Example 2.2. Consider a PBN in Fig. 1, and the update rules are given as follows: ⎧ ⎪ x1 (t + 1) = x2 (t ), c1(1) = 1, ⎪ ⎪ ⎪ ⎪ x3 (t ) ∧ x7 (t ), c1(2) = 0.7, ⎪ ⎪ x (t + 1) = 2 ⎪ ⎪ x3 (t ), c2(2) = 0.3, ⎪ ⎪ ⎪ ⎪ ⎪ x3 (t + 1) = x1 (t ) ∨ x2 (t ), c1(3) = 1, ⎪ ⎪ ⎪ ⎪ x1 (t ) ∨ x5 (t ) ∨ x7 (t ), c1(4) = 0.5, ⎪ ⎪ x (t + 1) = ⎨ 4 x1 (t ), c2(4) = 0.5, ⎪ x5 (t + 1) = ¬x4 (t ), c1(5) = 1, ⎪ ⎪ ⎪ ⎪ x6 (t + 1) = x6 (t ), c1(6) = 1, ⎪ ⎪ ⎪ ⎪ x7 (t + 1) = x6 (t ), c1(7) = 1, ⎪ ⎪ ⎪ ⎪ ⎪ x (t ) ∨ x9 (t ), c1(8) = 0.6, ⎪ ⎪ x8 (t + 1) = 7 ⎪ ⎪ x7 (t ) ∧ x9 (t ), c2(8) = 0.4, ⎪ ⎪ ⎩ x9 (t + 1) = x6 (t ), c1(9) = 1. The node set χ = {x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 , x9 } can be partitioned into 3 blocks χ i, i = 1, 2, 3, where χ 1 = {x1 , x2 , x3 }, χ 2 = {x4 , x5 }, and χ 3 = {x6 , x7 , x8 , x9 }. Reorder the blocks, and define χ 1 = χ 3, χ 2 = χ 3 ∪ χ 1, χ 3 = χ . By Definition 2.2, the aggregation is cascading. 3. Steady-state distribution of PBNs In this section, the definitions of Markov chains and steady-state distribution as well as some properties are presented firstly. Then the method to calculate the steady-state distributions of PBNs under cascading aggregation is given in two parts. Part 1 provides two necessary and sufficient conditions for the existence of the steady-state distributions of PBNs. Part 2 gives a method of calculating the steady-state distributions of PBNs. Definition 3.1. [30] Assume that X = {Xn , n = 0, 1, . . .} is a discrete random process which is defined in (, F, P ), and its state space S is a countable set or finite set. If for every nonnegative integer n, and arbitrary state i0 , i1 , . . . , in+1 ∈ S, Prob(X0 = i0 , X1 = i1 , . . . , Xn = in ) > 0, and Prob(Xn+1 = in+1 |X0 = i0 , . . . , Xn = in ) = Prob(Xn+1 = in+1 |Xn = in ), then we say that X is a discrete Markov chain, here Prob(Xn+1 = in+1 |Xn = in ) has no relationship with n. For a given Markov chain, if its state space is finite, then its transition probability can be expressed by a probabilistic matrix P, which is called the probabilistic transition matrix, Pi j = Prob(Xn+1 = j | Xn = i).
(9)
If the Markov chain is homogeneous, then the k step probabilistic transition matrix is Pk . Definition 3.2. [30] Given a finite homogeneous Markov chain with state space S = {1, 2, . . . , M}, we say the state space has a steady-state distribution, if there exists a probability distribution π = (π1 , π2 , . . . , πM )
3000
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
such that, πj =
M
πi Pi j .
i=1
Based on Lemma 2.1, the jth context of (1) can be converted to x(t + 1) = L j x(t ). Define the probabilistic transition matrix of system (1) as follows: L = c1 L 1 + c2 L 2 + · · · + cQ L Q , where Q = ni=1 q(i), the probability of choosing the jth context is cj . From the construction of matrix L, it is obvious to see that (L) ji = Prob(x(t + 1) = δ2jn |x(t ) = δ2i n ), hence the probabilistic transition matrix of the Markov chain associated with Eq. (1) is LT . Based on the results on Markov chains in [30], it is obviously to have the following results. Theorem 3.3. Assume that LT is the transition probability matrix of Eq. (1), then (i) there exists at least one steady-state distribution if and only if the geometric multiplicity of eigenvalue 1 of L is greater than or equal to 1. (ii) there exists a unique steady-state distribution if and only if the geometric multiplicity of eigenvalue 1 of L is equal to 1. Corollary 3.4. Let L be the probabilistic transition matrix of a finite homogeneous Markov chain, (i) there exists at least one steady-state distribution if and only if rank(L − I2n ) 2n − 1. (ii) there exists a unique steady-state distribution if and only if rank(L − I2n ) = 2n − 1. 3.1. Steady-state distribution of PBNs In the above discussions, two conditions are derived to judge the existence of the steadystate distribution based on the probabilistic transition matrix. However, when the dimensions of the matrix become bigger, the methods are not applicable. In the following process, we propose a method to compute the steady-state distributions for PBNs under network aggregation. Assume PBN (1) can be partitioned to μ blocks 1 , 2 , . . . , μ by cascading aggregation, and the blocks are ordered according to Definition 2.2. Thus 1 is a root block, and the rest blocks have inputs. For PBN 1 , assume its probabilistic transition matrix is ⎡ 1 ⎤ 1 1 l11 l12 . . . l12 n1 1 1 1 ⎢ l21 ⎥ l22 . . . l22 n1 ⎢ ⎥ L1 = ⎢ ⎥. .. ⎣ ⎦ . l21n1 1 l21n1 2 . . . l21n1 2n1 rank(L 1 − I2n1 ) = 2n1 − 1, and the right eigenvectors corresponding to eigenvalue 1 with l1 norm equal to 1 is S1 = [s11 , s21 , . . . , s21n1 ]T , where 0 ≤ si1 ≤ 1. Based on Definition 3.2, we can see that the steady-state distribution of PBN 1 is S1 .
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3001
For system 2 , which is a probabilistic Boolean control network that contains m2 inputs and n2 nodes, assume its probabilistic transition matrix is L2 , i.e., E (x 2 (t + 1)) = L 2 E (u2 (t )x 2 (t )), 2 2 where x 2 (t ) = nj=1 x2 j (t ), u2 (t ) = m k=1 u2k (t ). blocks as,
L2 can be partitioned to
2m2 equal size
L 2 = [L12 , L22 , . . . , L22m2 ]. Let
⎡
2 li11 2 ⎢ li21 ⎢ Li2 = ⎢ ⎣ li22 n2 1
2 li12 2 li22
... ... .. .
li22 n2 2
...
⎤ 2 li12 n2 2 ⎥ li22 n2 ⎥ ⎥. ⎦ li22 n2 2n2 2
Extend L2 to a higher dimension matrix L , such that 2
E (x 2 (t + 1)) = L E (x 1 (t )x 2 (t )), 2
n +n2
1 where x 1 (t ) = ni=1 x1i (t ). Then L ∈ R2 2 ×2 1
2
2
2
n
can be partitioned to 2n1 equal size blocks
2
L = [L 1 , L 2 , . . . , L 2 n1 ]. Next we will devote to finding the steady-state distributions of PBNs 1 × 2 . The prob2 abilistic transition matrix of 1 × 2 can be derived as (L 1 (I2n1 1T2n2 )) ∗ L . Assume PBN 1 × 2 has a steady-state distribution and is written as S2 = [s12 , s22 , . . . , s22n1 +n2 ]T . According to Definition 3.2, it is obvious to see that S2 is the right eigenvector of (L 1 (I2n1 2 1T2n2 )) ∗ L corresponding to eigenvalue 1. Partition S2 to 2n1 blocks, T
T
T
[S 21 , S 22 , . . . , S 22n1 ]T , then every block contains ⎧ 1 2 1 2 ⎪ l11 L 1 S 21 +l12 L 2 S 22 ⎪ ⎪ ⎪ 2 1 2 ⎨ l 1 L S 21 +l22 L 2 S 22 21 1 ⎪ ⎪ ⎪ ⎪ ⎩1 2 1 2 l2n1 1 L 1 S 21 +l12 L 2 S 22
2n2 elements, the sum of the elements in S 2i is si1 , and 2
+... +...
1 +l12 n1 L 2n1 S 22n1 2 1 +l22n1 L 2n1 S 22n1 .. .
+...
+l21n1 2n1 L 2n1 S 22n1
2
= S 21 , = S 22 ,
(10)
= S 22n1 .
Adding both sides of Eq. (10), we can get 21 n
2 L 1 S 21
+
2 L 2 S 22
+ ··· +
2 L 2n1 S 22n1
=
S 2i .
i=1
Moreover, in
2
L , the blocks corresponding to u2 = δ2jm2 are all equal to L 2j , hence, the 2
number of blocks in L that equal L 2j is 2n1 −m1 . Assume the set of states x1 corresponding to u2 = δ2i m2 is δ2n1 {i1 , i2 , . . . , i2n1 −m2 }, the set of states x1 corresponding to u2 = δ2jm2 is δ2n1 { j1 , j2 , . . . , j2n1 −m2 }, then it is obvious to see Xi2 = S 2i1 + S 2i2 + · · · + S 2i2n1 −m2 ,
3002
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
and L12 X12 + L22 X22 + · · · + L22m2 X22m2 = X12 + X22 + · · · + X22m2 .
(11)
Xi2
Solving Eq. (11), the least square solution of can be obtained. Adding the 2n1 −m2 equations that correspond to Xi2 in (10), we have ⎧ 2 2 LY +L22Y122 +... +L22m2 Y122 m2 = X12 , ⎪ ⎪ ⎪ 12 112 2 2 2 2 ⎨ L1 Y21 +L2 Y22 +... +L2m2 Y22m2 = X22 , .. ⎪ . ⎪ ⎪ ⎩ 2 2 L1 Y2m2 1 +L22Y22m2 2 + . . . +L22m2 Y22m2 2m2 = X22m2 .
(12)
Yi2j = (li11 j1 + li12 j1 + · · · + l21n1 −m1 j1 )S 2 j1 + (li11 j2 + · · · + l21n1 −m1 j2 )S 2 j2 + · · · + (li11 j n1 −m1 + li12 j n1 −m1 + · · · + l21n1 −m1 j n 2
2 1 −m1
2
)S 22n1 −m1 ,
(13)
Based on Eq. (12), the least square solution of Yi2j can be gotten. Then the least square solution of S 2i can be achieved, i.e., the steady-state distribution of 1 × 2 is obtained. For system e , 2 < e < μ, assume its probabilistic transition matrix is Le , then we have E (x e (t + 1)) = L e E (ue (t )x e (t )). Partition Le to 2me equal size blocks, L e = [L1e , L2e , . . . , L2i me ], where ⎡
l ej11 ⎢ l ej21 ⎢ L ej = ⎢ ⎣ l ej2ne 1
... ... .. . ...
l ej12 l ej22 l ej2ne 2
⎤ l ij12ne l ij22ne ⎥ ⎥ ⎥, ⎦ l ej2ne 2ne
Assume the steady-state distribution of 1 × 2 × · · · × e−1 is Se−1 = e e T e 1 [s1e−1 , s2e−1 , . . . , se−1 Let E (x (t + 1)) = L E (x (t )x 2 (t ) . . . x e (t )), then L can be
e−1 ] . n partitioned to 2 e
L =
2e−1i=1 i=1 ni
i
blocks,
e e e −1 n ]. [L 1 , L 2 , . . . , L 2 ei=1 i
If PBN 1 × 2 × · · · × e has a steady-state distribution, and it is written as Se = [s1e , s2e , . . . , s2e ei=1 ni ]T , partition it to 2 Se =
e−1 i=1
ni
equal size blocks,
T T T T −1 n ] . [S e1 , S e2 , . . . , S e2 ei=1 i
Similar to the method of calculating steady-state distributions of PBN 1 × 2 , the following equations can be obtained.
e−1
e L 1 S e1
+
e L 2 S e2
e
−1 n S + · · · + L 2 ei=1 i
e2
e−1 n i=1 i
=
n i=1 i 2
i=1
S ei ,
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3003
e
and the sum of elements in S ei is se−1i . In matrix L , the blocks corresponding to ue = δ2jme e is L ej , hence, in L the number of blocks that equal to L ej is 2n1 +n2 +···+ne−1 −me . Therefore, the above equation can be converted to L1e X1e + L2e X2e + · · · + L2eme x2eme = X1e + X2e + · · · + X2eme ,
(14) Xie
based on the above equation, the least square solution of can be achieved. Moreover, assume the set of states in x 1 x 2 . . . x e−1 that corresponds to ue = δ2i me is δ
2
e−1 n k=1 k
{i1 , i2 , . . . , i
2
e−1 n −me k=1 k
},
the set of states in x 1 x 2 . . . x e−1 that corresponds to ue = δ2jme is δ
2
e−1 n k=1 k
{ j1 , j2 , . . . , j
2
e−1 n −me k=1 k
},
therefore, Xie = S ei1 + S ei2 + · · · + S ei e−1 n −m . 2
k=1 k
e
e−1
Adding the 2 i=1 ni −me equations that correspond to Xie , we have ⎧ e e LY +L2eY12e +... +L2eme Y12e me = X1e ⎪ ⎪ ⎪ 1e 11e e e 2 e ⎨ L1Y21 +L2Y22 +... +L2me Y22me = X2e .. ⎪ . ⎪ ⎪ ⎩ e e L1Y2me 1 +L2eY2eme 2 + . . . +L2eme Y2eme 2me = X2em2e .
(15)
i Assume the set of xα (2 < α < e − 1) that corresponds to ue = δ2i me is δ2n1 Pα1 , the set of j j α e n x (2 < α < e − 1) that corresponds to u = δ2me is δ2 1 Pα2 . Let i Pα1 = {iα1 , iα2 , . . . , iατiα }, j Pα2 = { jα1 , jα2 , . . . , jατ jα },
then it is obvious that τiα = τ jα . Hence, Yij can be expressed as follows: i i Yiej = L 1 (P11 , j11 )Yiej11 + L 1 (P11 , j12 )Yiej22 e τ i +··· + L 1 (P11 , j1τ j1 )Yi j1 j1 .
(16) k1
j1k
When x 1 = δ2n 1 , assume the corresponding input to 2 is u2 = δ21m2 , hence, i i Yiej1k1 = Lk21 (P21 , j21 )Yiej21 + Lk21 (P21 , j22 )Yiej22 1 1 e τ i +··· + Lk21 (P21 , j2τ j2 )Yi j2 j2 .
(17)
1
j1k
j2k
k2
If we assume the input of 3 that corresponds to x 1 = δ2n11 , x 2 = δ2n22 is u3 = δ22m3 , therefore, i i Yiej2k2 = Lk32 (P31 , j31 )Yiej31 + Lk32 (P31 , j32 )Yiej32 2 2 e τ i +··· + Lk32 (P31 , j3τ j3 )Yi j3 j3 . (18) 2
3004
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Similar to the above analysis, assume the input to e−1 that corresponds to x 1 = δ2j1n1k1 , k e−2
j(e−2)k
j2k
x 2 = δ2n22 ,..., x e−2 = δ2ne−2 e−2 is δ2em−2e−1 , then the following result can be obtained, i e 1 i −2 Yiej(eke−2) = Lkee−1 Lkee−1 −2 (P(e−1)1 , j (e−1)1 )Yi j(e−1) + −2 (P(e−1)1 , j (e−1)2 ) · e−2 e−2 e τ j(e−1) i Yiej(e2 −1) + · · · + Lkee−1 −2 (P(e−1)1 , j (e−1)τ j(e−1) )Yi j(e−1) . Let
t
θ=1
e−2
(19)
t
nθ = N , thus
Yiej(es −1) = S e( e−2 ( j γ =1
γ kγ
−1)·2N
e−1 −N γ
)
.
(20)
Based on Eqs. (16)–(20), the sum of elements in Yiej can be obtained. By solving Eq. (15), the least square solution of Yiej can be gotten. In a similar manner, the sum of the elements in Yiej1k1 can be gained, and the least square solution of Yiej1k1 can be received by solving Eq. (16). And so on, the least square solution of Yiej(es −1) which meets condition that the sum e−1 of elements in Yiej(es −1) is s
can be achieved by solving equation. e−2 N e−1 −N γ γ =1 ( jγ kγ
−1)·2
Through the above analysis, an algorithm to calculate the steady-state distributions of PBN can be achieved. Algorithm. Step 1: If rank(L 1 − I2n1 ) 2n1 − 1, calculate the nonnegative right eigenvectors that corresponds to eigenvalue 1 with l1 norm equal to 1, i.e., the steady-state distributions of system 1 . For each eigenvector obtained in this step, proceed the following process. Else, the PBN has no steady-state distribution. Step 2: Based on the steady-state distributions obtained in Step 1, solve Eqs. (11)–(13), we can obtain the steady-state distributions of PBN 1 × 2 . If it exists, go to next step. Else, the PBN has no steady-state distribution. Repeat the above process. Step μ − 1: On the basis of the steady-state distributions obtained in Step μ − 2, and the above method, we can obtain the steady-state distributions of PBN 1 × 2 × · · · × μ−1 . If it exists, go to next step. Else, the PBN has no steady-state distribution. Step μ: Based on the steady-state distributions obtained in Step μ − 1, as well as the method presented above, we can obtain the steady-state distributions of PBN 1 × 2 × · · · × μ . If it exists, the steady-state distribution of the whole network is achieved. Else, the PBN has no steady-state distribution. 4. Numerical example Consider the following PBN, which contains 11 nodes x 1, x 2, . . . , x 11, its network aggregation graph is shown in Fig. 2. This graph only presents the relation between blocks. Let χ1 = {x 1, x 2, x 3, x 4}, χ2 = {x 5, x 6, x 7}, χ3 = {x 8, x 9, x 10, x 11}, and calculate its steady-state distribution. The probabilistic transition matrices of 1 , 2 and 3 can be found in [31]. For PBN 1 , with nodes set χ 1 , its steady-state distribution S1 can be computed as S1 = [0.0677, 0.0678, 0.0639, 0.0766, 0.0600, 0.0616, 0.0701, 0.0502, 0.0502, 0.0548, 0.0650, 0.0696, 0.0474, 0.0841, 0.0504, 0.0607]T . For system 2 , its nodes set is χ 2 , input is x4 , and its probabilistic transition matrix is L 2 = [L12 , L22 ].
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3005
Fig. 2. Example of aggregation of a network with 11 nodes.
Through Eq. (11), we have L12 X12
+ L22 X22 = X12 + X22 ,
8 1 1 2 and the sum of elements in X12 is 8j=1 s2· j−1 , the sum of elements in X2 is j=1 s2· j , hence, the least square solutions of X12 and X22 can be gained. X12 = [0.0775, 0.0849, 0.0834, 0.0370, 0.0398, 0.0206, 0.0643, 0.0673]T , X22 = [0.0705, 0.0419, 0.0612, 0.0958, 0.0667, 0.0799, 0.0609, 0.0484]T . Through Eq. (12), we have 2 2 L1 Y11 + L22Y122 = X12 , L12Y212 + L22Y222 = X22 , By solving Eq. (21), least square solutions Yi2j , i, j ∈ {1, 2} are obtained, where Y112 = [0.0288, 0.0364, 0.0286, 0.0252, 0.0234, 0.0231, 0.0266, 0.0258]T , Y122 = [0.0321, 0.0343, 0.0318, 0.0326, 0.0289, 0.0401, 0.0291, 0.0282]T , Y212 = [0.0306, 0.0278, 0.0303, 0.0308, 0.0314, 0.0321, 0.0306, 0.0304]T , Y222 = [0.0354, 0.0341, 0.0358, 0.0351, 0.0358, 0.0327, 0.0357, 0.0363]T . Let S12 = [S 21 , S 23 , S 25 , S 27 , S 29 , S 211 , S 213 , S 215 ]T , S22 = [S 22 , S 24 , S 26 , S 28 , S 210 , S 212 , S 214 , S 216 ]T , based on Eq. (13), it is obvious to see that Y112 = [0.3321, 0.5721, 0.3956, 0.5011, 0.6023, 0.4197, 0.5425, 0.4940]S12 , Y212 = [0.6679, 0.4279, 0.6044, 0.4989, 0.3797, 0.5803, 0.4575, 0.5060]S12 ,
(21)
3006
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Y122 = [0.4971, 0.5790, 0.3935, 0.5133, 0.5012, 0.4737, 0.3237, 0.5290]S22 , Y222 = [0.5029, 0.4210, 0.6065, 0.4867, 0.4988, 0.5263, 0.6763, 0.4710]S22 , then the steady-state distribution S2 of PBN 1 × 2 can be gotten via the above equation. S2 = [0.0085, 0.0083, 0.0076, 0.0090, 0.0066, 0.0083, 0.0062, 0.0065, 0.0086, 0.0097, 0.0078, 0.0064, 0.0070, 0.0088, 0.0106, 0.0077, 0.0082, 0.0096, 0.0077, 0.0098, 0.0082, 0.0085, 0.0073, 0.0073, 0.0087, 0.0100, 0.0077, 0.0066, 0.0071, 0.0089, 0.0104, 0.0079, 0.0084, 0.0082, 0.0075, 0.0089, 0.0065, 0.0082, 0.0061, 0.0064, 0.0086, 0.0096, 0.0079, 0.0064, 0.0070, 0.0088, 0.0107, 0.0077, 0.0084, 0.0075, 0.0073, 0.0084, 0.0057, 0.0079, 0.0055, 0.0060, 0.0086, 0.0098, 0.0078, 0.0064, 0.0070, 0.0088, 0.0106, 0.0077, 0.0084, 0.0072, 0.0073, 0.0082, 0.0054, 0.0078, 0.0052, 0.0058, 0.0081, 0.0090, 0.0076, 0.0059, 0.0065, 0.0084, 0.0106, 0.0072, 0.0086, 0.0072, 0.0074, 0.0082, 0.0053, 0.0079, 0.0053, 0.0058, 0.0095, 0.0111, 0.0080, 0.0074, 0.0079, 0.0095, 0.0103, 0.0088, 0.0084, 0.0078, 0.0074, 0.0086, 0.0061, 0.0080, 0.0058, 0.0062, 0.0081, 0.0090, 0.0076, 0.0059, 0.0065, 0.0084, 0.0106, 0.0072, 0.0083, 0.0076, 0.0073, 0.0084, 0.0059, 0.0079, 0.0056, 0.0060, 0.0081, 0.0089, 0.0076, 0.0058, 0.0064, 0.0084, 0.0107, 0.0071]T . For system 3 , its nodes set is χ 3 , input is x6 , and the probabilistic transition matrix is L 3 = [L13 , L23 ]. Through the above analysis, we have L13 X13 + L23 X23 = X13 + X23 ,
(22)
and the sum of elements in Xi3 is the sum of elements in S2 that corresponds to x6 = δ2i . The least square solutions of X13 , X23 can be received. X13 = [0.0297, 0.0290, 0.0282, 0.0312, 0.0362, 0.0223, 0.0323, 0.0320, 0.0371, 0.0341, 0.0347, 0.0341, 0.0338, 0.0376, 0.0345, 0.0289]T , X23 = [0.0278, 0.0272, 0.0258, 0.0285, 0.0330, 0.0210, 0.0296, 0.0290, 0.0343, 0.0320, 0.0331, 0.0319, 0.0334, 0.0365, 0.0336, 0.0279]T .
According to Eq. (15), next we will solve the following equation: L13Y113 + L23Y123 = X13 , L13Y213 + L23Y223 = X23 ,
(23)
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3007
Its least square solution is shown as follows: Y113 = [0.0129, 0.0148, 0.0147, 0.0170, 0.0159, 0.0092, 0.0134, 0.0149, 0.0156, 0.0174, 0.0164, 0.0164, 0.0195, 0.0166, 0.0169, 0.0174]T , Y123 = [0.0102, 0.0169, 0.0122, 0.0181, 0.0125, 0.0143, 0.0162, 0.0127, 0.0222, 0.0159, 0.0174, 0.0163, 0.0122, 0.0169, 0.0176, 0.0167]T , Y213 = [0.0199, 0.0162, 0.0156, 0.0157, 0.0160, 0.0174, 0.0165, 0.0154, 0.0185, 0.0182, 0.0154, 0.0186, 0.0139, 0.0190, 0.0176, 0.0128]T , Y223 = [0.0153, 0.0096, 0.0173, 0.0107, 0.0197, 0.0118, 0.0120, 0.0133, 0.0155, 0.0151, 0.0139, 0.0170, 0.0183, 0.0158, 0.0174, 0.0137]T . Because the following process is tedious, there is no need to enumerate all the discussions here. We only write the finial result which is achieved by the method presented in this paper, i.e., the steady-state distribution S3 of PBN 1 × 2 × 3 . 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 S3 = [γ11 , γ21 , φ11 , φ21 , γ31 , γ41 , φ31 , φ41 , γ51 , γ61 , φ51 , φ61 , γ71 , γ81 , φ71 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ81 , γ12 , γ22 , φ12 , φ22 , γ32 , γ42 , φ32 , φ42 , γ52 , γ62 , φ52 , φ62 , γ72 , γ82 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ72 , φ82 , γ13 , γ23 , φ13 , φ23 , γ33 , γ43 , φ33 , φ43 , γ53 , γ63 , φ53 , φ63 , γ73 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ83 , φ73 , φ83 , γ14 , γ24 , φ14 , φ24 , γ34 , γ44 , φ34 , φ44 , γ54 , γ64 , φ54 , φ64 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ74 , γ84 , φ74 , φ84 , γ15 , γ25 , φ15 , φ25 , γ35 , γ45 , φ35 , φ45 , γ55 , γ65 , φ55 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ65 , γ75 , γ85 , φ75 , φ85 , γ16 , γ26 , φ16 , φ26 , γ36 , γ46 , φ36 , φ46 , γ56 , γ66 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ56 , φ66 , γ76 , γ86 , φ76 , φ86 , γ17 , γ27 , φ17 , φ27 , γ37 , γ47 , φ37 , φ47 , γ57 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ67 , φ57 , φ67 , γ77 , γ87 , φ77 , φ87 , γ18 , γ28 , φ18 , φ28 , γ38 , γ48 , φ38 , φ48 , 3 3 3 3 3 3 3 3 T γ58 , γ68 , φ58 , φ68 , γ78 , γ88 , φ78 , φ88 ] .
In order to save paper, all the elements of S3 that can be found in [31] are omitted. 5. Conclusion In this paper, we have considered the steady-state distributions of PBNs with special structure. To solve this problem, the network has been partitioned into several blocks by cascading aggregation. Then finding steady-state distributions of the whole networks has been converted to computing several corresponding equations. Two necessary and sufficient conditions for the existence of the steady-state distributions for PBNs have been obtained. An algorithm has been presented to calculate the steady-sate distributions for the whole networks. In fact, our method can deal with large-scale PBNs theoretically, but due to the complexity of computation, we give an example with 11 nodes. Acknowledgment This work was supported by National Natural Science Foundation under Grants 61773371 and 61877036, and China Postdoctoral Science Foundation under Grant 2016M-602143.
3008
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
References [1] S. Kim, S. Imoto, S. Miyano, Dynamic Bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data, Biosystems 75 (2004) 57–65. [2] S. Kauffman, Metabolic stability and epigenesis in randomly constructed genetic nets, J. Theor. Biol. 244 (4) (1969) 670–679. [3] L. Steggles, R. Banks, A. Wipat, Modelling and analysing genetic networks: from Boolean networks to petri nets 4210 (2006) 127–141. [4] W. Ching, M. Michael, S. Eric, T. Akutsu, On construction of stochastic genetic networks based on gene expression sequences, Int. J. Neural Syst. 15 (4) (2005) 297–310. [5] S. Zhang, W. Ching, N. Tsing, H. Leung, D. Guo, A new multiple regression approach for the construction, Artif. Intell. Med. 48 (2010) 153–160. [6] I. Shmulevich, E.R. Dougherty, S. Kim, W. Zhang, Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks, Bioinformatics 18 (2002) 261–274. [7] D. Cheng, H. Qi, Z. Li, Analysis and Control of Boolean Networks: A Semi-tensor Product Approach, Springer, 2011. [8] M. Meng, L. Liu, G. Feng, Stability and L1 gain analysis of Boolean networks with Markovian jump parameters, IEEE Trans. Autom. Control 62 (8) (2017) 4222–4228. [9] H. Li, Y. Wang, Z. Wu, Further results on feedback stabilization control design of Boolean control networks, Automatica 83 (2017) 303–308. [10] Y. Yu, J. Feng, J. Pan, Block decoupling of Boolean control networks, IEEE Trans. Autom. Control (2018), doi:10.1109/TAC.2018.2880411. [11] S. Wang, J. Feng, J. Zhao, J. Xia, Controllability decomposition of dynamic-algebraic Boolean control networks, Int. J. Control (2018), doi:10.1080/00207179.2018.1527040. [12] Y. Zhao, D. Cheng, On controllability and stabilizability of probabilistic Boolean control networks, Sci. China 57 (1) (2014) 1–14. [13] Y. Liu, H. Chen, J. Lu, B. Wu, Controllability of probabilistic Boolean control networks based on transition probability matrices, Automatica 52 (2015) 340–345. [14] K. Zhang, L. Zhang, Controllability of probabilistic Boolean control networks with time-variant delays in states, in: Proceedings of the 52th IEEE Conference on Decision and Control, 2013, pp. 7211–7216. [15] M. Meng, J. Lam, J. Feng, K.C. Cheung, Stability and stabilization of Boolean networks with stochastic delays, IEEE Trans. Autom. Control (2018), doi:10.1109/TAC.2018.2835366. [16] Z. Li, J. Song, J. Yang, Partial stability of probabilistic Boolean network, in: Proceedings of the 26th Chinese Control and Decision Conference, 2014, pp. 1952–1956. [17] J. Zhao, Z. Liu, Observability of probabilistic Boolean networks, in: Proceedings of the 34th Chinese Control Conference, 2015, pp. 183–186. [18] H. Chen, J. Liang, J. Lu, J. Qiu, Synchronization for the realization-dependent probabilistic Boolean networks, IEEE Trans. Neural Netw. Learn. Syst. 29 (4) (2018) 819–831. [19] J. Lu, J. Zhong, L. Li, D.W.C. Ho, J. Cao, Synchronization analysis of master-slave probabilistic Boolean networks, Sci. Rep. 5 (2015) 13437. [20] Z. Liu, Y. Wang, Optimal finite-horizon control problem of context-sensitive probabilistic Boolean networks with perturbation, in: Proceedings of the 31th Chinese Control Conference, 2012, pp. 140–145. [21] M. Yang, R. Li, T. Chu, A new method and application for controlling steady-state probability distributions of probabilistic Boolean networks, in: Proceedings of the 2014 IEEE Congress on Evolutionary Computation, 2014, pp. 1490–1495. [22] X. Cheng, Y. Qiu, W. Hou, W.K. Ching, A semi-tensor product approach for probabilistic Boolean networks, in: Proceedings of the 8th International Conference on System Biology, 2014, pp. 85–90. [23] R. Pal, Context-sensitive probabilistic Boolean networks: steady-state properties, reduction, and steady-state approximation, IEEE Trans. Signal Process. 58 (2) (2010) 879–890. [24] I. Shmulevich, I. Gluhovsky, R.F. Hashimoto, E.R. Dougherty, W. Zhang, Steady-state analysis of genetic regulatory networks modelled by probabilistic Boolean networks, Compar. Funct. Genom. 4 (2003) 601–608. [25] W.K. Ching, S. Zhang, M.K. Ng, T. Akutsu, An approximation method for solving the steady-state probability distribution of probabilistic Boolean networks, Bioinformatics 23 (12) (2007) 1511–1518. [26] I. Apostolopoulou, D. Marculescu, Tractable learning and inference for large-scale probabilistic Boolean networks, IEEE Trans. Neural Netw. Learn. Syst. (2019), doi:10.1109/TNNLS.2018.2886207.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3009
[27] Y. Zhao, J. Kim, M. Filippone, Aggregation algorithm towards large-scale Boolean network analysis, in: Proceedings of the IEEE Transactions on Automatic Control, 58, 2013, pp. 1976–1985. [28] Y. Zhao, B.K. Ghosh, D. Cheng, Control of large-scale Boolean networks via network aggregation, IEEE Trans. Neural Netw. Learn. Syst. 27 (7) (2016) 1527–1536. [29] K. Zhang, Observability and reconstructibility of large-scale Boolean control networks via network aggregations, (2017) arXiv preprint, arXiv:1704.03231. [30] J.G. Kemeny, J.L. Snell, Finite Markov Chains, Springer, 1983. [31] J. Pan, Control and Applications of Boolean Networks via Algebraic Method, 2018 Ph.D thesis.
Journal of the Franklin Institute 356 (2019) 2994–3009 www.elsevier.com/locate/jfranklin
Steady-state analysis of probabilistic Boolean networks Jinfeng Pan a, Jun-e Feng a,∗, Min Meng b b School
a School of Mathematics, Shandong University, Jinan 250100, China of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798, Singapore
Received 20 April 2018; received in revised form 18 January 2019; accepted 27 January 2019 Available online 1 February 2019
Abstract This paper investigates steady-state distributions of probabilistic Boolean networks via cascading aggregation. Under this approach, the problem is converted to computing least square solutions to several corresponding equations. Two necessary and sufficient conditions for the existence of the steady-state distributions for probabilistic Boolean networks are given firstly. Secondly, an algorithm for finding the steady-state distributions of probabilistic probabilistic Boolean networks is given. Finally, a numerical example is given to show the effectiveness of the proposed method. © 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
1. Introduction Developing models to analyse gene regulatory networks is one of the aims in systems biology. Numbers of mathematical models have been proposed to study gene regulatory networks, such as Bayesian networks [1], Boolean networks (BNs) [2], Petri nets [3], probabilistic Boolean networks (PBNs) [4], linear regression model [5]. Among them, BNs and PBNs, which were originally introduced by Kauffman [2], are some of the most attractive models. One BN is defined by a set of nodes (genes) {x1 , x2 , . . . , xn }, xi ∈ {0, 1}, and a set of Boolean functions { f1 , f2 , . . . , fn }, with the update rules as xi (t + 1) = fi (x j (t ), j ∈ Ni ). Hence, for a ∗
Corresponding author. E-mail addresses: [email protected] (J. Pan), [email protected] (J.-e. Feng), mengminmath@ gmail.com (M. Meng). https://doi.org/10.1016/j.jfranklin.2019.01.039 0016-0032/© 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2995
BN with n nodes, there are 2n states. On the basis of update rules, the state transition graph of a BN can be uniquely determined, and in this graph, the out-degree of any state is 1. In addition, due to the uncertainty of gene regulation and errors of experimental noise, it is appropriate that one update rule of a node in a BN is randomly chosen at each time among the candidates of Boolean functions, which results in the appearance of PBNs [6]. As an extension of BNs, PBNs have been widely studied in recent years. A PBN can be viewed as a collection of BNs. For a PBN with n nodes {x1 , x2 , . . . , xn }, the update rule f(i) of xi is chosen from { fi1 , fi2 , . . . , fiq(i) } in probability. Then for the whole network, there are ni=1 q(i) BNs. In this paper, the vector-valued function fl ( fl(l ) , fl(2) , . . . , fl(n) ) is used to determine a context of the PBN. For each gene, at each time point, a BN is randomly chosen to determine the network’s transition. A PBN can be regarded as a homogeneous Markov chain. Therefore, some results of Markov chains can also be used to analyse PBNs. Recently, a new matrix product, semi-tensor product of matrices, is proposed by Cheng and his co-workers [7]. Based on it, BNs can be converted to algebraic forms, under which many interesting results were obtained [8–11]. Nowadays, the semi-tensor product is applied to solve problems for PBNs as well, such as the controllability [12–14], stability [15,16], observability [17], synchronization [18,19], optimal problems [20], steady-state distributions [21,22]. When considering the long-run behavior of a PBN, finding its steady-state distributions is an important aspect of analysing biological systems. This steady-state distribution can provide a first-order statistical information for PBNs. On this basis, the network and the impact of different genes can be understood. In addition, one can figure out how to control some genes in the network so that the whole system can evolve to the target state or the expected steadystate probability distribution. Then we can develop and research therapeutic gene intervention or gene control policy [6]. If a PBN has steady-state distribution, then in the long run, the probability that the chain is in a given state is stationary, independent of the initial state. But if the PBN has no steady-state distribution, such question need not be answered. The steady-state distributions can be obtained based on transition matrices. Up to now, several references were devoted to finding the steady-state distributions of PBNs [23–25]. Context-reduction mapping and an approximation method to calculate the steady-state distributions of context-sensitive PBNs were proposed in [23]. Monte Carlo methods [24] were used to calculate the steadystate distributions of PBNs. In [25], an approximation method for computing the steady-state distributions of PBNs was presented by neglecting some BNs with very small probabilities during the construction of the transition probabilistic matrices. In [26], the long run behavior of large-scale PBN was analyzed based on its stochastic conjunctive normal form network. In [27,28], an algorithm was proposed to find the attractors of large-scale Boolean networks by network aggregation, which reduces the computation significantly. Transition probabilistic matrix-based algorithm to find the steady-state distributions of PBNs is complicated even if the number of states is not large. Hence, in this paper, the problem to find the steady-state distributions of PBNs with special structures is also considered by using the cascading aggregation. Using semi-tensor product, the probabilistic transition matrix of each block can be obtained. The main contributions of this paper are concluded as follows: For PBNs, two equivalent conditions for the existence of the steady-state distribution are proposed on the basis of the probabilistic transition matrices. The network aggregation is firstly used to find steady-state distributions of PBNs. Based on cascading aggregation, the problem can be solved by calculating the least square solutions to several corresponding equations. By the method presented in this paper, the computational complexity is reduced.
2996
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
The rest of this paper is organized as follows. Section 2 introduces some fundamental definitions and notations. The definition and properties of Markov chains, and two necessary and sufficient conditions for the existence of the steady-state distributions for PBNs, as well as a method to compute steady-state distributions for PBNs, are derived in Section 3. In Section 4, a numerical example is given to illustrate the effectiveness of the proposed method. Section 5 makes a brief conclusion. 2. Preliminaries and problem formulation In this section, we present some properties about semi-tensor product of matrices that will be used in this paper, and the definition of network aggregation of PBNs. 2.1. Probabilistic Boolean networks A PBN can be regarded as an extension of BN with a probabilistic setting. Consider the following PBN:
⎧ x1 (t + 1) = f (1) ([x j (t )], j ∈ N1 ), ⎪ ⎪ ⎪ ⎨x2 (t + 1) = f (2) ([x j (t )], j ∈ N2 ), .. ⎪ . ⎪ ⎪ ⎩ xn (t + 1) = f (n) ([x j (t )], j ∈ Nn ),
(1)
where [x1 , x2 , . . . , xn ]T ∈ D n is the state vector, D = {0, 1}. Ni is the index set of xi . Assume that there are q(i) candidates of f(i) as fl(i) ([x j (t )], j ∈ Nl(i) ), l = 1, 2, . . . , q(i),
(2)
where Nl(i) is the index set corresponding to fl(i) , and the probability of choosing fl(i) as the update function is cl(i) , i.e., cl(i) = Prob( f (i) = fl(i) ), l = 1, 2, . . . , q(i).
(3)
The following relation q(i)
cl(i) = 1
(4)
l=1
must be satisfied. Ni , i = 1, 2, . . . , n, can be obtained as Ni :=
q(i)
Nl(i) .
(5)
l=1
Due to the stochasticity in genetic networks, the next state of gene i is determined by q(i) (i) (i) Boolean functions f1(i) , f2(i) , . . . , fq(i) , with probabilities c1(i) , c2(i) , . . . , cq(i) , respectively. For a n PBN of n genes, there are i=1 q(i) possible Boolean networks, each of which is a possible realization of the genetic network and can be considered as a context of the PBN. For the (1) (2) (n) jth context, assume the network function is given by f j = ( f j(1) , f j(2) , . . . , f j(n) ), then the n (i) probability that the jth context is selected is c j = i=1 c j(i) .
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2997
Example 2.1. Consider a PBN with three nodes. The update rules are given as follows: f = 1
f1(1) = x1 (t ) ∨ x2 (t ), c1(1) = 0.8, f2(1) = x1 (t ), c2(1) = 0.2,
f 2 = f1(2) = x1 (t ), c1(2) = 1, (3) f = x3 (t ), c1(3) = 0.5, f 3 = 1(3) f2 = ¬x3 (t ), c2(3) = 0.5. Then we have q(1) = 2, q(2) = 1, q(3) = 2, N1(1) = {1, 2}, N2(1) = {1}, N1 = N1(1) ∪ N2(1) = {1, 2}, N2 = N1(2) = {1}, N3 = N1(3) ∪ N2(3) = {3}, and Prob( f 1 = f1(1) , f 2 = f1(2) , f 3 = f1(3) ) = 0.4. 2.2. Algebraic forms of Boolean functions For ease of presentation, we first give some notations, many of which are standard in literature [7]. Rm × n is the set of m × n real matrices. For a matrix M ∈ Rm × n , the (i, j)th entry is denoted by Mij . M > 0 means Mij > 0 for all i, j. coli (M) (rowi (M )) is the ith column (row) of matrix M. col(M) (row(M )) is the set
of columns (rows) of M. For a given vector x = (x1 , x2 , . . . , xn ), the l1 norm of x is xl1 = ni=1 |xi |. δni is the ith column of the identity matrix In . n = {δni |i = 1, 2, . . . , n}. Assume that L ∈ Lm×n , L = [δmi1 δmi2 . . . δmin ], and its shorthand form is written as L = δm [i1 i2 . . . in ]. Given a set S = {c1 , c2 , . . . , cl }, 1 ≤ ci ≤ 2n , ci = cj , δ2n S is used to represent the set {δ2c1n , δ2c2n , . . . , δ2cln }. Given a matrix L ∈ Rm × n , assume
l P = { p1 , p2 , . . . , pl }, pi ∈ {1, 2, . . . , m}, we denote L(P, s) to represent k=1 L pk s , s ∈ {1, 2, . . . , n}. Assume A ∈ Rm × r , B ∈ Rn × r , the Khatri–Rao product of A and B is defined as follows: A ∗ B = [col1 (A ) col1 (B), col2 (A ) col2 (B), . . . , colr (A ) colr (B)]. To use matrix expression of logic, we identify T = 1 ∼ δ21 , F = 0 ∼ δ22 , thus D ∼ 2 , where ∼ means the logical equivalence. Subsequently, every Boolean function f : D n → D, can be equivalently expressed as a mapping f : n2 → 2 . Lemma 2.1. ([7]) Consider a Boolean function f (x1 , x2 , . . . , xn ) with xi ∈ 2 . Denote x = x1 x2 · · · xn , then there exists a unique matrix M f ∈ L2×2n , called the structural matrix of f, such that f (x1 , x2 , . . . , xn ) = M f ni=1 xi = M f x.
(6)
2.3. Aggregation of PBNs In this section, we first introduce the network aggregation, which plays a vital role in studying the steady-state distributions of PBNs. Denote χ = {x1 , x2 , . . . , xn }, then χ can be partitioned into μ blocks χ = χ1 ∪ χ2 ∪ · · · ∪ χμ , where χ i = ∅, χ i ⊂ χ , and χi ∩ χ j = ∅ for i
= j. Let χi = {xi1 , xi2 , . . . , xini }, where ni is the number of state nodes in the ith block, and μ i=1 ni = n.
2998
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Fig. 1. Example of aggregation of a network with 9 nodes into 3 Boolean control networks.
If j ∈ Ni , we say that there is a directed edge from xj to xi . All the starting nodes of the incoming and outgoing edges of block χ i constitute sets Ui = {ui1 , . . . , uimi } and Yi = {yi1 , . . . , yipi }, respectively. Set Ui (Yi ) may be empty, when there are no incoming (outgoing) edges of block χ i . Let C be the set of starting nodes whose edges cut by the aggregation, then we have C = ∪μi=1 Ui = ∪μi=1 Yi . Therefore the subnetwork i with nodes in χ i and inputs in Ui is a probabilistic Boolean control network, which can be described as follows: i :
xi j (t + 1) = f (i j) (xi1 (t ), xi2 (t ), . . . , xini (t ), ui1 (t ), ui2 (t ), . . . , uimi (t )),
(7)
for i = 1, 2, . . . , μ and j = 1, 2, . . . , ni .cl(i j) = Prob( f (i j) = fl(i j) ), l = 1, 2, . . . , q(i j). Definition 2.2. [27] The aggregation is said to be cascading if (with possible order) the groups χ i = ∪ij=1 χ j , i = 1, 2, . . . , μ,
(8)
have zero in-degree. Definition 2.3. [27] The aggregation is said to be an acyclic aggregation, if its aggregation graph contains no cycle. Lemma 2.4. [29] An aggregation is a cascading aggregation if and only if it is an acyclic aggregation. The method to find the acyclic aggregation of a network has been given in literature [27]. In what follows, an example is given to show the aggregation of PBNs.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
2999
Example 2.2. Consider a PBN in Fig. 1, and the update rules are given as follows: ⎧ ⎪ x1 (t + 1) = x2 (t ), c1(1) = 1, ⎪ ⎪ ⎪ ⎪ x3 (t ) ∧ x7 (t ), c1(2) = 0.7, ⎪ ⎪ x (t + 1) = 2 ⎪ ⎪ x3 (t ), c2(2) = 0.3, ⎪ ⎪ ⎪ ⎪ ⎪ x3 (t + 1) = x1 (t ) ∨ x2 (t ), c1(3) = 1, ⎪ ⎪ ⎪ ⎪ x1 (t ) ∨ x5 (t ) ∨ x7 (t ), c1(4) = 0.5, ⎪ ⎪ x (t + 1) = ⎨ 4 x1 (t ), c2(4) = 0.5, ⎪ x5 (t + 1) = ¬x4 (t ), c1(5) = 1, ⎪ ⎪ ⎪ ⎪ x6 (t + 1) = x6 (t ), c1(6) = 1, ⎪ ⎪ ⎪ ⎪ x7 (t + 1) = x6 (t ), c1(7) = 1, ⎪ ⎪ ⎪ ⎪ ⎪ x (t ) ∨ x9 (t ), c1(8) = 0.6, ⎪ ⎪ x8 (t + 1) = 7 ⎪ ⎪ x7 (t ) ∧ x9 (t ), c2(8) = 0.4, ⎪ ⎪ ⎩ x9 (t + 1) = x6 (t ), c1(9) = 1. The node set χ = {x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 , x9 } can be partitioned into 3 blocks χ i, i = 1, 2, 3, where χ 1 = {x1 , x2 , x3 }, χ 2 = {x4 , x5 }, and χ 3 = {x6 , x7 , x8 , x9 }. Reorder the blocks, and define χ 1 = χ 3, χ 2 = χ 3 ∪ χ 1, χ 3 = χ . By Definition 2.2, the aggregation is cascading. 3. Steady-state distribution of PBNs In this section, the definitions of Markov chains and steady-state distribution as well as some properties are presented firstly. Then the method to calculate the steady-state distributions of PBNs under cascading aggregation is given in two parts. Part 1 provides two necessary and sufficient conditions for the existence of the steady-state distributions of PBNs. Part 2 gives a method of calculating the steady-state distributions of PBNs. Definition 3.1. [30] Assume that X = {Xn , n = 0, 1, . . .} is a discrete random process which is defined in (, F, P ), and its state space S is a countable set or finite set. If for every nonnegative integer n, and arbitrary state i0 , i1 , . . . , in+1 ∈ S, Prob(X0 = i0 , X1 = i1 , . . . , Xn = in ) > 0, and Prob(Xn+1 = in+1 |X0 = i0 , . . . , Xn = in ) = Prob(Xn+1 = in+1 |Xn = in ), then we say that X is a discrete Markov chain, here Prob(Xn+1 = in+1 |Xn = in ) has no relationship with n. For a given Markov chain, if its state space is finite, then its transition probability can be expressed by a probabilistic matrix P, which is called the probabilistic transition matrix, Pi j = Prob(Xn+1 = j | Xn = i).
(9)
If the Markov chain is homogeneous, then the k step probabilistic transition matrix is Pk . Definition 3.2. [30] Given a finite homogeneous Markov chain with state space S = {1, 2, . . . , M}, we say the state space has a steady-state distribution, if there exists a probability distribution π = (π1 , π2 , . . . , πM )
3000
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
such that, πj =
M
πi Pi j .
i=1
Based on Lemma 2.1, the jth context of (1) can be converted to x(t + 1) = L j x(t ). Define the probabilistic transition matrix of system (1) as follows: L = c1 L 1 + c2 L 2 + · · · + cQ L Q , where Q = ni=1 q(i), the probability of choosing the jth context is cj . From the construction of matrix L, it is obvious to see that (L) ji = Prob(x(t + 1) = δ2jn |x(t ) = δ2i n ), hence the probabilistic transition matrix of the Markov chain associated with Eq. (1) is LT . Based on the results on Markov chains in [30], it is obviously to have the following results. Theorem 3.3. Assume that LT is the transition probability matrix of Eq. (1), then (i) there exists at least one steady-state distribution if and only if the geometric multiplicity of eigenvalue 1 of L is greater than or equal to 1. (ii) there exists a unique steady-state distribution if and only if the geometric multiplicity of eigenvalue 1 of L is equal to 1. Corollary 3.4. Let L be the probabilistic transition matrix of a finite homogeneous Markov chain, (i) there exists at least one steady-state distribution if and only if rank(L − I2n ) 2n − 1. (ii) there exists a unique steady-state distribution if and only if rank(L − I2n ) = 2n − 1. 3.1. Steady-state distribution of PBNs In the above discussions, two conditions are derived to judge the existence of the steadystate distribution based on the probabilistic transition matrix. However, when the dimensions of the matrix become bigger, the methods are not applicable. In the following process, we propose a method to compute the steady-state distributions for PBNs under network aggregation. Assume PBN (1) can be partitioned to μ blocks 1 , 2 , . . . , μ by cascading aggregation, and the blocks are ordered according to Definition 2.2. Thus 1 is a root block, and the rest blocks have inputs. For PBN 1 , assume its probabilistic transition matrix is ⎡ 1 ⎤ 1 1 l11 l12 . . . l12 n1 1 1 1 ⎢ l21 ⎥ l22 . . . l22 n1 ⎢ ⎥ L1 = ⎢ ⎥. .. ⎣ ⎦ . l21n1 1 l21n1 2 . . . l21n1 2n1 rank(L 1 − I2n1 ) = 2n1 − 1, and the right eigenvectors corresponding to eigenvalue 1 with l1 norm equal to 1 is S1 = [s11 , s21 , . . . , s21n1 ]T , where 0 ≤ si1 ≤ 1. Based on Definition 3.2, we can see that the steady-state distribution of PBN 1 is S1 .
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3001
For system 2 , which is a probabilistic Boolean control network that contains m2 inputs and n2 nodes, assume its probabilistic transition matrix is L2 , i.e., E (x 2 (t + 1)) = L 2 E (u2 (t )x 2 (t )), 2 2 where x 2 (t ) = nj=1 x2 j (t ), u2 (t ) = m k=1 u2k (t ). blocks as,
L2 can be partitioned to
2m2 equal size
L 2 = [L12 , L22 , . . . , L22m2 ]. Let
⎡
2 li11 2 ⎢ li21 ⎢ Li2 = ⎢ ⎣ li22 n2 1
2 li12 2 li22
... ... .. .
li22 n2 2
...
⎤ 2 li12 n2 2 ⎥ li22 n2 ⎥ ⎥. ⎦ li22 n2 2n2 2
Extend L2 to a higher dimension matrix L , such that 2
E (x 2 (t + 1)) = L E (x 1 (t )x 2 (t )), 2
n +n2
1 where x 1 (t ) = ni=1 x1i (t ). Then L ∈ R2 2 ×2 1
2
2
2
n
can be partitioned to 2n1 equal size blocks
2
L = [L 1 , L 2 , . . . , L 2 n1 ]. Next we will devote to finding the steady-state distributions of PBNs 1 × 2 . The prob2 abilistic transition matrix of 1 × 2 can be derived as (L 1 (I2n1 1T2n2 )) ∗ L . Assume PBN 1 × 2 has a steady-state distribution and is written as S2 = [s12 , s22 , . . . , s22n1 +n2 ]T . According to Definition 3.2, it is obvious to see that S2 is the right eigenvector of (L 1 (I2n1 2 1T2n2 )) ∗ L corresponding to eigenvalue 1. Partition S2 to 2n1 blocks, T
T
T
[S 21 , S 22 , . . . , S 22n1 ]T , then every block contains ⎧ 1 2 1 2 ⎪ l11 L 1 S 21 +l12 L 2 S 22 ⎪ ⎪ ⎪ 2 1 2 ⎨ l 1 L S 21 +l22 L 2 S 22 21 1 ⎪ ⎪ ⎪ ⎪ ⎩1 2 1 2 l2n1 1 L 1 S 21 +l12 L 2 S 22
2n2 elements, the sum of the elements in S 2i is si1 , and 2
+... +...
1 +l12 n1 L 2n1 S 22n1 2 1 +l22n1 L 2n1 S 22n1 .. .
+...
+l21n1 2n1 L 2n1 S 22n1
2
= S 21 , = S 22 ,
(10)
= S 22n1 .
Adding both sides of Eq. (10), we can get 21 n
2 L 1 S 21
+
2 L 2 S 22
+ ··· +
2 L 2n1 S 22n1
=
S 2i .
i=1
Moreover, in
2
L , the blocks corresponding to u2 = δ2jm2 are all equal to L 2j , hence, the 2
number of blocks in L that equal L 2j is 2n1 −m1 . Assume the set of states x1 corresponding to u2 = δ2i m2 is δ2n1 {i1 , i2 , . . . , i2n1 −m2 }, the set of states x1 corresponding to u2 = δ2jm2 is δ2n1 { j1 , j2 , . . . , j2n1 −m2 }, then it is obvious to see Xi2 = S 2i1 + S 2i2 + · · · + S 2i2n1 −m2 ,
3002
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
and L12 X12 + L22 X22 + · · · + L22m2 X22m2 = X12 + X22 + · · · + X22m2 .
(11)
Xi2
Solving Eq. (11), the least square solution of can be obtained. Adding the 2n1 −m2 equations that correspond to Xi2 in (10), we have ⎧ 2 2 LY +L22Y122 +... +L22m2 Y122 m2 = X12 , ⎪ ⎪ ⎪ 12 112 2 2 2 2 ⎨ L1 Y21 +L2 Y22 +... +L2m2 Y22m2 = X22 , .. ⎪ . ⎪ ⎪ ⎩ 2 2 L1 Y2m2 1 +L22Y22m2 2 + . . . +L22m2 Y22m2 2m2 = X22m2 .
(12)
Yi2j = (li11 j1 + li12 j1 + · · · + l21n1 −m1 j1 )S 2 j1 + (li11 j2 + · · · + l21n1 −m1 j2 )S 2 j2 + · · · + (li11 j n1 −m1 + li12 j n1 −m1 + · · · + l21n1 −m1 j n 2
2 1 −m1
2
)S 22n1 −m1 ,
(13)
Based on Eq. (12), the least square solution of Yi2j can be gotten. Then the least square solution of S 2i can be achieved, i.e., the steady-state distribution of 1 × 2 is obtained. For system e , 2 < e < μ, assume its probabilistic transition matrix is Le , then we have E (x e (t + 1)) = L e E (ue (t )x e (t )). Partition Le to 2me equal size blocks, L e = [L1e , L2e , . . . , L2i me ], where ⎡
l ej11 ⎢ l ej21 ⎢ L ej = ⎢ ⎣ l ej2ne 1
... ... .. . ...
l ej12 l ej22 l ej2ne 2
⎤ l ij12ne l ij22ne ⎥ ⎥ ⎥, ⎦ l ej2ne 2ne
Assume the steady-state distribution of 1 × 2 × · · · × e−1 is Se−1 = e e T e 1 [s1e−1 , s2e−1 , . . . , se−1 Let E (x (t + 1)) = L E (x (t )x 2 (t ) . . . x e (t )), then L can be
e−1 ] . n partitioned to 2 e
L =
2e−1i=1 i=1 ni
i
blocks,
e e e −1 n ]. [L 1 , L 2 , . . . , L 2 ei=1 i
If PBN 1 × 2 × · · · × e has a steady-state distribution, and it is written as Se = [s1e , s2e , . . . , s2e ei=1 ni ]T , partition it to 2 Se =
e−1 i=1
ni
equal size blocks,
T T T T −1 n ] . [S e1 , S e2 , . . . , S e2 ei=1 i
Similar to the method of calculating steady-state distributions of PBN 1 × 2 , the following equations can be obtained.
e−1
e L 1 S e1
+
e L 2 S e2
e
−1 n S + · · · + L 2 ei=1 i
e2
e−1 n i=1 i
=
n i=1 i 2
i=1
S ei ,
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3003
e
and the sum of elements in S ei is se−1i . In matrix L , the blocks corresponding to ue = δ2jme e is L ej , hence, in L the number of blocks that equal to L ej is 2n1 +n2 +···+ne−1 −me . Therefore, the above equation can be converted to L1e X1e + L2e X2e + · · · + L2eme x2eme = X1e + X2e + · · · + X2eme ,
(14) Xie
based on the above equation, the least square solution of can be achieved. Moreover, assume the set of states in x 1 x 2 . . . x e−1 that corresponds to ue = δ2i me is δ
2
e−1 n k=1 k
{i1 , i2 , . . . , i
2
e−1 n −me k=1 k
},
the set of states in x 1 x 2 . . . x e−1 that corresponds to ue = δ2jme is δ
2
e−1 n k=1 k
{ j1 , j2 , . . . , j
2
e−1 n −me k=1 k
},
therefore, Xie = S ei1 + S ei2 + · · · + S ei e−1 n −m . 2
k=1 k
e
e−1
Adding the 2 i=1 ni −me equations that correspond to Xie , we have ⎧ e e LY +L2eY12e +... +L2eme Y12e me = X1e ⎪ ⎪ ⎪ 1e 11e e e 2 e ⎨ L1Y21 +L2Y22 +... +L2me Y22me = X2e .. ⎪ . ⎪ ⎪ ⎩ e e L1Y2me 1 +L2eY2eme 2 + . . . +L2eme Y2eme 2me = X2em2e .
(15)
i Assume the set of xα (2 < α < e − 1) that corresponds to ue = δ2i me is δ2n1 Pα1 , the set of j j α e n x (2 < α < e − 1) that corresponds to u = δ2me is δ2 1 Pα2 . Let i Pα1 = {iα1 , iα2 , . . . , iατiα }, j Pα2 = { jα1 , jα2 , . . . , jατ jα },
then it is obvious that τiα = τ jα . Hence, Yij can be expressed as follows: i i Yiej = L 1 (P11 , j11 )Yiej11 + L 1 (P11 , j12 )Yiej22 e τ i +··· + L 1 (P11 , j1τ j1 )Yi j1 j1 .
(16) k1
j1k
When x 1 = δ2n 1 , assume the corresponding input to 2 is u2 = δ21m2 , hence, i i Yiej1k1 = Lk21 (P21 , j21 )Yiej21 + Lk21 (P21 , j22 )Yiej22 1 1 e τ i +··· + Lk21 (P21 , j2τ j2 )Yi j2 j2 .
(17)
1
j1k
j2k
k2
If we assume the input of 3 that corresponds to x 1 = δ2n11 , x 2 = δ2n22 is u3 = δ22m3 , therefore, i i Yiej2k2 = Lk32 (P31 , j31 )Yiej31 + Lk32 (P31 , j32 )Yiej32 2 2 e τ i +··· + Lk32 (P31 , j3τ j3 )Yi j3 j3 . (18) 2
3004
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Similar to the above analysis, assume the input to e−1 that corresponds to x 1 = δ2j1n1k1 , k e−2
j(e−2)k
j2k
x 2 = δ2n22 ,..., x e−2 = δ2ne−2 e−2 is δ2em−2e−1 , then the following result can be obtained, i e 1 i −2 Yiej(eke−2) = Lkee−1 Lkee−1 −2 (P(e−1)1 , j (e−1)1 )Yi j(e−1) + −2 (P(e−1)1 , j (e−1)2 ) · e−2 e−2 e τ j(e−1) i Yiej(e2 −1) + · · · + Lkee−1 −2 (P(e−1)1 , j (e−1)τ j(e−1) )Yi j(e−1) . Let
t
θ=1
e−2
(19)
t
nθ = N , thus
Yiej(es −1) = S e( e−2 ( j γ =1
γ kγ
−1)·2N
e−1 −N γ
)
.
(20)
Based on Eqs. (16)–(20), the sum of elements in Yiej can be obtained. By solving Eq. (15), the least square solution of Yiej can be gotten. In a similar manner, the sum of the elements in Yiej1k1 can be gained, and the least square solution of Yiej1k1 can be received by solving Eq. (16). And so on, the least square solution of Yiej(es −1) which meets condition that the sum e−1 of elements in Yiej(es −1) is s
can be achieved by solving equation. e−2 N e−1 −N γ γ =1 ( jγ kγ
−1)·2
Through the above analysis, an algorithm to calculate the steady-state distributions of PBN can be achieved. Algorithm. Step 1: If rank(L 1 − I2n1 ) 2n1 − 1, calculate the nonnegative right eigenvectors that corresponds to eigenvalue 1 with l1 norm equal to 1, i.e., the steady-state distributions of system 1 . For each eigenvector obtained in this step, proceed the following process. Else, the PBN has no steady-state distribution. Step 2: Based on the steady-state distributions obtained in Step 1, solve Eqs. (11)–(13), we can obtain the steady-state distributions of PBN 1 × 2 . If it exists, go to next step. Else, the PBN has no steady-state distribution. Repeat the above process. Step μ − 1: On the basis of the steady-state distributions obtained in Step μ − 2, and the above method, we can obtain the steady-state distributions of PBN 1 × 2 × · · · × μ−1 . If it exists, go to next step. Else, the PBN has no steady-state distribution. Step μ: Based on the steady-state distributions obtained in Step μ − 1, as well as the method presented above, we can obtain the steady-state distributions of PBN 1 × 2 × · · · × μ . If it exists, the steady-state distribution of the whole network is achieved. Else, the PBN has no steady-state distribution. 4. Numerical example Consider the following PBN, which contains 11 nodes x 1, x 2, . . . , x 11, its network aggregation graph is shown in Fig. 2. This graph only presents the relation between blocks. Let χ1 = {x 1, x 2, x 3, x 4}, χ2 = {x 5, x 6, x 7}, χ3 = {x 8, x 9, x 10, x 11}, and calculate its steady-state distribution. The probabilistic transition matrices of 1 , 2 and 3 can be found in [31]. For PBN 1 , with nodes set χ 1 , its steady-state distribution S1 can be computed as S1 = [0.0677, 0.0678, 0.0639, 0.0766, 0.0600, 0.0616, 0.0701, 0.0502, 0.0502, 0.0548, 0.0650, 0.0696, 0.0474, 0.0841, 0.0504, 0.0607]T . For system 2 , its nodes set is χ 2 , input is x4 , and its probabilistic transition matrix is L 2 = [L12 , L22 ].
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3005
Fig. 2. Example of aggregation of a network with 11 nodes.
Through Eq. (11), we have L12 X12
+ L22 X22 = X12 + X22 ,
8 1 1 2 and the sum of elements in X12 is 8j=1 s2· j−1 , the sum of elements in X2 is j=1 s2· j , hence, the least square solutions of X12 and X22 can be gained. X12 = [0.0775, 0.0849, 0.0834, 0.0370, 0.0398, 0.0206, 0.0643, 0.0673]T , X22 = [0.0705, 0.0419, 0.0612, 0.0958, 0.0667, 0.0799, 0.0609, 0.0484]T . Through Eq. (12), we have 2 2 L1 Y11 + L22Y122 = X12 , L12Y212 + L22Y222 = X22 , By solving Eq. (21), least square solutions Yi2j , i, j ∈ {1, 2} are obtained, where Y112 = [0.0288, 0.0364, 0.0286, 0.0252, 0.0234, 0.0231, 0.0266, 0.0258]T , Y122 = [0.0321, 0.0343, 0.0318, 0.0326, 0.0289, 0.0401, 0.0291, 0.0282]T , Y212 = [0.0306, 0.0278, 0.0303, 0.0308, 0.0314, 0.0321, 0.0306, 0.0304]T , Y222 = [0.0354, 0.0341, 0.0358, 0.0351, 0.0358, 0.0327, 0.0357, 0.0363]T . Let S12 = [S 21 , S 23 , S 25 , S 27 , S 29 , S 211 , S 213 , S 215 ]T , S22 = [S 22 , S 24 , S 26 , S 28 , S 210 , S 212 , S 214 , S 216 ]T , based on Eq. (13), it is obvious to see that Y112 = [0.3321, 0.5721, 0.3956, 0.5011, 0.6023, 0.4197, 0.5425, 0.4940]S12 , Y212 = [0.6679, 0.4279, 0.6044, 0.4989, 0.3797, 0.5803, 0.4575, 0.5060]S12 ,
(21)
3006
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
Y122 = [0.4971, 0.5790, 0.3935, 0.5133, 0.5012, 0.4737, 0.3237, 0.5290]S22 , Y222 = [0.5029, 0.4210, 0.6065, 0.4867, 0.4988, 0.5263, 0.6763, 0.4710]S22 , then the steady-state distribution S2 of PBN 1 × 2 can be gotten via the above equation. S2 = [0.0085, 0.0083, 0.0076, 0.0090, 0.0066, 0.0083, 0.0062, 0.0065, 0.0086, 0.0097, 0.0078, 0.0064, 0.0070, 0.0088, 0.0106, 0.0077, 0.0082, 0.0096, 0.0077, 0.0098, 0.0082, 0.0085, 0.0073, 0.0073, 0.0087, 0.0100, 0.0077, 0.0066, 0.0071, 0.0089, 0.0104, 0.0079, 0.0084, 0.0082, 0.0075, 0.0089, 0.0065, 0.0082, 0.0061, 0.0064, 0.0086, 0.0096, 0.0079, 0.0064, 0.0070, 0.0088, 0.0107, 0.0077, 0.0084, 0.0075, 0.0073, 0.0084, 0.0057, 0.0079, 0.0055, 0.0060, 0.0086, 0.0098, 0.0078, 0.0064, 0.0070, 0.0088, 0.0106, 0.0077, 0.0084, 0.0072, 0.0073, 0.0082, 0.0054, 0.0078, 0.0052, 0.0058, 0.0081, 0.0090, 0.0076, 0.0059, 0.0065, 0.0084, 0.0106, 0.0072, 0.0086, 0.0072, 0.0074, 0.0082, 0.0053, 0.0079, 0.0053, 0.0058, 0.0095, 0.0111, 0.0080, 0.0074, 0.0079, 0.0095, 0.0103, 0.0088, 0.0084, 0.0078, 0.0074, 0.0086, 0.0061, 0.0080, 0.0058, 0.0062, 0.0081, 0.0090, 0.0076, 0.0059, 0.0065, 0.0084, 0.0106, 0.0072, 0.0083, 0.0076, 0.0073, 0.0084, 0.0059, 0.0079, 0.0056, 0.0060, 0.0081, 0.0089, 0.0076, 0.0058, 0.0064, 0.0084, 0.0107, 0.0071]T . For system 3 , its nodes set is χ 3 , input is x6 , and the probabilistic transition matrix is L 3 = [L13 , L23 ]. Through the above analysis, we have L13 X13 + L23 X23 = X13 + X23 ,
(22)
and the sum of elements in Xi3 is the sum of elements in S2 that corresponds to x6 = δ2i . The least square solutions of X13 , X23 can be received. X13 = [0.0297, 0.0290, 0.0282, 0.0312, 0.0362, 0.0223, 0.0323, 0.0320, 0.0371, 0.0341, 0.0347, 0.0341, 0.0338, 0.0376, 0.0345, 0.0289]T , X23 = [0.0278, 0.0272, 0.0258, 0.0285, 0.0330, 0.0210, 0.0296, 0.0290, 0.0343, 0.0320, 0.0331, 0.0319, 0.0334, 0.0365, 0.0336, 0.0279]T .
According to Eq. (15), next we will solve the following equation: L13Y113 + L23Y123 = X13 , L13Y213 + L23Y223 = X23 ,
(23)
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3007
Its least square solution is shown as follows: Y113 = [0.0129, 0.0148, 0.0147, 0.0170, 0.0159, 0.0092, 0.0134, 0.0149, 0.0156, 0.0174, 0.0164, 0.0164, 0.0195, 0.0166, 0.0169, 0.0174]T , Y123 = [0.0102, 0.0169, 0.0122, 0.0181, 0.0125, 0.0143, 0.0162, 0.0127, 0.0222, 0.0159, 0.0174, 0.0163, 0.0122, 0.0169, 0.0176, 0.0167]T , Y213 = [0.0199, 0.0162, 0.0156, 0.0157, 0.0160, 0.0174, 0.0165, 0.0154, 0.0185, 0.0182, 0.0154, 0.0186, 0.0139, 0.0190, 0.0176, 0.0128]T , Y223 = [0.0153, 0.0096, 0.0173, 0.0107, 0.0197, 0.0118, 0.0120, 0.0133, 0.0155, 0.0151, 0.0139, 0.0170, 0.0183, 0.0158, 0.0174, 0.0137]T . Because the following process is tedious, there is no need to enumerate all the discussions here. We only write the finial result which is achieved by the method presented in this paper, i.e., the steady-state distribution S3 of PBN 1 × 2 × 3 . 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 S3 = [γ11 , γ21 , φ11 , φ21 , γ31 , γ41 , φ31 , φ41 , γ51 , γ61 , φ51 , φ61 , γ71 , γ81 , φ71 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ81 , γ12 , γ22 , φ12 , φ22 , γ32 , γ42 , φ32 , φ42 , γ52 , γ62 , φ52 , φ62 , γ72 , γ82 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ72 , φ82 , γ13 , γ23 , φ13 , φ23 , γ33 , γ43 , φ33 , φ43 , γ53 , γ63 , φ53 , φ63 , γ73 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ83 , φ73 , φ83 , γ14 , γ24 , φ14 , φ24 , γ34 , γ44 , φ34 , φ44 , γ54 , γ64 , φ54 , φ64 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ74 , γ84 , φ74 , φ84 , γ15 , γ25 , φ15 , φ25 , γ35 , γ45 , φ35 , φ45 , γ55 , γ65 , φ55 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ65 , γ75 , γ85 , φ75 , φ85 , γ16 , γ26 , φ16 , φ26 , γ36 , γ46 , φ36 , φ46 , γ56 , γ66 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 φ56 , φ66 , γ76 , γ86 , φ76 , φ86 , γ17 , γ27 , φ17 , φ27 , γ37 , γ47 , φ37 , φ47 , γ57 , 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 γ67 , φ57 , φ67 , γ77 , γ87 , φ77 , φ87 , γ18 , γ28 , φ18 , φ28 , γ38 , γ48 , φ38 , φ48 , 3 3 3 3 3 3 3 3 T γ58 , γ68 , φ58 , φ68 , γ78 , γ88 , φ78 , φ88 ] .
In order to save paper, all the elements of S3 that can be found in [31] are omitted. 5. Conclusion In this paper, we have considered the steady-state distributions of PBNs with special structure. To solve this problem, the network has been partitioned into several blocks by cascading aggregation. Then finding steady-state distributions of the whole networks has been converted to computing several corresponding equations. Two necessary and sufficient conditions for the existence of the steady-state distributions for PBNs have been obtained. An algorithm has been presented to calculate the steady-sate distributions for the whole networks. In fact, our method can deal with large-scale PBNs theoretically, but due to the complexity of computation, we give an example with 11 nodes. Acknowledgment This work was supported by National Natural Science Foundation under Grants 61773371 and 61877036, and China Postdoctoral Science Foundation under Grant 2016M-602143.
3008
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
References [1] S. Kim, S. Imoto, S. Miyano, Dynamic Bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data, Biosystems 75 (2004) 57–65. [2] S. Kauffman, Metabolic stability and epigenesis in randomly constructed genetic nets, J. Theor. Biol. 244 (4) (1969) 670–679. [3] L. Steggles, R. Banks, A. Wipat, Modelling and analysing genetic networks: from Boolean networks to petri nets 4210 (2006) 127–141. [4] W. Ching, M. Michael, S. Eric, T. Akutsu, On construction of stochastic genetic networks based on gene expression sequences, Int. J. Neural Syst. 15 (4) (2005) 297–310. [5] S. Zhang, W. Ching, N. Tsing, H. Leung, D. Guo, A new multiple regression approach for the construction, Artif. Intell. Med. 48 (2010) 153–160. [6] I. Shmulevich, E.R. Dougherty, S. Kim, W. Zhang, Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks, Bioinformatics 18 (2002) 261–274. [7] D. Cheng, H. Qi, Z. Li, Analysis and Control of Boolean Networks: A Semi-tensor Product Approach, Springer, 2011. [8] M. Meng, L. Liu, G. Feng, Stability and L1 gain analysis of Boolean networks with Markovian jump parameters, IEEE Trans. Autom. Control 62 (8) (2017) 4222–4228. [9] H. Li, Y. Wang, Z. Wu, Further results on feedback stabilization control design of Boolean control networks, Automatica 83 (2017) 303–308. [10] Y. Yu, J. Feng, J. Pan, Block decoupling of Boolean control networks, IEEE Trans. Autom. Control (2018), doi:10.1109/TAC.2018.2880411. [11] S. Wang, J. Feng, J. Zhao, J. Xia, Controllability decomposition of dynamic-algebraic Boolean control networks, Int. J. Control (2018), doi:10.1080/00207179.2018.1527040. [12] Y. Zhao, D. Cheng, On controllability and stabilizability of probabilistic Boolean control networks, Sci. China 57 (1) (2014) 1–14. [13] Y. Liu, H. Chen, J. Lu, B. Wu, Controllability of probabilistic Boolean control networks based on transition probability matrices, Automatica 52 (2015) 340–345. [14] K. Zhang, L. Zhang, Controllability of probabilistic Boolean control networks with time-variant delays in states, in: Proceedings of the 52th IEEE Conference on Decision and Control, 2013, pp. 7211–7216. [15] M. Meng, J. Lam, J. Feng, K.C. Cheung, Stability and stabilization of Boolean networks with stochastic delays, IEEE Trans. Autom. Control (2018), doi:10.1109/TAC.2018.2835366. [16] Z. Li, J. Song, J. Yang, Partial stability of probabilistic Boolean network, in: Proceedings of the 26th Chinese Control and Decision Conference, 2014, pp. 1952–1956. [17] J. Zhao, Z. Liu, Observability of probabilistic Boolean networks, in: Proceedings of the 34th Chinese Control Conference, 2015, pp. 183–186. [18] H. Chen, J. Liang, J. Lu, J. Qiu, Synchronization for the realization-dependent probabilistic Boolean networks, IEEE Trans. Neural Netw. Learn. Syst. 29 (4) (2018) 819–831. [19] J. Lu, J. Zhong, L. Li, D.W.C. Ho, J. Cao, Synchronization analysis of master-slave probabilistic Boolean networks, Sci. Rep. 5 (2015) 13437. [20] Z. Liu, Y. Wang, Optimal finite-horizon control problem of context-sensitive probabilistic Boolean networks with perturbation, in: Proceedings of the 31th Chinese Control Conference, 2012, pp. 140–145. [21] M. Yang, R. Li, T. Chu, A new method and application for controlling steady-state probability distributions of probabilistic Boolean networks, in: Proceedings of the 2014 IEEE Congress on Evolutionary Computation, 2014, pp. 1490–1495. [22] X. Cheng, Y. Qiu, W. Hou, W.K. Ching, A semi-tensor product approach for probabilistic Boolean networks, in: Proceedings of the 8th International Conference on System Biology, 2014, pp. 85–90. [23] R. Pal, Context-sensitive probabilistic Boolean networks: steady-state properties, reduction, and steady-state approximation, IEEE Trans. Signal Process. 58 (2) (2010) 879–890. [24] I. Shmulevich, I. Gluhovsky, R.F. Hashimoto, E.R. Dougherty, W. Zhang, Steady-state analysis of genetic regulatory networks modelled by probabilistic Boolean networks, Compar. Funct. Genom. 4 (2003) 601–608. [25] W.K. Ching, S. Zhang, M.K. Ng, T. Akutsu, An approximation method for solving the steady-state probability distribution of probabilistic Boolean networks, Bioinformatics 23 (12) (2007) 1511–1518. [26] I. Apostolopoulou, D. Marculescu, Tractable learning and inference for large-scale probabilistic Boolean networks, IEEE Trans. Neural Netw. Learn. Syst. (2019), doi:10.1109/TNNLS.2018.2886207.
J. Pan, J.-e. Feng and M. Meng / Journal of the Franklin Institute 356 (2019) 2994–3009
3009
[27] Y. Zhao, J. Kim, M. Filippone, Aggregation algorithm towards large-scale Boolean network analysis, in: Proceedings of the IEEE Transactions on Automatic Control, 58, 2013, pp. 1976–1985. [28] Y. Zhao, B.K. Ghosh, D. Cheng, Control of large-scale Boolean networks via network aggregation, IEEE Trans. Neural Netw. Learn. Syst. 27 (7) (2016) 1527–1536. [29] K. Zhang, Observability and reconstructibility of large-scale Boolean control networks via network aggregations, (2017) arXiv preprint, arXiv:1704.03231. [30] J.G. Kemeny, J.L. Snell, Finite Markov Chains, Springer, 1983. [31] J. Pan, Control and Applications of Boolean Networks via Algebraic Method, 2018 Ph.D thesis.