A Boolean network control algorithm guided by forward dynamic programming

Control problem in a biological system is the problem of finding an interventional policy for changing the state of the biological system from an undesirable state, e.g. disease, into a desirable healthy state. Boolean networks are utilized as a mathematical model for gene regulatory networks. This paper provides an algorithm to solve the control problem in Boolean networks. The proposed algorithm is implemented and applied on two biological systems: T-cell receptor network and Drosophila melanogaster network. Results show that the proposed algorithm works faster in solving the control problem over these networks, while having similar accuracy, in comparison to previous exact methods. Source code and a simple web service of the proposed algorithm is available at http://goliaei.ir/net-control/www/.

Detecting a set of perturbations which cause the desirable changes in cellular behavior has many applications such as cancer treatment and drug discovery [5,[19][20][21][22][23]. This highlights the necessity of developing a control theory for the gene regulatory networks. Using GRN control model is considered as a key method to design the experimental control policies [23].
The control problem includes finding a sequence of interventions to be applied on the system, which changes state of the system from an undesirable state of the network to a desirable one [24][25][26]. The undesirable state in a gene regulatory network may express a disease such as cancer, and the desirable state can express the wellness, for example as induction of apoptosis a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 large size networks. In this paper, we have presented a new algorithm to solve the control problem, and we investigate the applicability of the provided algorithm on two GRNs: T-cell receptor network [40] and Drosophila melanogaster network [41], and compare the efficiency with other algorithms.
In another completely different view, a mathematical formulation of the controllability problem of Boolean networks is provided. These works formulate the dynamics of Boolean networks as linear systems. They defined an extension of matrix multiplication to achieve this goal. They applied this technique on classical Boolean network control problem and showed that this formulation exactly models the desired problem. Le et al proved theorems that specify some necessary and sufficient conditions for a Boolean network to be controllable [42]. The same technique is applied to model a Boolean network with time delays. In a Boolean network with time delay, state of nodes may affect other nodes after some predefined delay. Some necessary and sufficient conditions for controllability of these networks is provided [43]. An extension to the network control problem, which has not only one initial and one desired states, but also a set of states as initial and desired ones (called trajectory), is provided recently [44]. With similar algebraic techniques, they proved theorems on controllability for trajectory problems.

General idea of proposed algorithm
The proposed algorithm is an exact algorithm, i.e. an algorithm that always finds the correct answer, but it is using some properties of the network to enumerate a smaller number of states and perform faster than previously provided algorithms. The general idea of the proposed algorithm is to divide the network to partially-dependent parts and then consider all the possible states the network may be in. To discover these independent subnetworks, we consider strongly connected components and categorize them to find independent parts. Strongly connected components induce an order of dependency between genes. Thus, we can handle them with the induced order, one by one. For each strongly connected component, we consider all possible states. If possible states of a subnetwork (that forms a strongly connected component), is independent of some other parts, we can enumerate these two parts independently, and by this trick, the number of states that we should consider would be reduced. More precisely, if there are k independent parts having number of status #S 1 , . . ., #S k , number of status that we consider reduces from #S 1 × . . . × #S k to #S 1 + . . . + #S k . The multiplication is the number of states that the Datta algorithm enumerates. With this trick, we reduce the running time of our algorithm.
As it could be seen, size of strongly connected components of a network could be a good estimation of how faster our algorithm is, in comparison of previously provided dynamic programming algorithms (e.g. Datta algorithm). In a recent study, GRN of E. coli is analyzed. In this network, among 1807 genes, there are 202 transcription factors. This network does not contain any strongly connected component with more than 5 nodes [45]. Also, a manually curated GRN of mouse consists of 274 genes with 176 transcription factors. Its larges strongly connected component contains only 80 genes, in which only 29 are transcription factors [45]. These observations show that real regulatory networks have small strongly connected components which are the cases that our algorithm works better for them. By taking advantages of this observation we provided an algorithm in this paper.

Background on Boolean network control
In a Boolean network, each node represents a gene, and each edge represents a regulatory effect of one gene expression on another one, which may cause an increase or decrease in the gene expression [46].
A Boolean network is modelled as a directed graph G = (V, F), which includes a set of n nodes V = {v 1 , v 2 , . . ., v n }, and a set of Boolean functions F = {f 1 , f 2 , . . ., f n }. Boolean network is a discrete time model. Each node v i has a state variable v i (t) 2 {0, 1} representing the state of node v i at time t, where 0 (1) indicates lack of (existence of) gene expression. Also, each node v i has a Boolean function f i , representing how to obtain v i (t + 1) from the state of the incoming nodes to v i at time step t by applying basic Boolean operations (and, or, not). The network state at time t, is defined as vector v t = [v 1 (t), v 2 (t), . . ., v n (t)], which describes the state of the nodes in time step t.
An example of a Boolean network is represented in Fig 1(a). In Fig 1( In the control problem of Boolean networks, a Boolean network G = (V, F), initial state and a desirable network state v τ is given. The set of nodes V is called internal nodes. A set of control nodes {u 1 , . . ., u m } are added to the network, known also as external nodes, which are used to influence internal nodes to attain the desirable state. The external nodes have no incoming edges, and their values are specified externally. The problem is to find a sequence of state values u 0 , . . ., u τ for external nodes, which leads the network to be in desirable state v τ in time step τ. If there exists no such control sequence, this fact should be announced as the output. In gene regulatory networks, the desirable state of the network represents a healthy state of the system, and external nodes represent potential medicines affecting network behavior. Thus, finding control strategies has applications in various medical areas, including medical protocol design for example in cancer treatment [19,20].
An example of a control problem on a Boolean network is represented in Fig 2. In this example, {v 1 , . . ., v 6 } is the set of internal nodes and {u 1 , u 2 , u 3 } is the set of external nodes. The initial state of the network is v 0 = [1, 0, 1, 1, 0, 0], and the desirable state is v 3 = [0, 1, 1, 1, 0, 1]. Thus, we are looking to find a control sequence u 0 , . . ., u 3 in such a way that the network would be in state

The proposed algorithm
The idea of our proposed algorithm is to reduce backtracking space by separating dependencies and merging information from some sub-parts of the graphs, which we call them components. In addition, we categorize components to find some easily solvable components, nonbranching and in some cases the branching components, to reduce the number of states to enumerate. In our proposed algorithm, we compute an intermediate variable U to be used in the computation of the final result. For each node, at each time step, we calculate answer of two questions: is it possible for this node to be in state 1 (0), and store result of this question in an intermediate variable U. For this computation, we design the following steps.
For network node v i , Boolean variable b 2 {0, 1}, and time step t, we define variable U v i b ðtÞ, which is true if and only if it is possible to assign values to external nodes in such a way that it cause v i (t) to get value b, and it is false otherwise. In the other words, U v i 1 ðtÞ and U v i 0 ðtÞ represent the possibility for node v i in time step t to have state value 1 and 0, respectively.
Let v i 1 . . . v i k be the input nodes to node v i . According to the definition of U we have U v i 1 ðt þ 1Þ ¼ true if, and only if, there exists ½b To obtain the final result, we should check the existence of a control sequence which leads to U v i v t ½i� ðtÞ ¼ true for all the nodes. Also, to specify and output the desired control sequence, the regression technique may be used [39].
Pseudocode of our proposed algorithm is presented in Fig 3. The source codes are available as S1 Codes.

Strongly connected components.
In order to compute U values, we partition the network into strongly connected components.
A strongly connected component (SCC) of a network, is a maximal subset of network nodes, in which every node is reachable from every other one. We partition the network nodes into strongly connected components, using SCC algorithm [47]. This algorithm is one of the classical graph theory algorithms based on two depth-first searches on graph and reverse of the graph, respectively. A topological order on strongly connected components of a network, is an order on its components such that for every directed edge xy from component x to component y, x appears before y in the ordering.
We find a topological order on the components using topological sort algorithm [47]. Then, based on the size of the component and the number of outputs edges, we categorize components into three categories: 1) non-branching single node components, 2) branching single node components, 3) multi-node components. A branching node is a node with at least two outgoing edges and a non-branching node is a node with at most one outgoing edge.

Non-branching single node components.
In this case, current component consists of one node v i , and v i has at most one outgoing edge. We simply find U v i 0 ðtÞ and U v i 1 ðtÞ for 0 � t � τ from state values of incoming nodes to v i in last time step.
For example, node v 2 in Fig 2 is a Non-Branching Single Node Component for which U v 2 0 ðtÞ and U v 2 1 ðtÞ for 0 � t � 2 should be calculated based on state values of their incoming nodes v 1 and u 2 in the previous time step. Note that, interestingly, for node v 2 in time step t = 2 it is possible to be in state 1 and state 0. Another example with four non-branching single node components and its U values is present in Fig 4. In accordance with dynamic programming arrays, we fill U to desirable time step τ. Since the nodes are met based on a topological sort, and all the states of the incoming nodes of the related node till the time step τ are calculated before, this calculation is possible. In time step τ, check if the related node has attained the desirable state or not. If the answer was negative, announce that there is no control sequence.

Branching single node components.
In this case, current component consists of one node v i , and v i has at least two outgoing edges. Same as the previous case, we can simply find values U v i 0 ðtÞ and U v i 1 ðtÞ. The difference is that since v i has more than two outgoing edges, there are at least two other nodes in the network which their state values depend on the state value of v i in the previous time. Thus, in some networks, the state value of v i which is required by its outgoing neighbors may be inconsistent. An example of a branching node that is not able to attain both states 0 and 1 in a time step is present in Fig 5. Suppose that initial state is 101100 (as it is the case in Fig 2) and we are supposed to reach to a state with v 6 = 0 and v 3 = 1 in time step t = 4. If we treat this node same as Non-Branching Single Node Components, we would face a problem. This algorithm announces that in time step t = 4, we can have v 6 = 0 and v 3 = 1, but there is no control sequence reaches to a state with v 6 = 0 and v 3 = 1 in t = 4. The problem is that both v 3 and v 5 in t = 3 are able to attain values 0 and 1, but, there is a dependency between them. This dependency does not let v 6 = 0 and v 3 = 1 in t = 4, simultaniously.
To deal with this problem, we fix some value for the node and recursively consider nodes for next SCCs. In this manner, we remove the dependency between the following nodes by considering different states one by one, recursively.
In other words, as we detect a dependency between states, we enumerate recursively all possible states. In the case of the branching node, if the branching node has more than one state in time t, we go through backtracking. Otherwise, we can branch-out some states. If the node has one possible state, there is no need to go through backtracking. We can simply fill the next components sequentially.

Multi-node components.
In this case, the current component consists of at least two nodes. We apply Datta algorithm in this case [24]. Datta et al fill an array , t], for t = τ to t = 0 according to the following procedure: if there exists ðv t ; uÞsuch that D½v 1 ðtÞ; v 2 ðtÞ; . . . ; v n ðtÞ; t� ¼ 1 and After applying Datta algorithm, we check if it is possible to attain the desirable state in time step t = τ. Then, we fill U array accordingly. For example, node v i in time step t may attain the desired state in time step τ with both states 0 and 1, thus, we will have U v i 0 ðtÞ ¼ true and U v i 1 ðtÞ ¼ true. Then we partition nodes with edges going out of the component into two categories, nodes with at most one outgoing edge and nodes with at least two outgoing edges (branching node). Nodes with at most one outgoing edge are easy to handle, however, we treat branching nodes which may attain both states 0 and 1 in the same time step, like branching single node components. In both cases, we fix values for this component and go for the next component recursively.

Complexity analysis.
Let C 1 , . . ., C #C be the strongly connected components of G, and σ t (C i ) be the number of valid states of component C i at time step t. Our algorithm enumerates all the network states, but non-branching single node components and branching single node components.
Datta algorithm finds all potentially desirable network states backward in time. To find all desirable states in time step t, it generates all potential states for the network, which takes O(2 n ) time, for an n node network. Then, for each potential state, it calculates next state and checks whether it leads to the desirable state by searching it within desirable states of the time step t + 1. Let T(G) be the time required to calculate the next state of a network's state. The amount of computation for time step t, for each potential state is T(G) for calculating the next step in addition to searching it, which takes lg σ t+1 (G). Datta performs these two operations for all 2 n potential states. Therefore, Datta algorithm's running time is O(∑ t 2 n (T(G) + lg σ t (G))) = O(τ2 n T(G)). Note that, calculation of T(G) requires computation of Boolean functions for each node, and if functions are given as input, their calculation asymptotic time is not more than the size of the input, which is acceptable.
On the other hand, our algorithm proceeds component by component. For each component, we consider all the σ t (C i ) states for all time steps t. However, in the following three cases, we do not enumerate all the states: • For non-branching single node components, we compress states. In other words, we check further components not by considering two possible states for these components recursively. However, since network states of following components do not depend on different values of these components, we make our enumerations by half, for each of these components. Thus, if we have k such components, our running time will be reduced by a factor of 1 2 k , in comparison to normal dynamic programming. • Leaf components with out-degree 0 are independent, and not required to be checked recursively, but sequentially.
Considering above notes, our running time is O(S t P i:I(t) σ t (C i )T(C i )). Note that since 2 n � σ t (G) � P i σ t (C i ) and T(G) � ∑ i T(C i ), our running time is always better than running time of dynamic programming.

Example
Consider the network shown in Fig 6. This network consists of one external node u and four internal nodes v 1 , v 2 , v 3 , and v 4 . Transition function and initial/desired state for the network is shown in Fig 6. Based on the proposed algorithm, we first partition network nodes to SCCs. SCCs of the network of First we process SCC S 1 = {u}. Since node u is an external node, in every time step t, u can attain both states 0 and 1, thus, U u b ðtÞ ¼ true for all 0 � t � 3 and b 2 {0, 1}. This SCC is a branching single node component. Thus, we fix U u b ðtÞ for all 0 � t � 3 for all different combinations of assignment values, as it is represented in Table 1. In the recursion, if correct solution is found, the algorithm reports the solution and halts. Otherwise, the algorithm fixes another combination and goes for next SCC, recursively. For this example, suppose that the algorithm fixes row #2 of the Table 1 for component S 1 , and then goes through next SCCs.
Next SCC, according to the topological order of SCCs, is S 2 = {v 1 , v 2 }. The logic behind behavior of S 2 is that, if u(t) = 1, then, we have v 1 (t + 1) = 1 and v 2 (t + 1) = 0, and if u = 0, v 1 and v 2 exchange their state at the next time step. S 2 is a multi-node component. Thus, the algorithm first executes a Datta algorithm on this SCC (line 28 of Fig 3). The result of Datta algorithm is shown in Table 2. As it is shown in Table 2, since we already fixed one state over all time steps for SCC S 1 , there is only one possible state for the S 2 . For SCC S 2 , after calculating U table, since S 2 is a multi-node component, the algorithm fixes this state and goes through next SCCs.
Next SCC, S 3 = {v 3 }, is a single node non-branching component. Thus, the algorithm calculates U v 3 b ðtÞ for all 1 � t � 3 and b 2 {0, 1}. By the definition of transition function (Fig 6), . Then, we proceed to next SCC. Next SCC is S 4 = {v 4 }. So far, the value of U table is shown in Table 3. Like SCC S 3 , SCC S 4 is a single node non-branching component. Thus, the algorithm first calculates U table. According to the transition function (Fig 6), . Now, the algorithm processed all the SCCs, thus, it reports the result.
Note that, in any of the recursion steps, if the algorithm does not find an answer that matches to the desired state, it traces back to the SCCs that are not single node non-branching, tries the next state from the state space. In the case of the current example, the algorithm will

Proof of the correctness of the proposed algorithm
In a general view, the proposed algorithm enumerates all possible states, but, removes some useless states. There are two extreme approaches to enumerate the states. One approach is to recursively enumerate them node-by-node. The other extreme, like Datta algorithm, is to fill an array that represents reachability of the states, time step-by-time step. Our algorithm is somewhat in-between, it fills arrays for SCCs, and fix them recursively.
We provide a stronger result for the proposed algorithm. We let logical formulas that describe states of the nodes in time step t + 1, to not only include variables of the time step t, but also variables of time steps 0 to t. We name these networks as extended networks. Also, we name single node non-branching SCCs as simple and the other SCCs as hard components.
Consider all the non-branching single node components that proceed an SCC S in the topological order. All easy components may have 0 or more input edges, but only one output edge. Note that, if there is a single node non-branching SCC without output edge, it could be removed from the network. After finishing the process of the network, these nodes could be checked for its consistency with the desired state at the end. These simple component nodes make a rooted tree which is directed toward the root. A schematic view is shown in Fig 7. We can replace them in all logical formulas by the state of their input nodes in earlier time steps, or in some cases by nodes' initial states. As an example, consider the schematic view of the network shown in Fig 8. In this network, we can replace v 3 in all transition formulas by its inputs and change the network accordingly, as it is shown in Fig 8. The resulting network is a simple component-free extended network.
By this technique, we will have an extended network without any simple SCC. We provide a proof of correctness of our algorithm on these extended networks. While we proved the correctness of the provided algorithm, on these extended networks, as a direct result, the correctness of our algorithm on hard components follows immediately. Also, it is shown that simple components are derived deterministically from the state of nodes within hard components. This shows the correctness of the proposed algorithm.
Considering a simple component-free extended network, the proposed algorithm processes it SCC by SCC. The algorithm first generates all possible states for each SCC, and then enumerates them one by one and goes for the next SCC, recursively. In other words, instead of enumerating all the states for the network node by node, it makes groups of nodes (SCCs) and enumerates states for them. It is clear in this case that, the algorithm enumerates states just like a brute force algorithm. Since the brute force algorithm enumerates and checks all the states, its correctness is obvious from its definition. Note that, the proposed algorithm is not just one brute force, but, we reduced the correctness proof of our algorithm to the correctness of a brute force algorithm.

Drosophila Melanogaster's network
In this section, we present the steps of applying the proposed algorithm on Drosophila melanogaster subnetwork [41], which contains 15 nodes. As it can be seen in Fig 9, three external nodes U1, U2 and U3 are added to this network. Evolving functions regarding the added control sequences are shown in Table 4.
First, we partitioned the network's graph into strongly connected components. As it could be seen in Fig 9, graph has only one multi-node component, which is shown through dotted lines. Then, nodes are ordered based on topological sort and the result is shown as numberings beside each node in Fig 9. We treat graph nodes one by one according to their topological ordering. First, "SLP" node is met. This node is a none-branching single node component; therefore, we compute U for this node till t = τ. Since this node is a constant node (a node without entering edge), it will always keep its initial state. The U concerning this node will be filled to this constant value. Then, for time step τ, we check that if this node has attained the desirable state or not. If not, it would be announced that there exist no desired control sequence, and the algorithm terminates.
Next, the second node in topological sort that is the "en" node will be visited. This node is also a none-branching single node component, thus, it is processed the same as "SLP" node.
Then the "EN" node is considered. It is a branching node, but since it has no control node before itself, it can attain only one state in every time step. Therefore this node is treated like "SLP" and "EN" nodes. For "ci" and "CI" nodes, the very procedure will go. Now we reach a multi-node component. For this strongly connected component, the dynamic programming is computed, then it is checked if there is a possibility to attain the desired state in time step τ for the nodes inside this component or not. If we can attain the desirable state within the time step τ, the U for each node would be filled. Then, branching nodes which have outgoing edges from the component, would be handled as "CIR" node to check if these nodes are possible to attain two states of 0 and 1 in a time step. If they could, a state sequence that we can attain the desirable state, considered and the algorithm would be recalled for the network's remaining nodes.
Then nodes "wg", "WG", "hh" and "HH" are considered. these nodes are none-branching single node components. Note that, if for example "CIR" with two state sequences can attain the desirable state, and for the first state sequence, it would not be possible for the next nodes (e.g. "wg") to attain the desirable state, the algorithm will return and for the sequences of the second state of "CIR", it would check the remaining nodes. In the case that all the nodes were possible to attain the desirable state, it will be announced that there exists a control sequence that causes attaining the desirable state in time step τ.

T-Cell receptor kinetics' network
We apply the proposed algorithm to the Boolean network model of T-cell receptor kinetics [40]. The multi-node component and the order of nodes based on a topological sort are depicted in Fig 10. As it is obvious in this figure, this model of the Boolean network has 40 In T-cell receptor kinetics, the network which can be seen in Fig 10, arrows with pointed heads represent activation and dashed arrows with bar heads represent inhibition (the dashed arrows represent "not" which is related to that node). The network is represented as a network of "or"s of "and"s. Thus, large filled circles representing "and" of their inputs, while when edges are going to a node, state of the node is determined according to the "or" of its incoming edges. For example

Fynðt þ 1Þ ¼ ðCD45ðtÞ^TCRbindðtÞÞ _ ðCD45ðtÞ^LckðtÞÞ
Our suggested algorithm on the Boolean network model of Drosophila melanogaster was implemented for an initial state (in time step t = 0) and a desired state (in time step t = 6), which is depicted in Table 5. It is noteworthy to say that there exists a control sequence to reach the desirable state in time step t = 6. Also, the Datta et al algorithm was implemented in this dataset, for comparison. Note that, since this dataset has 22 edges and 15 internal nodes, the algorithm presented by Akutsu et al requires a tree which has 8 fewer edges from this network (H = 8).

Comparison with previous methods
With initial and desired states that are mentioned in Table 6, we evaluated the proposed algorithm on the Boolean network of Drosophila melanogaster and compared it with previous algorithms. This dataset contains 22 edges and 15 internal nodes, the algorithm presented by Akutsu et al requires a tree which has 8 fewer edges from this network. Thus all combinations of 8 edges are to be removed from the graph that makes it very slow for this case.
The initial state and desired state for T-cell receptor model is shown in Fig 10. This dataset has 40 internal nodes (n = 40) and 3 external nodes (m = 3). We compared the proposed algorithm with the algorithm of Datta et al and Akutsu et al on these initial and desired conditions. Note that there is a control sequence for the above mentioned desirable state in time step t = 5. Akutsu et al [39] first provide a polynomial time algorithm for directed trees. However, their algorithm does not work on general graphs. For general graphs, they remove some edges from the graph to obtain a tree and then search for all possible solutions. Among all the solutions, if one is compatible with the removed edges, their algorithm reports the solution. Note that, since they have to obtain all solutions for the obtained trees, their algorithm on general graphs is not a polynomial time algorithm anymore.
The algorithm presented by Akutsu  In some applications, not only a valid control strategy is needed, but also among valid control strategies (that lead to the desired state), the one with some optimality is preferred. To address this requirement we added the ability to compare valid control strategies for the proposed algorithm. Actually, the proposed algorithm finds a valid control strategy that has minimum switching of control nodes, i.e. from 1 to 0 or 0 to 1. Adding this ability to the proposed algorithm makes it more like an optimization algorithm.
The results of the comparison of our algorithm and previous algorithms implemented on a PC with Dual-Core 2.5GHz CPU, 2G RAM are depicted in Table 7.

Conclusion
In this paper, we proposed an algorithm for network control problem that improves the running time of previously provided algorithms. The extent of improvements and efficiency of our proposed algorithm depends on the size of the multi-node components of the network and also on the positioning of the control nodes or more generally, on the accessibility of both states of 0 and 1 in the same time steps for branching nodes. Since in the proposed algorithm, the nodes are met based on a topological order, the states of the entering nodes of each node are calculated before visiting that node. As a result, in each node, it would be possible to handle states from the initial time step to the desirable time. This would cause early detection of the case that if a node is not possible to attain the desirable state, thus, there be no need to check all the nodes to time step τ, and the lack of a control sequence is reported early. Although the proposed algorithm shows to be efficient on real control networks, however, if all the nodes of the network are in one strongly connected components, the proposed algorithm is not performing well anymore. Handling this case could be future work.   (Table 6) and Drosophila melanogaster's network (Table 5).

Algorithm T-Cell Drosophila
Algorithm of Datta