Polynomial-Time Algorithm for Controllability Test of a Class of Boolean Biological Networks

In recent years, Boolean-network-model-based approaches to dynamical analysis of complex biological networks such as gene regulatory networks have been extensively studied. One of the fundamental problems in control theory of such networks is the problem of determining whether a given substance quantity can be arbitrarily controlled by operating the other substance quantities, which we call the controllability problem. This paper proposes a polynomial-time algorithm for solving this problem. Although the algorithm is based on a sufficient condition for controllability, it is easily computable for a wider class of large-scale biological networks compared with the existing approaches. A key to this success in our approach is to give up computing Boolean operations in a rigorous way and to exploit an adjacency matrix of a directed graph induced by a Boolean network. By applying the proposed approach to a neurotransmitter signaling pathway, it is shown that it is effective.


Introduction
Various approaches to modeling, analysis, and control synthesis of biological networks such as gene regulatory networks and metabolic networks have been recently developed in the control community as well as the theoretical biology community [1]. In these approaches, it is one of the final goals to develop systematic drug discovery and cancer treatment [2,3]. Biological networks in general can be expressed by ordinary/partial differential equations with high nonlinearity and high dimensionality. Since such complexities cause difficulties in analysis and control design, various simpler models such as Petri nets, Bayesian networks, Boolean networks, and hybrid systems have been proposed for dealing with complex and large-scale biological networks at the expense of rigorous analysis (see e.g., [4,5]). This paper discusses the controllability problem of biological networks. In gene regulatory networks, for example, the controllability problem is defined as the problem of determining whether expressions of genes of interest can be arbitrarily controlled by expressions of a specified set of the other genes. As far as we know, two approaches to the controllability analysis of such biological networks have been developed so far: a piecewise affine model-based approach and a Boolean network model-based approach. However, the former approach can be applied to only the class of relatively low-dimensional systems [6,7].
On the other hand, a Boolean network model, where binary state variables are assigned to nodes and the transition rules of the state are given by Boolean functions [8,9], will be more practical for analysis of large-scale biological networks thanks to its bold simplification. Akutsu et al. have recently discussed the controllability problem of Boolean networks with control nodes and controlled nodes and have proven that this problem is NP-hard in a general setting [10]. Furthermore, they have proposed a polynomialtime algorithm for the classes of networks including a tree structure or at most one loop, and an exponential-time algorithm for the other classes. Indeed there is a criticism that a Boolean network model is too simple as a model of biological networks, but for large-scale networks it will be able to provide some indication or clue towards further detailed analysis. Thus various approaches based on this model have been well-studied so far (see e.g., [11][12][13][14][15][16][17][18][19]).

EURASIP Journal on Bioinformatics and Systems Biology
Motivated by the theoretical results in [10], this paper also focuses on the controllability problem of Boolean networks with control nodes and controlled nodes and proposes a sufficient condition for the Boolean network to be controllable, which can be easily verified by a polynomialtime algorithm. Our standing point is to give up computing complex Boolean operations in a rigorous way and to focus on deriving an easily-checkable sufficient condition for controllability so as to be applied to large-scale networks. The obtained algorithm is based on simple operations on an adjacency matrix of a directed graph induced by a Boolean network. This is a remarkable point of our approach, different from the method in [10], and enables us to apply our approach to a wider class of Boolean networks including nontree structures.
First, after the definition of controllability of Boolean network models with control nodes and controlled nodes is described, a sufficient condition for the controllability is derived in the form of an algorithm. Next, the computational complexity for the algorithm is discussed to show that it is a polynomial-time algorithm. In addition, PC-based numerical experiments show that the obtained algorithm is applicable to a class of Boolean networks with at least 1000 nodes. Finally, as an illustrative example, the proposed algorithm is applied to the Boolean network model of a neurotransmitter signaling pathway [20], which expresses an interaction pathway between the glutamatergic and dopaminergic receptors. Note that the polynomial-time algorithm proposed in [10] cannot be always applied to this problem. This Boolean network model consists of 16 nodes, and the problem of simultaneously controlling two important nodes among them, that is, concentration of exocytosis and phospholipase C, is discussed based on the proposed algorithm. As a result, we show that for example, they can be simultaneously controlled by keeping substance concentration at the other 4 nodes constant with appropriate values. Notation 1. Let N denote the set of nonnegative integers and {0, 1} m×n the set of m × n matrices consisting of elements 0 and 1. We also denote by I n and 0 m×n the n × n identity matrix and the m × n zero matrix, respectively. For simplicity of notation, we sometimes use the symbol 0 instead of 0 m×n and the symbol I instead of I n . Let M express the transpose of the matrix M.

Boolean Network Models
This section provides a brief review on a Boolean network model [8,9]. A Boolean network model consists of a set of nodes and a set of regulation rules for nodes, where each node expresses a gene, a molecule, or an event in the genetic network. The state variable ξ i at node i takes a Boolean value of 0 or 1 representing "inactive" or "active" status of the node, respectively. A regulation rule for each node is given in terms of a Boolean function, and each node state changes synchronously.
As an example, we consider a very simple and interesting Boolean network model of an apoptosis network in Figure 1  given by where ¬, ∧, and ∨ denote logical NOT, AND, and OR, respectively, k ∈ N denotes the discrete time, the concentration level (high or low) of the tumor necrosis factor (TNF, a stimulus) is denoted by ξ 1 , the concentration level of the inhibitor of apoptosis proteins (IAP) by ξ 2 , the concentration level of the active caspase 3 (C3a) by ξ 3 , and the concentration level of the active caspase 8 (C8a) by ξ 4 .
Here if the binary variable ξ i has the value of "1", then the concentration of a certain reactant gets larger than a prescribed threshold (i.e., it is active), otherwise less than that. In addition, logical NOT corresponds to inhibition of gene expressions. Since the caspase C3a is responsible for cleaving or breaking many other proteins, a high-level of the C3a concentration, that is, ξ 3 = 1 implies cell near-death; otherwise, cell survival. As seen in (1), if the concentration of IAP is high (ξ 2 = 1) or the concentration of the caspase C8a is low (ξ 4 = 0), then the concentration of C3a gets low, that is, ξ 3 = 0. On the other hand, ξ 2 and ξ 4 at the next time depend on the value of ξ 3 as well as ξ 1 . In this way, some dynamical interactions exist. See [21,22] for further details.

Problem Formulation
In a Boolean network model (2), the state ξ(k) is uniquely determined by giving the initial state ξ(0) = ξ 0 ∈ {0, 1} l , which implies that (2) is an autonomous system and has no control nodes.
On the other hand, this paper will consider the Boolean network model with control (i.e., input) nodes and controlled (i.e., output) nodes to discuss the outputcontrollability of this model. This model is given by where each element of u ∈ {0, 1} m denotes the state of the control node whose value can be arbitrarily given as an external control input in the Boolean network, each element of x ∈ {0, 1} n denotes the state of the node except for the control nodes in the Boolean network, and each element of y ∈ {0, 1} r denotes the state of the node to be controlled as an output in the network. Note here that y does not imply a measured output. Hereafter according to control theory, x, u, and y are called a "state", "control input" and "output", respectively. In addition, f : {0, 1} n × {0, 1} m → {0, 1} n is a Boolean function, and C ∈ {0, 1} r×n is the output matrix satisfying for each element Furthermore, the product of C and x in y = Cx expresses a product operation on matrices/vectors of the real number field. Thus the above condition on C guarantees that the output is the state variable itself, that is, for each i there exists j such that y i = x j holds. The case of y = x is also included here. This condition on C will not be restrictive in analyzing controllability of biological networks such as gene regulatory networks, since the relation on regulation among genes/molecules will be mainly discussed there.
For the system Σ of (3), the notion of output-controllability is defined as follows.
Definition 1. Suppose that for the system Σ of (3), the finite time T ∈ N and the initial state x(0) = x 0 ∈ {0, 1} n are given. Then the system Σ is said to be T-output-controllable at x 0 if for every y f ∈ {0, 1} r , there exists a control input sequence u(k) ∈ {0, 1} m , k = 0, 1, . . . , T − 1, such that y(T) = y f . Furthermore, the system Σ is said to be T-outputcontrollable if it is T-output-controllable at every x 0 .
The above notion of controllability comes from the fact that, for example, in control of genetic networks we often would like to determine if expressions of certain gene of interest (corresponding to y) will be able to be inhibited (or activated) by means of appropriately adjusting the expressions of a given set of genes (corresponding to u). It is remarked that we assume that the control time T is explicitly specified in the above definition.
Let us get back to the Boolean network model (1) of an apoptosis network. As discussed in [21,22], we consider ξ 1 (TNF) itself as a control input. So by ignoring the dynamics on ξ 1 , that is, ξ 1 (k + 1) = ξ 1 (k), we suppose in (1) that x(k) = [ξ 2 (k) ξ 3 (k) ξ 4 (k)] and u(k) = ξ 1 (k), which yields (3) of the form where x i (k) denotes the i-th element of x(k). As for the output y = Cx, either case of can be treated by assumption. Then let us verify the T-output controllability of the system (5). As discussed in Section 2, x 2 (= ξ 3 ) = 1 expresses cell near-death, and x 2 (= ξ 3 ) = 0 expresses cell survival. So we would like to know if the system is T-output-controllable with respect to the output y = x 2 .
Suppose that x 0 = [0 0 0] (i.e., the initial states of IAP, C3a and C8a are all low-level), C = [0 1 0] (i.e., y = x 2 ), and T = 2. Then since y(2) = 0 holds independently of u by simple calculation, we see that system (5) is not 2-outputcontrollable at x 0 , which implies that we cannot control the system from the state "cell survival" within 2 time steps no matter how the control value of u is given. On the other hand, suppose in (1) that x(k) = [ξ 1 (k) ξ 2 (k) ξ 3 (k)] and u(k) = ξ 4 (k). Then we obtain (3) of the form where ξ 4 This implies we can simultaneously control the value of x 2 and x 3 at T = 2.
In this way, the proposed controllability enables us to verify the existence of a control input sequence such that the output has the desired value in a given finite time, and the obtained result indicates how to give the value of a control input sequence.

EURASIP Journal on Bioinformatics and Systems Biology
Next, we will explain our basic strategy for deriving the controllability condition. Let us consider a Boolean network expressed as the state equation which is given by [10]. Although this model is very simple, it provides significant clues to address this problem. For the Boolean network model (9), we can consider three possible specifications, choosing either ξ 1 (k), ξ 2 (k), or ξ 3 (k) to be the control input for the system.
Note that for (11) with y = x, we see that the controllability property does not hold due to the fact that y(T) directly depends on x(0). On the other hand, for (13) with y = x, y 1 (2)(= x 1 (2)) is adjacent to u(1) and u(0) in the Boolean network, which implies that y 1 (2) is arbitrarily given by u(1) and u(0). In a similar way, y 2 (2)(= x 2 (2)) is adjacent to u(1). However, (y 1 (2), y 2 (2)) = (1, 1) cannot be realized by u(0) and u(1) because y 1 (2) = 0 always holds when y 2 (2) = 1. These examples are very important in discussing the controllability in a Boolean network, that is, if the Boolean function of y i (T) includes an initial state x(0), or includes the same input in the outputs at the same time, then the system in question is not T-output-controllable. In the following section, by motivating the above discussion, we will consider to derive a controllability condition. Remark 1. In the above example, we assume that when some genes are identified as control inputs, the original dynamics of the corresponding genes can be ignored. However, in the case that the corresponding gene has a strong interaction with other genes, this assumption may not be suitable. One of methods for coping with such a case is to add a new gene (node) that works as the control input [10], where it is called an external control node. Our approach below can be also applied to this case.

Preliminaries.
This section presents a sufficient condition for the system (3) to be T-output-controllable in the form of an algorithm.
Consider a simple example given by This system has the following relation: Similarly, we see that y 2 (T) = 0, T ≥ 2, hold identically. In Boolean functions, identical equations are in general given by where h(·) is any Boolean function of a vector of binary variables. Obviously such identities on x i or u i affect the controllability in a Boolean network (note that even if Let us consider again the Boolean network model (5) of an apoptosis network. If we suppose that x(0) = x 0 = [0 0 0] , C = [0 1 0] (i.e., y = x 2 ), and T = 2, then by a simple calculation, we obtain the following identity: So in Boolean biological networks, there exists the case that identities are appeared. However, identities may not be appeared in the real biological relevance. The reasons why such identities are appeared are that the state is binarized and that a time-delay of the state is ignored. To overcome the latter point, a temporal Boolean network model ξ(k + 1) = f a (ξ(k), ξ(k − 1), . . . , ξ(k − T)) has been proposed in [23]. However, identities may appear even in a temporal Boolean network. The output-controllability condition proposed below can be similarly applied to a temporal Boolean network model.
Thus first of all, we will focus on finding such identities in y(T) before discussing a kind of initial condition and a kind of input-independency. This will require the introduction for several symbols.
The following assumption is made.  h(a, b) can be rewritten as h(a) = a. Any given Boolean function can be changed so as to satisfy Assumption 1: after it is transformed into an appropriate canonical form (e.g., Reed-Muller canonical form (polynomials over the finite field GF(2))), it is easy to eliminate redundant variables by expanding based on four operations over GF (2). Also in the identification of Boolean network models (e.g., see [24]), since the correlations between variables are checked, the Boolean function f in (3) will satisfy Assumption 1 in many cases. By Assumption 1, it is guaranteed that the Boolean function f itself does not include any identities, although y(T) may include some identities. Let p denote the number of the logical NOT appeared in (3), where the logical NOT operators are distinguished when the corresponding terms are different even if the corresponding variables are the same. In addition, consider the fictitious inputs v i (k) = 1, i = 1, 2, . . . , p, which have one-to-one correspondence with the variables operated by the logical NOT, that is, ¬x i or ¬u i in (3). Then the system (3) can be equivalently rewritten as the following system: where the Boolean function f v does not include the logical NOT, and For example, system (17) is rewritten as subject to v(k) = 1.
Next, consider the adjacency matrix Φ ∈ {0, 1} (n+m+p)×(n+m+p) for the directed graph induced by the Boolean network of the system (21). For example, the adjacency matrix for the system (23) is given by where if there exists an arc from node i to node j, then the (i, j)-th element of Φ is 1. Hereafter, without loss of generality, the i-th element of [x u v ] is assigned to node i in the directed graph, where i ∈ {1, 2, . . . , n + m + p}.
In the case of (24), x 1 , x 2 , x 3 , u, and v are assigned to nodes 1, 2, 3, 4, and 5, respectively. Then in Figure 2, which shows a temporal/spatial network of the system (17), we say that for example, there exists a path between x 2 (2) and u(0). Using the adjacency matrix Φ, we also compute the matrix Φ t C 0 C , t = 1, 2, . . . , T, where In the case of the system (17), we have where C = [0 2×1 I 2 ].

EURASIP Journal on Bioinformatics and Systems Biology
For the system (21), Φ t C 0 C expresses whether there exist paths between y(T) and x(T − t), y(T) and u(T − t), or y(T) and v(T − t) for any given T. In the case of (26), we see that y 2 (2)(= x 3 (2)) is adjacent to x 1 (1) and x 2 (1). In other words, Φ t C 0 C expresses which elements of x(T − t), u(T − t), and v(T − t) are variables of a Boolean function representing y i (T). However, note here that from Φ t C 0 C , we cannot specify an explicit form of the Boolean function in question.
Furthermore, the following symbol is used: where X t ∈ N n×r , U t ∈ N m×r , and V t ∈ N p×r . Let also X t ix, jx , U t iu, ju , and V t iv, jv denote each element of X t , U t , V t , respectively. If X t ix, jx ≥ 1 holds, then there exist X t ix, jx paths between y jx (T) and x ix (T − t). For the state x ix , i x = 1, 2, . . . , n, of the system (3), let P x express the index set of elements of x ix operated by the logical NOT as ¬x ix . In a similar way, for the control input u iu , i u = 1, 2, . . . , m, of the system (3), let P u express the index set of elements of u iu operated by the logical NOT as ¬u iu . Here, p = |P x | + |P u | holds. In addition, there is a one-to-one correspondence between each element of P x , P u and the index i v of v. Let ν(i x ) and ν(i u ) express the index i v of v corresponding to i x ∈ P x and i u ∈ P u , respectively. In the case of the system (23), P x = ∅, P u = {1} hold, and for i u = 1, ν(i u ) = 1 holds. Finally, we define the following matrices: where

Proposed
Algorithm. Now we are in a position to propose a T-output-controllability test algorithm. Since this kind of problem is NP-hard [10], we pay our attention on deriving a sufficient condition for the controllability. Although this sufficient condition is given in the form of an algorithm, it is somewhat complex. Thus before describing an algorithm, we describe the outline of the algorithm. First, we consider a necessary condition for y(T) to include identical equations. From Figure 2 of the example (23), we see that y 2 (2)(= x 3 (2)) in (18), which has no identities, has two paths from u(0), and that v(0) is Figure 2: Temporal/spatial network of the system (23).
connected to some node on the paths. In this way, if some identical equation exists in y j (T), there always exist more than 2 paths from y j (T) to some state and also the logical-NOT operations exist on the paths, which is a necessary condition and not necessarily a sufficient condition. Since it will spend huge time to rigorously specify the existence of identities for a large network, we consider here to exclude the cases satisfying the above necessary condition, that is, we do not determine here the controllability in such cases. Next, for the system that includes no identical equations, we use a kind of input-independency to determine the controllability. For example, consider the case that neither identity on u nor x exists in y(T) and that y(T) is expressed by y 1 (T) = h 1 (u 1 (0), u 2 (3)), as a result of recursive calculation (see Section 6 for such an example), where h 1 , h 2 are some Boolean functions. This system is obviously T-output controllable because each y j (T) is expressed by different u i (k) and no x 0 exists in y j (T). From the viewpoint of adjacency relation, this implies that there exists no path between x(0) and y(T), there exists at least one path from each y j (T) to some u i (k), and each u i (k) has a path with only one y j (T) or has no path to any y j (T). This can be easily found from the adjacency matrix, although it is a sufficient condition for the controllability. This is a rough story of our approach. The proposed algorithm is given as follows.

Algorithm 1 (T-output-controllability test algorithm).
Part A: Check of the Existence of Identical Equations.
Step 2. If T = 1, go to Step 6. Otherwise set t = t + 1. Compute X t , U t , and V t .
Step 3. If there exists (i x , j x ) such that X t ix, jx ≥ 2 or (i u , j u ) such that U t iu, ju ≥ 2, denote them by (i * x , j * x ) or (i * u , j * u ), respectively, and go to Step 4. Otherwise, go to Step 2 if t < T and go to Step 6 if t = T.
Step 4. If there exists i * x such that i * x ∈ P x or i * u such that i * u ∈ P u , and V t holds, go to Step 8. Otherwise, go to Step 5.
Substep 5.2. If any element of j * x -th column or j * u -th column in V j is greater than or equal to 1, go to Step 8. Otherwise, go to Substep 5.3.

Part B: Check of the Independence of Each y(T).
Step 6. If the following conditions hold for the matrices X 0 and U in (28), system (3) is T-output-controllable, or else if only condition (i) does not hold, then go to Step 7. Otherwise go to Step 8.
(i) X 0 = 0 n×r holds; (ii) each column vector of U is a nonzero vector; (iii) each row vector of U is a zero vector, or has only one element with a nonzero value.
Step 8. This algorithm cannot determine whether the system (3) is T-output-controllable or not (at x 0 ).
The above algorithm allows us to determine the Toutput-controllability of the system as follows.
First, noting that the identical equations have the form in (19), and x(T) is obtained recursively from (21), we see that the identical equations appeared in x(T) always have the form where w(k) denotes either variable of x(k) or u(k), and I 1 , I 2 are some subsets of the index set {(i, j)|i = 1, 2, . . . , p; j = 0, 1, . . . , T − 1}. Then using the forms of (31) and (32), the following lemma on Part A of Algorithm 1 is obtained.
Proof. In Step 3, from X t ix, jx ≥ 2 for some t, i x = i * x , and j x = j * x , we see that more than 2 paths from y jx (T) to x ix (T − t) exist, which is necessary for the identity on x ix (T − t) to exist (similarly for the case of U t iu, ju ≥ 2). Thus we next focus on the existence of logical NOT (i.e., v i ) in these paths in Step 4 and Step 5.
Consider the case that the logical NOT (i.e., v i ) corresponding to Step 3 exists in (21), in other words, either i * x ∈ P x or i * u ∈ P u holds. Then the condition V t is included in the paths in question, which is a necessary condition for the existence of the identity in y j * x (T). Thus we exclude this case (Step 4). (similarly for the case u i * u (T − t)).
In the other case, from (31), (32), for v(T − j), some j ∈ {1, 2, . . . , t − 1}, to exist in the paths in question is necessary for the existence of identities. If any element of the j * x -column or the j * u -column of V j is greater than or equal to 1, some element of v(T − j) exists in the paths in question. Thus we exclude this case (Substep 5.2). Therefore, it follows that y(T) includes no identities in Step 6.
From Lemma 1, we see that the case that y(T) includes the identities that have the form of (31) or (32) is excluded from the viewpoint of a necessary condition for the identity to exist in y(T). Thus we obtain the following theorem. Theorem 1. For a given T, the following statements hold.

(i) the system (3) is T-output-controllable if conditions (i), (ii), and (iii) in Step 6 hold subject to Part A of Algorithm 1,
(ii) for a given x 0 ∈ {0, 1} n , the system (3) is T-outputcontrollable at x 0 if condition (iv) in Step 7 holds subject to Part A and Step 6.
Proof. First, the statement (i) is proven for the system satisfying the condition that y(T) includes neither identities of (31) nor identities of (32). From Lemma 1, this condition is satisfied in Step 6. Then condition (i) in Step 6 implies that there exists no path between each element of x(0) and each element of y(T), since the (i, j)-th element of X 0 expresses if a path from x i (0) to y j (T) exists or not. On the other hand, note that (mh + i, j)-th element of U expresses if a path from u i (T − h − 1) to y j (T) exists or not (h = 0, 1, . . . , T − 1). Thus condition (ii) in Step 6 implies that there exists at least one path from each element of y(T) to some u i (k). Furthermore, condition (iii) in Step 6 means that the input u i (k) for each i ∈ {1, 2, . . . , m} and k ∈ {0, 1, . . . , T−1} has a path connected to only one element of y(T) or has no path to any element of y(T). From these conditions, it follows that each u i (k) affects at most one y j (T) and not the other y h (T), h / = j. Hence the value of y j (T) can be 8 EURASIP Journal on Bioinformatics and Systems Biology independently specified by the corresponding u i (k), which implies that system (21) is T-output-controllable. Next, the statement (ii) is proven. Since condition (i) in Step 6 does not hold, in this case, there exists a path between some element of x(0) and some element of y(T). On the other hand, condition (iv) in Step 7 guarantees that there exists no path between constant elements of x(1) = f v (x 0 , u(0), v(0)) and elements of y(T). Thus y(T) is not affected by the value of x 0 . Therefore, from (ii)-(iv), it follows that system (21) is T-output-controllable at x 0 . This completes the proof.
As for identical equations, the proposed algorithm excludes the case of ¬h(a) ∧ ¬h(a) as well as (19). This is a weak point of this algorithm. Furthermore, consider the following system: This system is T-output-controllable for T = 1. However, the proposed algorithm cannot determine whether this system is 1-output-controllable or not; thus there exists a class of systems such that the proposed algorithm cannot determine the controllability. Needless to say, it will not be so easy to cope with various cases stated above due to high nonlinearity of Boolean functions.
While the proposed algorithm includes such disadvantages, one of the main advantages of the algorithm is that the computational complexity of the above algorithm is very small. This will be discussed in the following section.

Computational Complexity Analysis
In this section, we discuss the computational complexity of the algorithm proposed in the previous section.
First, let us recall the definition of the symbols used here. The number of the state, the control input, and the output in (3) are denoted by n, m, and r, respectively. The number of the logical NOT appeared in (3) is expressed by p. In addition, T ∈ N expresses the control time. Then the following result is obtained.

Lemma 2. The computational complexity of the proposed algorithm is O(
Proof. The computation of the proposed algorithm consists of (a) checking each condition of Part A, and (b) checking whether conditions (i) to (iv) hold or not.
First, (b) is considered. The computational complexity to compute Φ 2 is O((n + m + p) 3 ). So the computational complexity to compute Φ T and Φ is given by both O((n + m + p) 3 (T − 1)). Further, the computational complexity to compute the product of Φ and C 0 C is O((n + m + p)nrT). So by simple calculation, the computational complexity of U is obtained as O((n + m + p) 3 (T − 1) + (n + m + p)nrT). The computational complexity of generating X 0 is obviously less than the case of U. Therefore, the computational complexity to compute X 0 and U is O((n + m + p) 3 (T −1)+(n+m+ p)nrT), which also includes the computational complexity to check conditions (i) to (iv) in Steps 6 and 7 for given X 0 and U.
Next, (a) is considered. The matrices X t , U t , V t are obtained directly from ΦC 0 C , and the computational complexity of Step 5 is O(prT). As a result, since the computational complexity of each checking in Part A is O((n + m + p) 2 (T − 1)) + O(prT), the computational complexity of Part A is less than O((n + m + p) 3 (T − 1) + (n + m + p)nrT).
Therefore, the computational complexity of the proposed algorithm is given by O((n + m + p) 3

(T−1)+(n+m+ p)nrT).
From Lemma 2, we see that the proposed algorithm is a polynomial-time algorithm. Furthermore, the computational time for performing the proposed algorithm is evaluated by numerical experiments, where the total computational time in Part B is measured because from the proof of Lemma 2 we see that the computational complexity of Part B is dominant. So the adjacency matrices to be evaluated are generated randomly for each l(= n + m), where n = m = l/2, p = 0 are given. The results are shown in Table 1, where MATLAB on the computer with the Intel Core 2 Duo CPU 3.0 GHz and the 2 GB memory is used. In Table 1, the worst computational time implies the worst value among 100 cases randomly selected for each l. From Table 1, we see that the proposed algorithm can be applied to relatively large-scale Boolean network models. COMT  In the case of dim u(k) = 4 and T = 6, we can also find controllable control inputs among 14 C 4 (= 1001) input-combinations. For example, we obtain as one of combinations of x(k) ∈ {0, 1} 12 and u(k) ∈ {0, 1} 4 that make the system 6-output-controllable It is remarked that the polynomial-time algorithm proposed in [10] cannot be applied to the system with the state (38) and the input (39) because the network includes the two loops, that is, the loop of ξ 2 , ξ 4 , ξ 6 , and ξ 5 , and the loop of ξ 7 , ξ 11 , ξ 12 , and ξ 14 . Furthermore, based on the above result the Boolean function of y(6) can be derived as which implies that the value of y(6) can be freely given by control inputs, for example, (a) y(6) = [0 0] for u 1 (4) = 0, u 1 (1) = 0, u 2 (2) = 0, u 3 (0) = 1, u 4 (2) = 0, and (b) y(6) = [1 1] for u 1 (4) = 1, u 1 (1) = 1, u 2 (2) = 0, u 3 (0) = 1, u 4 (2) = 0. Finally, we discuss the control input sequence realizing the desired output values. One of criticisms in control of Boolean networks is to assume that the value of the control input can be arbitrarily given at each time. In many biological systems, this assumption is not always satisfied, and input constraints are frequently imposed. One of input constraints is that the value of the control input is given as a constant within a certain sufficiently long time period. Although it is one of future works to explicitly deal with such an input constraint, based on the proposed algorithm, we may also find a constant-valued sequence of control inputs for which the desired values of outputs are obtained. For example, in (40) and (41), let us consider to find a control input sequence satisfying y 1 (6) = 0 and y 2 (6) = 1. Since u 1 (4) = 0, u 1 (1) = 0, u 2 (2) = 1, u 3 (0) = 0, and u 4 (2) = 1 are obtained as one of solutions, it is remarked that the following control inputs are given as any binary value: u 1 (0), u 1 (2), u 1 (3), u 1 (5), and u 2 (0), u 2 (1), u 2 (3), u 2 (4), u 2 (5), and u 3 (1), u 3 (2), u 3 (3), u 3 (4), u 3 (5), and u 4 (0), u 4 (1), u 4 (3), u 4 (4), u 4 (5). This allows us to give the values of the control input sequences as a constant, that is, u 1 (k) = 0, u 2 (k) = 1, u 3 (k) = 0, u 4 (k) = 1, k = 0, 1, . . . , 5. Thus the proposed algorithm helps us to find a practically useful control input sequence. Furthermore, this kind of degree of freedom in control inputs may be used for the optimal control problem. Once we can determine control input variables by our algorithm, we can use a tool for finding optimal control input sequences, which have been developed in hybrid control theory, for example, [25,26]. It is expected that such an analysis will provide one of guidelines in experimental approaches to the control problem of biological networks.

Conclusion
In this paper, the controllability analysis for biological networks expressed by a Boolean network model with control nodes (inputs) and controlled nodes (outputs) has been discussed. First, a sufficient condition for the Boolean network model to be output-controllable has been derived by exploiting an adjacency matrix of its network graph. The obtained condition, which is given in the form of an algorithm, can be checked in polynomial time with respect to the state/input dimensions and the control time period; it will be one of the powerful tools that can provide some clues for finding effective control inputs to control a large-scale biological network. Next, by PC-based numerical experiments, it has been shown that the proposed method is applicable to large-scale Boolean networks with at least 1000 nodes. Finally, the proposed method has been applied to the Boolean network model expressing a neurotransmitter signaling pathway, and has shown that it is controllable with respect to both exocytosis and phospholipase C when appropriate control inputs are used.
There are many interesting open problems to be addressed in the future. It is one of the most important issues to characterize a class of Boolean networks to which our algorithm can be applied as a necessary and sufficient condition. In addition, extensions to the case of systems with input constraints and uncertainty are also one of the significant topics.