A Feedback Vertex Set-Based Approach to Simplifying Probabilistic Boolean Networks

Koichi KOBAYASHI

Summary :

A PBN is well known as a mathematical model of complex network systems such as gene regulatory networks. In Boolean networks, interactions between nodes (e.g., genes) are modeled by Boolean functions. In PBNs, Boolean functions are switched probabilistically. In this paper, for a PBN, a simplified representation that is effective in analysis and control is proposed. First, after a polynomial representation of a PBN is briefly explained, a simplified representation is derived. Here, the steady-state value of the expected value of the state is focused, and is characterized by a minimum feedback vertex set of an interaction graph expressing interactions between nodes. Next, using this representation, input selection and stabilization are discussed. Finally, the proposed method is demonstrated by a biological example.

• Publication: IEICE TRANSACTIONS on Fundamentals Vol.E107-A No.5 pp.779-785
• Publication Date: 2024/05/01
• Publicized: 2023/09/26
• Online ISSN: 1745-1337
• DOI: 10.1587/transfun.2023MAP0004
• Type of Manuscript: Special Section PAPER (Special Section on Mathematical Systems Science and its Applications)

1.  Introduction

Control of complex network systems such as gene regulatory networks is one of the central problems in control theory. To simplify such systems, approximation by discrete models is frequently effective. In analysis and control of gene regulatory networks, a Boolean network (BN) is well known as one of the discrete models [1]. In a BN, expression of a gene (node) is given by a binary value (ON/OFF), and interactions between genes are modeled by Boolean functions. Furthermore, a BN is extended to a probabilistic Boolean network (PBN). In a PBN, the candidates of Boolean functions are given in advance, and one is probabilistically chosen from the candidates [2].

In the last decade, several control methods for PBNs have been developed (see, e.g., [3]-[11]). These methods have drawbacks. For example, in both the Markov chain-based method [7]-[9] and the semi-tensor product method [5], [6], [11], \(2^n \times 2^n\) matrices must be manipulated, where \(n\) is the number of nodes. In the methods of [3], [4], an integer programming problem or a polynomial optimization problem must be solved. It is important to develop a simple and useful method for analysis and control of a PBN.

In this paper, we study a new approach to simplifying a PBN focusing on interactions between nodes.

First, we propose a simplified representation, which can be derived from the polynomial representation proposed in [4]. The polynomial representation represents the expected value of the state of a given PBN. However, in [4], the optimal control problem has been mainly focused, and simplification of the polynomial representation has not been studied. To analyze the behavior of a PBN, it is important to consider a simplified representation that some properties are preserved. The simplified representation proposed here preserves the steady-state properties. To derive it, we use a minimum feedback vertex set of an interaction graph expressing interactions between nodes, where each vertex of an interaction graph corresponds to each node (see Sect. 2 for further details). Using the proposed representation, the steady-state value of the expected value of the state can be characterized by only vertices in a minimum feedback vertex set. In general, the number of vertices in a minimum feedback vertex set is much smaller than the number of vertices in an interaction graph. Hence, the dynamics of a given PBN can be simplified. For a BN, such a method has been proposed in e.g., [12], [13]. In addition, the significance of feedback vertex sets in analysis of BNs has been pointed out in [14]. In [12]-[14], only a BN has been focused, and a PBN has not been considered. The proposed method may be regarded as an extension of these methods to a PBN.

Next, the simplified representation can be applied to the input selection problem. In many cases, gene regulatory networks have no (external) control inputs. As the first step of studying control of such systems, it is important to choose control nodes from nodes. In the control node, we assume that the initial binary value can be arbitrarily set (i.e., we consider the constant control input). In the proposed method, the vertices in a minimum feedback vertex set are regarded as control nodes (see also [15], [16]). By a simplified representation, stability can be easily analyzed, and the constant control inputs can be derived.

Finally, we demonstrate a derivation method of constant control inputs by an example of a neurotransmitter signaling pathway [17].

This paper is organized as follows. First, the outline of PBNs is explained, and interaction graphs and minimum feedback vertex sets are defined. Next, after a polynomial representation of PBNs is summarized, a simplified representation is proposed. The input selection problem using it is also explained. Finally, a biological example is presented.

Notation: Let \(\{ 0,1 \}^{n}\) denote the set of \(n\)-dimensional vectors, which consists of elements \(0\) and \(1\). For the \(n\)-dimensional vector \(x = [ x_1~~x_2~~\cdots\) \(x_n ]^{\top}\) and the index set \({\mathcal I} = \{ i_1, i_2, \dots, i_m \} \subseteq \{ 1,2,\dots,n\}\), define \([x_i]_{i \in {\mathcal I}} := [ x_{i_1}~~x_{i_2}~~\cdots~~x_{i_m} ]^{\top}\).

2.  Probabilistic Boolean Networks

In this section, first, a PBN is briefly explained.

In a PBN representing a complex network network with \(n\) nodes, each node has a binary value \(x_i \in \{ 0,1 \}\), \(i \in \{ 1,2,\dots,n \}\). The dynamics for \(x_i\) are given by

\[\begin{align} x_i(k+1) &= f^{(i)}([x_j(k)]_{j \in {\mathcal N}^{(i)}}) \nonumber \\ &= \begin{cases} f^{(i)}_1([x_j(k)]_{j \in {\mathcal N}^{(i)}_1}), \mbox{Prob. $c^{(i)}_1$,} \\ f^{(i)}_2([x_j(k)]_{j \in {\mathcal N}^{(i)}_2}), \mbox{Prob. $c^{(i)}_2$,} \\ \vdots \\ f^{(i)}_{q(i)}([x_j(k)]_{j \in {\mathcal N}^{(i)}_{q(i)}}), \mbox{Prob. $c^{(i)}_{q(i)}$,} \end{cases} \tag{1} \end{align}\]

where \(k = 0,1,2,\dots\) is the discrete time, the set \({\mathcal N}^{(i)}_l \subseteq \{ 1,2,\dots,n \}\), \(l=1,2,\dots,q(i)\) is a given index set, and the set \({\mathcal N}^{(i)}\), \(i=1,2,\dots,n\) is defined by

\[{\mathcal N}^{(i)} := \bigcup_{l=1}^{q(i)} {\mathcal N}^{(i)}_{l}.\]

The function \(f^{(i)}_l : \{ 0,1 \}^{|{\mathcal N}^{(i)}_l|} \rightarrow \{ 0,1 \}^1\) is a given Boolean function consisting of logical operators such as AND (\(\wedge\)), OR (\(\vee\)), and NOT (\(\neg\)). If \({\mathcal N}^{(i)}_l = \emptyset\) holds, then the value (\(0\) or \(1\)) of the Boolean function \(f^{(i)}_l\) is uniquely determined as \(0\) or \(1\). The probability \(c^{(i)}_l\), \(l=1,2,\dots,q(i)\) is formally defined by

\[c^{(i)}_l := {\rm Prob} \left( f^{(i)}=f^{(i)}_l \right).\]

For \(c^{(i)}_l\), the following relation:

\[\begin{equation*} \sum_{l=1}^{q(i)} c^{(i)}_l = 1 \tag{2} \end{equation*}\]

must be satisfied. We define the state \(x(k) := [ x_1(k) x_2(k) \cdots\) \(x_n(k) ]^{\top} \in \{ 0,1 \}^n\).

We present a simple example.

Example 1: Consider the PBN in which Boolean functions and probabilities are given by

\[\begin{aligned} f^{(1)} &= \begin{cases} f^{(1)}_1 = x_3(k), c^{(1)}_1 = 0.8, \\ f^{(1)}_2 = \neg x_3(k), c^{(1)}_2 = 0.2, \end{cases} \\ f^{(2)} &= f^{(2)}_1 = x_1(k) \wedge \neg x_3(k), c^{(2)}_1 = 1.0, \\ f^{(3)} &= \begin{cases} f^{(3)}_1 = x_1(k) \wedge \neg x_2(k), c^{(3)}_1 = 0.7, \\ f^{(3)}_2 = x_2(k), c^{(3)}_2 = 0.3, \end{cases} \end{aligned}\]

where \(q(1) = 2\), \(q(2) = 1\) and \(q(3) = 2\) hold, \({\mathcal N}^{(1)} = \{ 3 \}\), \({\mathcal N}^{(2)} = \{ 1,3 \}\), and \({\mathcal N}^{(3)} = \{ 1,2 \}\) hold, and we see that the relation (2) is satisfied. Next, consider the state trajectory. Then, for \(x(0) = [ 0~~ 0~~ 0 ]^{\top}\), we obtain

\[\begin{aligned} {\rm Prob} \left( x(1) = [ 0 ~~ 0 ~~ 0 ]^{\top} | x(0) = [ 0 ~~ 0 ~~ 0 ]^{\top} \right) &= 0.8, \\ {\rm Prob} \left( x(1) = [ 1 ~~ 0 ~~ 0 ]^{\top} | x(0) = [ 0 ~~ 0 ~~ 0 ]^{\top} \right) &= 0.2. \end{aligned}\]

In this example, the cardinality of the finite state set \(\{ 0,1 \}^3\) is given by \(2^3 = 8\), and we obtain the state transition diagram of Fig. 1 by computing the transition from each state. In Fig. 1, the number assigned to each node denotes \(x_1\), \(x_2\), \(x_3\), and the number assigned to each arc denotes the transition probability from some state to other state. Note here that for simplicity, the state transition from only \(x(k) = [ 0 ~~ 0 ~~ 0 ]^{\top}, \ [ 0 ~~ 0 ~~ 1 ]^{\top}, \ [ 0 ~~ 1 ~~ 0 ]^{\top}, \ [ 1 ~~ 1 ~~ 0 ]^{\top}\) is illustrated in Fig. 1. \(\tag*{◻}\)

Fig. 1  State transition diagram of the PBN in Example 1. For simplicity, the state transition from only \(x(k) = [ 0 ~~ 0 ~~ 0 ]^{\top}, \ [ 0 ~~ 0 ~~ 1 ]^{\top}, \ [ 0 ~~ 1 ~~ 0 ]^{\top}, \ [ 1 ~~ 1 ~~ 0 ]^{\top}\) is illustrated.

Next, we define an interaction graph for a given PBN as follows.

Definition 1: An interaction graph of a given PBN is defined by a directed graph \(G = ({\mathcal V}, {\mathcal E})\), where \({\mathcal V} = \{ 1,2,\dots,n \}\) is the set of vertices corresponding to \(x_i\), \(i \in \{ 1,2,\dots,n \}\), and \({\mathcal E} = \{ (j,i) \in \{ 1,2,\dots,n \} \times \{ 1,2,\dots,n \} | j \in {\mathcal N}^{(i)} \}\) is the set of arcs.

For a given interaction graph, a feedback vertex set is defined as follows (see, e.g., [18], [19]).

Definition 2: A set of vertices of an interaction graph is called a feedback vertex set if removal of vertices results an acyclic graph. In particular, a feedback vertex set is called a minimum feedback vertex set if the number of its elements is minimum.

The computational complexity of finding a minimum feedback vertex set is NP-complete [19], but an approximate algorithm of finding it has been developed (see, e.g., [18]).

We present two examples.

Example 2: Consider the PBN in Example 1 again. From \({\mathcal N}^{(1)} = \{ 3 \}\), \({\mathcal N}^{(2)} = \{ 1,3 \}\), and \({\mathcal N}^{(3)} = \{ 1,2 \}\), its interaction graph is given by Fig. 2. In addition, the minimum feedback vertex set is given by \(\{ 3 \}\). \(\tag*{◻}\)

Fig. 2  Interaction graph of the PBN in Example 1.

Example 3: Consider the interaction graph of an apoptosis network given by Fig. 3. See [20] for details. From Fig. 3, we see that one of the minimum feedback vertex sets is given by {NF\(\kappa\)B\(_{\rm nuc}\), C3a, TNF}. \(\tag*{◻}\)

Fig. 3  Interaction graph of an apoptosis network [20] in Example 3.

3.  Polynomial Representation of Probabilistic Boolean Networks

In this section, we explain a polynomial representation of a given PBN [4]. Using it, the expected value of the state (\(E[x_i(k)|\ast]\)) can be represented by a polynomial system.

As a preparation, the following lemma [21] is introduced.

Lemma 1: Consider two binary variables \(\delta_1\) and \(\delta_2\). Then the following relations hold.

(i) \(\neg \delta_1\) is equivalent to \(1 - \delta_1\).

(ii) \(\delta_1 \wedge \delta_2\) is equivalent to \(\delta_1 \delta_2\).

(iii) \(\delta_1 \vee \delta_2\) is equivalent to \(\delta_1 + \delta_2 - \delta_1 \delta_2\).

By Lemma 1, a given Boolean function can be transformed into a polynomial on the real number field. For example, the logical formula \(\delta_1 \vee \neg \delta_2\) is equivalently transformed into the polynomial \(\delta_1 + (1-\delta_2) - \delta_1 (1-\delta_2) = 1 - \delta_2 + \delta_1 \delta_2\).

Hereafter, we assume that the Boolean function \(f^{(i)}_j\) is expressed as a polynomial with binary variables. For simplicity of notation, the condition in \(E[x_i(k)|\ast]\) is omitted. Then, the following result has been obtained in [4].

Theorem 1: Suppose that for the PBN (1), the initial state \(x(0)=x_0\) is arbitrarily given. Then, the expected value of the state, \(E[x_i(k)] \in [0,1]\) is expressed as the following polynomial system:

\[\begin{align} E[x_i(k+1)] &= \tilde{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal N}^{(i)}}]), \tag{3} \\ \tilde{f}^{(i)} &= \sum_{l=1}^{q(i)} c^{(i)}_l f^{(i)}_l(E[[x_j(k)]_{j \in {\mathcal N}_j^{(i)}}]). \nonumber \end{align}\]

We present a simple example.

Example 4: Consider the PBN in Example 1 again. Then, the polynomial system for \(x_1\) is derived as

\[\begin{aligned} E[x_1(k+1)] &= 0.8 E[x_3(k)] + 0.2 (1-E[x_3(k)]) \\ &= 0.2 + 0.6 E[x_3(k)]. \end{aligned}\]

In a similar way, the polynomial system for \(x_2\) is derived as

\[\begin{aligned} E[x_2(k+1)] &= E[x_1(k)] - E[x_1(k)] E[x_3(k)]. \end{aligned}\]

Finally, the polynomial system for \(x_3\) is derived as

\[\begin{aligned} E[x_3(k+1)] =& \ 0.7 ( E[x_1(k)] - E[x_1(k)] E[x_2(k)] ) \\ & + 0.3 E[x_2(k)] \\ =& \ 0.7 E[x_1(k)] + 0.3 E[x_2(k)] \\ & - 0.7 E[x_1(k)] E[x_2(k)]. \end{aligned}\]

\(\tag*{◻}\)

We remark here that the polynomial system (3) is a class of positive systems under a binary initial state.

4.  Simplified Representation of Probabilistic Boolean Networks

In this section, we consider simplifying the polynomial representation (3) for a given PBN. Here, we focus on the steady state of the system (3).

As a preparation, let \({\mathcal W}\) denote a minimum feedback vertex set of the interaction graph \(G\). Let \(G_I = ({\mathcal V}, {\mathcal E}_I)\), \({\mathcal E}_I := \{ (i,j) | (j,i) \in {\mathcal E} \}\) denote the directed graph obtained by flipping the start and end nodes of the arc of the interaction graph \(G\). For \(G_I\), let \(d_i\) denote the length of the longest path/cycle from the node \(i\) to one vertex of \({\mathcal W}\). For example, in the case of the interaction graph in Fig. 2, since \({\mathcal W}\) is given by \({\mathcal W} = \{ 3 \}\), \(d_1 = 1\) holds. For the node \(2\), there are two paths (\(2 \rightarrow 1 \rightarrow 3\) and \(2 \rightarrow 3\)), and \(d_2 = 2\) holds. For the node \(3\), there are three paths (\(3 \rightarrow 1 \rightarrow 3\), \(3 \rightarrow 2 \rightarrow 3\), and \(3 \rightarrow 2 \rightarrow 1 \rightarrow 3\)), and \(d_3 = 3\) holds. We remark that even if the node \(i\) is included in \({\mathcal W}\), we must find paths.

In addition, let \({\mathcal W}_{i,l}\), \(l=0,1,\dots,d_i-1\) denote a subset of \({\mathcal W}\) such that there exist a path from its element to the node \(i\) with the length \(l+1\) over the graph \(G\). For example, in the case of the interaction graph in Fig. 2 (\({\mathcal W} = \{ 3 \}\)), we can obtain \({\mathcal W}_{1,0} = \{ 3 \}\), \({\mathcal W}_{2,0} = {\mathcal W}_{2,1} = \{ 3 \}\), \({\mathcal W}_{3,0} = \emptyset\), and \({\mathcal W}_{3,1} = {\mathcal W}_{3,2} = \{ 3 \}\).

The proposed procedure is given as follows.

Procedure for deriving a simplified representation:

Step 1: For a given interaction graph, find a minimum feedback vertex set \({\mathcal W}\).

Step 2: Using \(E[x_i(k-l+1)] = \tilde{f}^{(i)}(E[[x_j(k-l)]_{j \in {\mathcal N}^{(i)}}])\), \(l = 1,2,\dots,d_i-1\), rewrite the polynomial \(\tilde{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal N}^{(i)}}])\) as the following equivalent polynomial:

\[\begin{eqnarray*} \tilde{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal W}_{i,0}}], E[[x_j(k-1)]_{j \in {\mathcal W}_{i,1}}],\nonumber\\ \dots,E[[x_j(k-d_i+1)]_{j \in {\mathcal W}_{i,d_i-1}}]). \tag{4} \end{eqnarray*}\]

Step 3: For the obtained polynomial, replace \(E[[x_j(k-l)]_{j \in {\mathcal W}_{i,l}}]\), \(l = 1,2,\dots,d_i-1\) with \(E[[x_j(k)]_{j \in {\mathcal W}_{i,l}}]\). Let \(\hat{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal W}_i}])\) denote the polynomial obtained, where \({\mathcal W}_i = \cup_{l} {\mathcal W}_{i,l}\).

Step 4: Obtain the simplified representation as

\[\begin{equation*} E[\hat{x}_i(k+1)] = \hat{f}^{(i)}(E[[\hat{x}_j(k)]_{j \in {\mathcal W}_i}]), \tag{5} \end{equation*}\]

where \(\hat{x}_i\) is newly defined as a binary variable.

In Step 2, \(\tilde{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal N}^{(i)}}])\) and (4) are equivalent. Because the polynomial systems that time is shifted are substituted into \(\tilde{f}^{(i)}(E[[x_j(k)]_{j \in {\mathcal N}^{(i)}}])\). From the definition of minimum feedback vertex sets, Step 2 can be necessarily performed. In the proposed simplified representation, the dynamics are expressed by only \(E[[\hat{x}_j(k)]_{j \in {\mathcal W}_i}]\). As was explained in Sect. 2, the computational complexity of finding a minimum feedback vertex set is NP-complete. Hence, the computational complexity of Step 1 is also NP-complete. In Step 2, for each \(i\), substitutions of polynomials to certain terms of other polynomials are performed \(d_i-1\) times. In Step 3, \(k-1, k-2, \dots, k-d_i+1\) are replaced with \(k\). In Step 4, the binary variable in the simplified representation is newly defined to distinguish between the original PBN and the simplified PBN. If a minimum feedback vertex set is given, Steps 2-4 can be implemented by a suitable software such as MATLAB and Mathematica.

We present a simple example.

Example 5: Consider the polynomial system obtained in Example 4. In Example 2, the minimum feedback vertex set is given by \(\{ 3 \}\). From Fig. 2, \(d_1=1\), \(d_2=2\), and \(d_3=3\) holds. First, consider the polynomial system for \(x_1\). In this system, only \(E[x_3(k)]\) is included as a state variable (this can be confirmed from also \({\mathcal W}_{1,0} = \{ 3 \}\)). Then, from \(d_1=1\), we cannot simplify it, and we can obtain the following simplified representation:

\[\begin{align} E[\hat{x}_1(k+1)] &= \hat{f}^{(1)}(E[\hat{x}_3(k)]) \nonumber \\ &= 0.2 + 0.6 E[\hat{x}_3(k)]. \tag{6} \end{align}\]

Next, consider the polynomial system for \(x_2\). From \(d_2=2\), the dynamics at \(k+1\), \(k\), and \(k-1\) can be derived as

\[\begin{aligned} E[x_2(k+1)] =& \ E[x_1(k)] - E[x_1(k)] E[x_3(k)]. \\ =& \ (0.2 + 0.6 E[x_3(k-1)]) \\ & \ \times (1- E[x_3(k)]). \end{aligned}\]

From \({\mathcal W}_{2,0} = {\mathcal W}_{2,1} = \{ 3 \}\), we see that the function \(\tilde{f}^{(2)}\) can represented by \(E[x_3(k-1)]\) and \(E[x_3(k)]\). By replacing \(E[x_3(k-1)]\) with \(E[x_3(k)]\), the following simplified polynomial can be derived:

\[\begin{align} E[\hat{x}_2(k+1)] =& \ \hat{f}^{(2)}(E[\hat{x}_3(k)]) \nonumber \\ =& \ 0.2 + 0.4 E[\hat{x}_3(k)] \nonumber \\ & \ - 0.6 E[\hat{x}_3(k)]^2. \tag{7} \end{align}\]

Finally, consider the polynomial system for \(x_3\). From \(d_3=3\), the dynamics at \(k+1\), \(k\), \(k-1\), and \(k-2\) can be derived as

\[\begin{aligned} E[x_3(k+1)] =& \ 0.7 ( 0.2 + 0.6 E[x_3(k-1)] ) \\ & \ + 0.3 ( 0.2 + 0.6 E[x_3(k-2)] ) \\ & \ ~~~~~~~~~\times ( 1 - E[x_3(k-1)] ) \\ & \ - 0.7 ( 0.2 + 0.6 E[x_3(k-1)] ) \\ & \ ~~~~~~~~~\times ( 0.2 + 0.6 E[x_3(k-2)] ) \\ & \ ~~~~~~~~~\times ( 1 - E[x_3(k-1)] ). \end{aligned}\]

From \({\mathcal W}_{3,0} = \emptyset\) and \({\mathcal W}_{3,1} = {\mathcal W}_{3,2} = \{ 3 \}\), we see that the function \(\tilde{f}^{(3)}\) can represented by \(E[x_3(k-2)]\) and \(E[x_3(k-1)]\). Then, the following simplified polynomial can be derived:

\[\begin{align} E[\hat{x}_3(k+1)] =& \ \hat{f}^{(3)}(E[\hat{x}_3(k)]) \nonumber \\ =& \ 0.172 + 0.4 E[\hat{x}_3(k)] \nonumber \\ & \ - 0.264 E[\hat{x}_3(k)]^2+ 0.252 E[\hat{x}_3(k)]^3. \tag{8} \end{align}\]

From (6), (7), and (8), we see that \(E[\hat{x}_i(k+1)]\), \(i=1,2,3\) can be expressed by a polynomial with respect to only \(E[\hat{x}_3(k)]\). \(\tag*{◻}\)

Next, we consider analyzing the simplified representation obtained. The notion of fixed points is defined as follows.

Definition 3: For the system (3), \(x_F := [ x_1^F x_2^F \cdots x_n^F ]^{\top}\) is called a fixed point if \(x_i^F = E[x_i(k+1)] = E[x_i(k)]\), \(i \in \{ 1, 2, \dots, n \}\) hold. In a similar way, for the system (5), \(\hat{x}_F := [ \hat{x}_1^F ~ \hat{x}_2^F ~ \cdots ~ \hat{x}_n^F ]^{\top}\) is called a fixed point if \(\hat{x}_i^F = E[\hat{x}_i(k+1)] = E[\hat{x}_i(k)]\), \(i \in \{ 1, 2, \dots, n \}\) hold.

We remark that \(x_F\) and \(\hat{x}^F\) are not uniquely determined in general. Then, we have the following theorem.

Theorem 2: The fixed point \(x^F\) for the system (3) is in a one-to-one correspondence with the fixed point \(\hat{x}^F\) for the system (5).

Proof: When \(x^F\) is calculated, the polynomial (4) may be used. The fixed point \(x^F\) can be derived from (4) and \(x_i^F = E[x_i(k+1)] = E[x_i(k)] = E[x_i(k-1)] = \cdots = E[x_i(k-d_i+1)]\), \(i \in \{ 1, 2, \dots, n \}\). In the derivation of the simplified representation (5), the polynomial (4) is simplified based on \(E[x_i(k)] = E[x_i(k-1)] = \cdots = E[x_i(k-d_i+1)]\), \(i \in \{ 1, 2, \dots, n \}\). Hence, there is a one-to-one correspondence between \(x^F\) and \(\hat{x}^F\). \(\tag*{◻}\)

From this theorem, we see that the simplified representation 5 can be used for analysis and control of PBNs focusing on the fixed point (i.e., the long-time behavior).

5.  Discussion on Input Selection and Stabilization

In this section, we discuss input selection and stabilization as one of the methods for utilization of the simplified representation (5). In [15], [16], a control method using a minimum feedback vertex set has been proposed for the system represented by ordinary differential equations. In this method, it has been shown that stabilization or destabilization can be achieved by controlling nodes in a minimum feedback vertex set. Motivated by this method, in this paper, the following assumption, in which only the initial value is focused for simplicity, is imposed for (5).

Assumption 1: The polynomial system \(E[\hat{x}_i(k+1)] = \hat{f}^{(i)}(\cdot)\), \(i \in {\mathcal W}\) in (5) can be replaced with

\[\hat{x}_i(k+1) = \hat{x}_i(k) = \hat{x}_i^0, i \in {\mathcal W},\]

where the initial value \(\hat{x}_i^0 \in \{ 0,1 \}^{| {\mathcal W} |}\) can be arbitrarily given.

In this paper, a vertex included in \({\mathcal W}\) is called a control node. This assumption implies that binary values in control nodes are regarded as constant control inputs. In other words, finding a minimum feedback vertex set of an interaction graph implies selecting of constant control inputs. Using the simplified representation (5), the value of constant control inputs can be easily determined.

We present a simple example.

Example 6: Consider the PBN in Example 1 again. Under Assumption 1, its simplified representation is given by (6), (7), and \(\hat{x}_3(k+1) = \hat{x}_3(k) = \hat{x}_3^0 \in \{ 0,1 \}\). In the case of \(\hat{x}_3^0=0\), we can obtain \(E[\hat{x}_1(k+1)]=E[\hat{x}_2(k+1)]=0.2\). In the case of \(\hat{x}_3^0=1\), we can obtain \(E[\hat{x}_1(k+1)]=0.8\) and \(E[\hat{x}_2(k+1)]=0\). Figures 4 and 5 show the simulation results using the PBN in Example 1, where the initial state is given by \(x_1(0)=x_2(0)=1\), and the simulation was repeated 10000 times. From these figures, we see that the simplified representation expresses the long-time behavior. Furthermore, in the case of \(\hat{x}_3^0=1\), we guarantee that \(x_2(k)=0\), \(k \geq 1\) holds with the probability \(1\). We can confirm this fact from both the simplified representation and the simulation result. \(\tag*{◻}\)

Fig. 4  Time response of the average of 10000 samples for the PBN in Example 1, where \(x_3(k)=0\). Solid line: \(x_1\). Broken line: \(x_2\).

Fig. 5  Time response of the average of 10000 samples for the PBN in Example 1, where \(x_3(k)=1\). Solid line: \(x_1\). Broken line: \(x_2\).

In stability analysis and stabilization of PBNs, a semi-tensor product (STP) method is well known as one of the typical methods (see, e.g., [5]). In the STP method, a PBN with \(n\) nodes and \(m\) control nodes is represented by a matrix with the size of \(2^n \times 2^{n+m}\). Hence, calculations on stability analysis and stabilization may be difficult for large-scale PBNs.

The proposed simplified representation can express the behavior by using fewer variables, and can be used in stability analysis and stabilization of PBNs. A PBN is complex, and there may be the case that the control input is not given explicitly. In such cases, the simplified representation provides us useful information on input selection. If stabilization cannot be performed under Assumption 1, then the nodes that do not converge to the origin must be controlled directly.

The value \(\hat{x}_i^0\) of constant control inputs such that a given PBN satisfies a certain specification can be easily calculated by using e.g., MATLAB and Mathematica.

6.  Biological Example

In this section, we present a biological example. Consider the BN model of a neurotransmitter signaling pathway [17]. This model expresses an interaction pathway between the glutamatergic and dopaminergic receptors. The BN model of this system is given by

\[\begin{aligned} x_1(k+1) &= x_1(k), \\ x_2(k+1) &= x_1(k) \wedge \neg x_3(k), \\ x_3(k+1) &= x_2(k), \\ x_4(k+1) &= x_2(k), \\ x_5(k+1) &= x_2(k), \\ x_6(k+1) &= x_4(k) \wedge \neg x_5(k), \\ x_7(k+1) &= x_6(k), \\ x_8(k+1) &= x_5(k), \\ x_9(k+1) &= x_8(k) \vee x_{14}(k), \\ x_{10}(k+1) &= x_9(k), \\ x_{11}(k+1) &= x_7(k) \wedge \neg x_{10}(k), \\ x_{12}(k+1) &= \neg x_{11}(k), \\ x_{13}(k+1) &= x_{13}(k), \\ x_{14}(k+1) &= x_7(k) \wedge \neg x_{12}(k) \wedge x_{13}(k), \\ x_{15}(k+1) &= x_{14}(k), \\ x_{16}(k+1) &= x_{15}(k). \end{aligned}\]

The interaction graph is given by Fig. 6. Here, to add probabilistic behavior, the dynamics of \(x_6\) is artificially changed to

\[x_6(k+1) = \begin{cases} x_4(k) \wedge \neg x_5(k), \mbox{Prob. 0.7,} \\ x_4(k), \mbox{Prob. 0.3.} \end{cases}\]

Fig. 6  Interaction graph of the Boolean network model of a neurotransmitter signaling pathway [17].

Consider deriving the simplified representation. One of the minimum feedback vertex sets is given by \({\mathcal W} = \{ 1,2,10,13 \}\). Then, under Assumption 1, we can obtain the following simplified representation:

\[\begin{aligned} \hat{x}_1(k+1) &= \hat{x}_1(k), \\ \hat{x}_2(k+1) &= \hat{x}_2(k), \\ \hat{x}_3(k+1) &= \hat{x}_2(k), \\ \hat{x}_4(k+1) &= \hat{x}_2(k), \\ \hat{x}_5(k+1) &= \hat{x}_2(k), \\ E[\hat{x}_6(k+1)] &= 0.3 \hat{x}_2(k), \\ E[\hat{x}_7(k+1)] &= 0.3 \hat{x}_2(k), \\ \hat{x}_8(k+1) &= \hat{x}_2(k), \\ \hat{x}_9(k+1) &= \hat{x}_2(k), \\ \hat{x}_{10}(k+1) &= \hat{x}_{10}(k), \\ E[\hat{x}_{11}(k+1)] &= 0.3 \hat{x}_2(k) - 0.3 \hat{x}_2(k) \hat{x}_{10}(k), \\ E[\hat{x}_{12}(k+1)] &= 1 - 0.3 \hat{x}_2(k) + 0.3 \hat{x}_2(k) \hat{x}_{10}(k), \\ \hat{x}_{13}(k+1) &= \hat{x}_{13}(k), \\ E[\hat{x}_{14}(k+1)] &= 0.09 \hat{x}_2(k) \hat{x}_{13}(k), \\ & \ \ \ \ - 0.09 \hat{x}_2(k) \hat{x}_{10}(k) \hat{x}_{13}(k), \\ \hat{x}_{15}(k+1) &= \hat{x}_2(k), \\ \hat{x}_{16}(k+1) &= \hat{x}_2(k). \end{aligned}\]

We remark that the dynamics of some nodes except for nodes in \({\mathcal W}\) are deterministic. From this representation, the following two facts for the original system can be obtained.

1. By setting \(x_1(0)=x_2(0)=x_{10}(0)=x_{13}(0)=0\), \(x_i\), \(i \neq 12\) converge to \(0\), but only \(x_{12}\) converges to \(1\).
2. If it is desirable that \(x_{12}\) converges to a smaller value, then we must set \(x_2(0)=1\) and \(x_{10}(0)=0\).

We validate these facts by a numerical simulation. We consider the following two cases.

• Case 1:  \(x_i(0)=0\), \(i \in \{ 1,2,10,13 \}\) and \(x_j(0)=1\), \(j \neq i\).
• Case 2:  \(x_{10}(0)=0\) and \(x_j(0)=1\), \(j \neq 10\).

Figure 7 shows the simulation result in Case 1. From this figure, we see that \(x_{i}\), \(i \neq 12\) converge to \(0\), and \(x_{12}\) converges to \(1\). That is, to stabilize this system at the origin, \(x_{12}\) must be also controlled directly. From each sample, we see also that all \(x_i\) converge to either \(0\) or \(1\) with the probability \(1\). Figure 8 shows the simulation result in Case 2. From this figure, we see that the average of \(x_{12}\) converges to \(0.7\). This fact can also be obtained from the simplified representation. Thus, using the simplified representation, we can obtain useful information for controlling the long-time behavior of a given PBN.

Fig. 7  Time response of the average of 10000 samples for the PBN in Sect. 6 (Case 1). Green line: \(x_{12}\).

Fig. 8  Time response of the average of 10000 samples for the PBN in Sect. 6 (Case 2). Green line: \(x_{12}\).

7.  Conclusion

In this paper, we proposed a simplified representation based on the steady value of the expected value of the state. In the proposed simplified representation, the behavior of a given PBN is characterized by some nodes determined by a minimum feedback vertex set of an interaction graph. In numerical simulations, we used the BN model of a neurotransmitter signaling pathway. Through numerical simulations, we clarified that the simplified representation with four nodes can be obtained from the PBN with sixteen nodes. In addition, the obtained representation was applied to the input selection problem. By the proposed method, the behavior of complex network systems can be represented by a mathematical model with a smaller number of nodes. Hence, it is expected that we can consider the analysis/control problems for large-scale systems such that handling has been difficult.

Input selection was discussed based on the use of only a constant control input. In complex network systems such as gene regulatory networks, such an input is frequently effective, because there is the case where frequent switching a binary value of the control input is difficult. On the other hand, to improve the behavior, there is a case that switching a binary value of the control input depending on the state is important. One of the future efforts is to develop a derivation method of state-feedback controllers based on the proposed method.

Acknowledgments

This work was partly supported by JSPS KAKENHI Grant Numbers JP21H04558, JP22K04163, JP23H01430.

References

[1] S.A. Kauffman, “Metabolic stability and epigenesis in randomly constructed genetic nets,” Journal of Theoretical Biology, vol.22, pp.437-467, 1969.

[2] I. Shmulevich, E.R. Dougherty, S. Kim, and W. Zhang, “Probabilistic Boolean networks: A rule-based uncertainty model for gene regulatory networks,” Bioinformatics, vol.18, no.2, pp.261-274, 2002.

[3] K. Kobayashi and K. Hiraishi, “An integer programming approach to optimal control problems in context-sensitive probabilistic Boolean networks,” Automatica, vol.47, no.6, pp.1260-1264, 2011.

[4] K. Kobayashi and K. Hiraishi, “Optimal control of probabilistic Boolean networks using polynomial optimization,” IEICE Trans. Fundamentals, vol.E95-A, no.9, pp.1512-1517, Sept. 2012.

[5] R. Li, M. Yang, and T. Chu, “State feedback stabilization for probabilistic Boolean networks,” Automatica, vol.50, no.4, pp.1272-1278, 2014.

[6] H. Li, Y. Wang, and P. Guo, “State feedback based output tracking control of probabilistic Boolean networks,” Information Sciences, vols.349-380, pp.1-11, 2016.

[7] R. Pal, A. Datta, M.L. Bittner, and E.R. Dougherty, “Intervention in context-sensitive probabilistic Boolean networks,” Bioinformatics, vol.21, pp.1211-1218, 2005.

[8] R. Pal, A. Datta, M.L. Bittner, and E.R. Dougherty, “Optimal infinite-horizon control for probabilistic Boolean networks,” IEEE Trans. Signal Process., vol.54, no.6, pp.2375-2387, 2006.

[9] I. Shmulevich and E.R. Dougherty, Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks, Society for Industrial and Applied Mathematics, 2010.

[10] P. Trairatphisan, A. Mizera, J. Pang, A.A. Tantar, J. Schneider, and T. Sauter, “Recent development and biomedical applications of probabilistic Boolean networks,” Cell Commun. Signal., vol.11, no.46, 25 pages, 2013.

[11] S. Zhu, J. Lu, Y. Liu, T. Huang, and J. Kurths, “Output tracking of probabilistic Boolean networks by output feedback control,” Information Sciences, vol.483, pp.96-105, 2019.

[12] A. Veliz-Cuba, “Reduction of Boolean network models,” Journal of Theoretical Biology, vol.289, pp.167-172, 2011.

[13] K. Kobayashi, “Design of fixed points in Boolean networks using feedback vertex sets and model reduction,” Complexity, vol.2019, article ID 9261793, 9 pages, 2019.

[14] T. Akutsu, S. Kuhara, O. Maruyama, and S. Miyano, “A system for identifying genetic networks from gene expression patterns produced by gene disruptions and overexpressions,” Genome Informatics, vol.9, pp.151-160, 1998.

[15] B. Fiedler, A. Mochizuki, G. Kurosawa, and D. Saito, “Dynamics and control at feedback vertex sets. I: Informative and determining nodes in regulatory networks,” J. Dyn. Diff. Equat., vol.25, no.3, pp.563-604, 2013.

[16] A. Mochizuki, B. Fiedler, G. Kurosawa, and D. Saito, “Dynamics and control at feedback vertex sets. II: A faithful monitor to determine the diversity of molecular activities in regulatory networks,” Journal of Theoretical Biology, vol.335, pp.130-146, 2013.

[17] S. Gupta, S.S. Bisht, R. Kukreti, S. Jain, and S.K. Brahmachari, “Boolean network analysis of a neurotransmitter signaling pathway,” Journal of Theoretical Biology, vol.244, pp.469-469, 2007.

[18] G. Even, J. Naor, B. Schieber, and M. Suden, “Approximating minimum feedback sets and multicuts in directed graphs,” Algorithmica, vol.20, no.2, pp.151-174, 1998.

[19] R.M. Karp, “Reducibility among combinatorial problems,” Proc Symp. on Complexity of Computer Computations, pp.85-103, 1972.

[20] L. Tournier and M. Chaves, “Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics,” Journal of Theoretical Biology, vol.260, no.2, pp.196-209, 2009.

[21] H.P. Williams, Model Building in Mathematical Programming, 5th ed., Wiley, 2013.

Authors

Koichi KOBAYASHI
  Hokkaido University

Keyword

feedback vertex set,  polynomial representation,  probabilistic Boolean network,  simplified representation