 Proceedings
 Open Access
 Published:
On control of singleton attractors in multiple Boolean networks: integer programmingbased method
BMC Systems Biology volume 8, Article number: S7 (2014)
Abstract
Background
Boolean network (BN) is a mathematical model for genetic network and control of genetic networks has become an important issue owing to their potential application in the field of drug discovery and treatment of intractable diseases. Early researches have focused primarily on the analysis of attractor control for a randomly generated BN. However, one may also consider how anticancer drugs act in both normal and cancer cells. Thus, the development of controls for multiple BNs is an important and interesting challenge.
Results
In this article, we formulate three novel problems about attractor control for two BNs (i.e., normal cell and cancer cell). The first is about finding a control that can significantly damage cancer cells but has a limited damage to normal cells. The second is about finding a control for normal cells with a guaranteed damaging effect on cancer cells. Finally, we formulate a definition for finding a control for cancer cells with limited damaging effect on normal cells. We propose integer programmingbased methods for solving these problems in a unified manner, and we conduct computational experiments to illustrate the efficiency and the effectiveness of our method for our multipleBN control problems.
Conclusions
We present three novel control problems for multiple BNs that are realistic control models for gene regulation networks and adopt an integer programming approach to address these problems. Experimental results indicate that our proposed method is useful and effective for moderate size BNs.
Background
In the postgenome era, we observe rapid development in systems biology, a field focusing on interactions among the components of biological systems. Gene regulatory networks (genetic networks, in short) have been proposed to better understand the interaction of various kinds of genes, proteins and molecules. Many formalisms have been developed as models of genetic regulation processes, in particular, Boolean network (BN) has received substantial attention owing to its capacity to capture the switching behavior of genetic processes. Furthermore, its dynamic process is rich and complex and can provide insight into the global behavior of large genetic regulatory networks. A BN is a simple mathematical network: each gene (node) takes either 1 (active) or 0 (inactive), and the state of a gene is regulated by several genes called its input genes via its Boolean functions. It is to be noted that the states of genes can be updated synchronously and asynchronously. Thus synchronous BNs and asynchronous BNs have been proposed to model the two different behaviors. In this paper, synchronous BNs are considered since existing control studies of BNs are based on synchronous BNs, and detection of singleton attractor is equivalent for both models. Though synchronous BNs may be too simple compared with asynchronous BNs, they are effective to model and analyze some properties of real gene networks. Indeed, they have been used to analyze D. melanogaster embryo development, and the robustness of the genetic networks of S. cerevisiae, E. coli, B. subtilis, D. melanogaster and A. thaliana [1].
Many studies have been done on the analysis of attractors for randomly given BNs [2, 3]. The main reason is that different attractors can be regarded as different cell types [2]. In [4–7], several approaches have been developed to efficiently identify and/or enumerate attractors in BNs. It was reported that detection of a singleton attractor (i.e., a fixed state) is an NPhard problem [8]. Devloo et al. studied a method by transformation to a constraint satisfaction problem [4], while Irons proposed another method based on the use of small subnetworks [6]. Garg et al. [5] proposed a method based on Binary Decision Diagrams (BDDs). However, the averagecase or worstcase complexity was not analyzed theoretically in these studies. Zhang et al. [7] developed algorithms to enumerate singleton and small attractors, and also studied the time complexities of average cases for these algorithms. Furthermore, Akutsu et al. [9] developed algorithms with guaranteed worst case time complexity for the singleton attractor detection for BNs with limited Boolean functions. In addition, Datta et al. [10, 11] proposed methods for finding a control sequence for Probabilistic BNs (PBNs), a probabilistic extension of BNs. They defined the control problem with minimization of the sum of the total control cost and the cost of the final state. The cost of control is the cost to apply control inputs in some specified states, with the higher terminal costs corresponding to undesirable states. Their method is also applicable to finding actions of control for BNs, since BNs are special versions of PBNs. But, it is necessary for all these approaches to deal with 2^{n} × 2^{n} matrices, which makes it difficult to apply them to largescale BNs. Consequently, Chen et al. [13] and Kobayashi and Hiraishi [14] proposed integer programming based approaches to variants of the PBN control problem. Akutsu et al. [15] also proposed an integer programmingbased approach for constructing the attractor control problem for moderate size BNs. The integer linear programmingbased (ILP) approach is effective to determine the optimal solutions subject to a series of constraints. Though [15] have shown that attractor control is ${\text{\Sigma}}_{2}^{p}\text{hard}$, the ILPbased method can be applied to mediumsize BNs. All of these works focused their approaches on a single randomly generated BN (i.e., one cell type), leaving the analysis of multiple BNs (various cell types) unaddressed. However, in real situations, there are different kinds of cell types, so it is more realistic to perform attractor analysis for multiple BNs. Thus, the development of attractor control problems for multiple BNs is an important and interesting challenge. Here, we propose an integer linear programming approach for constructing novel attractor controls for multiple BNs, and we provide three novel formulations of the attractor control problem to model various realistic cases. We consider simultaneous attractor control for multiple BNs and, specifically, focus on analyzing attractor control for two BNs that each corresponds to normal cells or cancer cells. Our objective, motivated by the fact that anticancer drugs act in both normal and cancer cells, is to find the same control (i.e., letting some genes be always active (1) and some genes always inactive (0)) for both networks. We aim at finding a control that can damage cancer cells significantly, while causing only limited damage to normal cells.
As another variant of attractor control, we find singleton attractors for normal cells with a guaranteed damaging level for cancer cells. In other words, we try to investigate whether there exists a control set that can damage cancer cells, while at the same time looking for a singleton attractor for normal cells under this control.
It is also of interest to consider finding a singleton attractor for cancer cells with limited damage to normal cells, to determine if there exists a control set ensuring no damage to normal cells, and under such a control set, to search for a singleton attractor for cancer cells.
To utilize ILP, all instances of the original problem are converted into ILP, to be able to apply an existing solver. These problems are transformed into ILP in a similar and systematic way to be shown later. The experimental results, using artificially generated BNs and realistic BNs, have demonstrated both the efficiency and the effectiveness of our proposed method.
Problem formulations
In this section, we first give a brief introduction to BNs. Since the attractor control problem is based on the attractor detection problem (ATTRACTOR DETECTION), we begin with a brief introduction to that problem. We then define three variants of the attractor control problem: simultaneous attractor control for two BNs (cancer cell vs normal cell), attractor control for normal cells under the assumption of significant damaging cancer cells, attractor control for cancer cells under the assumption that normal cells are not damaged.
Boolean networks
A Boolean Network (BN) G(V, F ) is a set of vertices V = {v_{1}, v_{2}, ⋯, v_{ n }} and a list of Boolean functions F = {f_{1}, f_{2}, ⋯, f_{ n }} where f_{ i } : {0, 1}^{n} → {0, 1}. Define v_{ i }(t) to be the state (0 or 1) of vertex v_{ i } at time t. The rules for the regulatory interactions among the genes are then expressed by
where
is called the Gene Activity Profile (GAP). x ∧ y, x ∨ y, x ⊕ y, $\overline{x}$ is used to represent "AND" of x and y, "OR" of x and y, exclusive OR of x and y, and "NOT" of x, respectively. We denote IN (v_{ i }) as the set of relevant input nodes v_{i 1}, ⋯, v_{ ik } to v_{ i }. The number of relevant nodes incoming to v_{ i } is called as indegree of v_{ i }. K is used to represent the maximum indegree of a BN.
The following is an example of a BN.
Each gene v_{ i } is updated by a regulatory function f_{ i }. The corresponding truth table of a BN is given in Table 1. The dynamics of a BN and state transition diagram are shown in Figure 1. For example, the second row of the table means that if the state of BN is [0, 0, 1] at time t, then the state at time t + 1 will be [1, 1, 0]. Similarly, the arc from 001 to 110 signifies that if the state of BN is [0, 0, 1] at time t, then the state will be [1, 1, 0] at time t + 1. It is easily seen that a BN with n nodes has in total 2^{n} possible states. Thus, the state transition table and the state transition diagram have 2^{n} rows and 2^{n} vertices respectively.
Detection of an attractor
In a BN, the GAP v(t + 1) at time t + 1 is determined by the GAP v(t) at time t. When an initial GAP v(0) is given, a BN will finally converge to a set of global states (i.e., a directed cycle in the diagram of state transition). We call this set an attractor. If an attractor has only one state (i.e, v = f(v)), we call it singleton attractor, which is corresponding to a fixed point. Otherwise, if it consists of p distinct states, i.e.,
it is called a cyclic attractor of period p. We can see from Figure 1 that 010 and 101 are two singleton attractors.
Here, the problem of attractor detection is defined to find an attractor for a given BN. It should be noted that finding attractors with large periods is not easy. Thus we study attractor detection in which upper bound of the period is limited by some given threshold pmax. We define this problem as follows. Here we consider the case of p_{ max } = 1 (i.e., the singleton attractor detection problem).
Definition 1: [ATTRACTOR DETECTION]
Instance: a BN and p_{ max }, the maximum period.
Problem: Output an attractor whose period is at most p_{ max }. If such an attractor does not exist, "NULL" should be the output.
For a BN of Table 1, ATTRACTOR DETECTION should output either 010 or 101. Since we only consider the case of p_{ max } = 1 hereafter, we also use v_{ i } to denote the (steady) state of v_{ i }.
Simultaneous attractor control for multiple Boolean networks
A control problem for a single BN was proposed by Akutsu et al. [15]. We consider whether there exists a control for multiple BNs. We consider two BNs, N_{1}(V, F_{1}) and N_{2}(V, F_{2}), representing a cancer cell and a normal cell, respectively. Here, V is a set of genes (nodes) while F_{1} and F_{2} are sets of regulation functions (i.e., sets of Boolean functions along with input genes) for N_{1} and N_{2}, respectively. We try to find the same control (i.e., some genes always active (1) and some genes always inactive (0)) for both networks, considering that anticancer drugs act in both normal and cancer cells. Therefore, we set our endpoint as finding a control that causes limited damage to normal cells but significant damage to cancer cells. Suppose control nodes are chosen as ${v}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{v}_{{i}_{m}}$ and these nodes are assigned by Boolean values of ${b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{m}}$. Then, we call v as a singleton attractor if it satisfies the following conditions (denoted by COND1):

(i)
v_{ i } = b_{ i }; if i ∈ {i_{1}, ⋯, i_{ m }},

(ii)
v_{ i } = f_{ i }(v); otherwise.
Furthermore, in order to evaluate how appropriate each attractor state is, we need a score function h(v) for cancer cells and g(v) for normal cells from {0, 1}^{n} to the set of real numbers. For simplicity, we assume these functions are given in linear form as follows:
where w_{ i } = 1 (resp., w_{ i } = 0) holds if v_{ i } is selected (resp., not selected) as a control node. This means that we do not take the scores of the selected control nodes into account for the score functions (h(v) and g(v)). Here, α_{ i } and β_{ i } are real constants. For example, if α_{ i } = 1 for all i, then v(t) = (1, 1, ⋯, 1) is the state with the highest score (i.e., the desired state for the cancer cells). Similarly, if β_{ i } = 1 for all i, then v(t) = (0, 0, ⋯, 0) is the state with the highest score (i.e., the desired state for the normal cells). We define
as the final score function for these two networks. Since the singleton attractors may not be uniquely determined, considering the worst case, the minimum score of singleton attractors is necessary to be maximized. Because it is difficult to give a direct formulation of this problem, the problem of simultaneous control of attractors is defined with θ, which is a threshold of the minimum score, as follows.
Definition 2 : [SIMULTANEOUS ATTRACTOR CONTROL] (SAC)
Instance: 2 BNs, the number of control nodes m, a score function G, and a threshold θ.
Problem: Find m control nodes and 01 assignment for them where the minimum score of singleton attractors is greater than θ. If no such nodes exist, then "NULL" is the output.
We give an example about attractor control for a given BN, say, the BN described in Table 1. We assume α_{1} = 1, α_{2} = 2, α_{3} = 3 and m = 1. If v_{1} is a control node and 1 is assigned to this node, there exists one singleton attractor 101 with score 3. If 0 is assigned to v_{2}, two singleton attractors 000 and 101 exist. Note that the scores of 000 and 101 are 0 and 4, respectively, so the minimum score is 0. Then, if 0 is assigned to v_{3}, 010 is a singleton attractor with score 2. Though checking the other three cases is necessary, we can conclude that the first case (assigning 1 to v_{1}) gives the solution.
Attractor control for normal cells under the assumption of significant damage to cancer cells
We consider another attractor control variant for multiple networks. We try to investigate whether there exists a control that can guarantee significant damage to cancer cells (N_{1}) and a singleton attractor for normal cells (N_{2}). To ensure that the control can cause significant damage to the cancer cells, we introduce a threshold ξ_{1} for N_{1} to be given later. The score function G(v) becomes g(v). It is possible that there exist multiple singleton attractors for N_{2}, so we maximize the minimum score of the singleton attractors and give the threshold θ_{1} for the minimum score. The problem is formulated as follows:
Definition 3: [Attractor Control under Damaging Cancer cells substantially] (ACDC)
Instance: 2 BNs, a score function G, the number of control node m, and thresholds θ_{1} and ξ_{1},
Problem: Find 01 assignment to m control nodes where the minimum score of singleton attractors of N_{2} is no less than θ_{1}, and the score of the singleton attractor for N_{1} is greater than ξ_{1}, respectively. If such nodes do not exist, then "NULL" is the output.
Attractor control for cancer cells under the assumption of limited damage to normal cells
We also investigate whether there exists a control that ensures limited damage to the normal cells (N_{2}), and whether there still exists a singleton attractor for cancer cells (N_{1}). We give a threshold ξ_{2} for the normal cells that guarantees keeping the normal cells undamaged and convert the score function into G(v)=h(v). In order to obtain the unique singleton attractor for N_{1}, we consider maximizing the minimum score function integrating the threshold θ_{2} for the minimum score of N_{1}.
Definition 4: [Attractor Control under Keeping Normal cells undamaged] (ACKN)
Instance: 2 BNs, a score function G, the number of control node m, and thresholds θ_{2} and ξ_{2},
Problem: Find m nodes and a 01 assignment of these control nodes for which the minimum score of singleton attractors for N_{1} is no less than θ_{2} and the score of singleton attractors for N_{2} is greater than ξ_{2}, respectively. If such nodes do not exist, then "NULL" is the output.
Methods
Integer programming, in particular, integer linear programming (ILP) is to maximize (or minimize) a linear objective function with linear constraints (i.e., linear inequalities and linear equations) with all the variables taking integer values. From here on, either 0 or 1 is assigned to each variable, representing the Boolean values. Furthermore, we introduce x_{ i } and si to represent a 01 variable that corresponds to v_{ i } for N_{1} and N_{2}, respectively.
ILP for ATTRACTOR DETECTION
Here we review how ATTRACTOR DETECTION is formalized as ILP in [15]. Define δ_{ b }(x) by
Then if Boolean function has k inputs, we represent it by
Since Boolean constraints cannot be used in ILP directly, we define τ_{ b }(x) as
We introduce ${x}_{i,{b}_{{i}_{1}}\cdots {b}_{{i}_{k}}}$ in order to represent whether each 01 assignment for $\left({b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{k}}\right)$ with ${f}_{i}\left({b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{k}}\right)=1$ is satisfied. If ${f}_{i}\left({b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{k}}\right)=1$, we give constraints
in which the former constraint forces ${x}_{i,{b}_{{i}_{1}}\cdots {b}_{{i}_{k}}}$ to be 1 if ${\delta}_{{b}_{1}}\left({x}_{{i}_{1}}\right)\wedge \cdots \wedge {\delta}_{{b}_{k}}\left({x}_{{i}_{k}}\right)$ is satisfied, and the latter forces ${x}_{i,{b}_{{i}_{1}}\cdots {b}_{{i}_{k}}}$ to be 0 if it is not satisfied. If ${f}_{i}\left({b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{k}}\right)=0$, a constraint ${x}_{i,{b}_{{i}_{1}}\cdots {b}_{{i}_{k}}}=0$ is simply added. These constraints guarantee that
if and only if
Finally, for every x_{ i }, we add the following constraints
and
The above constraints are included so as to ensure that ${x}_{i}=f\left({x}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{x}_{{i}_{k}}\right)$ holds for each x_{ i }. Thus any feasible solution corresponds to a singleton attractor. A simple example is given to illustrate this method. Assume x_{1} is determined by ${x}_{1}={f}_{1}\left({x}_{2},\phantom{\rule{2.77695pt}{0ex}}{x}_{3}\right)={x}_{2}\vee {\overline{x}}_{3}$. Then f_{1}(x_{2}, x_{3}) can be represented as
Thus, this Boolean function can be converted into the following inequalities:
It is easily seen from this example that the transformation is correct. Integrating all the constraints, we can formulate the integer programming for singleton attractor detection as below.
maximize x _{2}
subject to
We note that ATTRACTOR DETECTION is a decision problem, not an optimization problem, so we do not need an objective function, but we will define an objective function in order to apply ILP. Here, we give the objective function simply as "Maximize x_{2}". We can extend the above method to detect a cyclic attractor whose period is at most p_{ max }, but for that purpose, we will define many more variables ${x}_{i,t},\phantom{\rule{2.77695pt}{0ex}}{x}_{i,t,{b}_{{i}_{1}},\cdots \phantom{\rule{0.3em}{0ex}},{b}_{{i}_{k}}}\text{for}\phantom{\rule{2.77695pt}{0ex}}t\in \left\{0,1,\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},{p}_{max}1\right\}$.
ILP for SIMULTANEOUS ATTRACTOR CONTROL
Consider that each variable vi has two possibilities, i.e,

(i)
v_{ i } is not chosen for a control node (i.e., v_{ i } corresponds to an internal node),

(ii)
v_{ i } is chosen for a control node (i.e., v_{ i } becomes an external node).
In order to specify the state of a variable, we introduce the following variables and constraints for N_{1}.
Let x_{ i } be
This function can be represented by
The above representation means that w_{ i } = 1 corresponds to the case that x_{ i } is selected as a control node, to which z_{ i } gives a 01 assignment. Similarly, for N_{2} of the normal cell, another variable s_{ i } is defined as
This function can also be represented by
In this representation, we can see that w_{ i } = 1 also corresponds to the case that si is selected as a control node, and has the same 01 assignment. This guarantees that these two different BNs are under the same control. For this attractor control problem, we assume that the score functions for N_{1} and N_{2} take the following forms, respectively:
and
We can assume without loss of generality that α_{ i } ≥ 0 and β_{ i } ≤ 0. Then, for the first BN, N_{1}, we try to maximize the score of the internal nodes, so we need to maximize h(x). The score function can be reduced to ${\text{\Sigma}}_{i}\phantom{\rule{2.77695pt}{0ex}}{\alpha}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{u}_{i}$ if we add the constraints u_{ i } ≤ x_{ i } and u_{ i } ≤ 1  w_{ i }, where u_{ i } is a binary variable. In terms of maximizing the score function g(s) for N_{2}, we introduce the additional binary variable γ_{ i } such that γ_{ i } = β_{ i }. Thus, the score function for g(s) becomes ${\text{\Sigma}}_{i}{\gamma}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}\left(1{w}_{i}\right)\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}\left(1{s}_{i}\right)$. Furthermore, the score function can be converted into ${\text{\Sigma}}_{i}\phantom{\rule{0.3em}{0ex}}{\gamma}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{r}_{i}$, if we add the constraints r_{ i } ≤ (1  w_{ i }) and r_{ i } ≤ (1  s_{ i }).
The original aim is that the minimum score of singleton attractors is maximized for these two networks. However, it is difficult to give a direct ILP formalization for this problem, because of the ${\text{\Sigma}}_{2}^{p}\text{hardness}$ of the problem, so as in [15], we firstly formalize an ILP to find a singleton attractor with the maximum score by choosing and controlling m nodes. Incorporating the inequalities mentioned above with the ILP formalization for ATTRACTOR DETECTION, the following ILP formalization for SAC is obtained.
maximize ${\sum}_{i}{\alpha}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{u}_{i}+{\gamma}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{r}_{i}$
subject to
We denote this ILP formulation as ILPA. Since singleton attractors are not always uniquely determined, considering the worst case, it is necessary to maximize the minimum score of the singleton attractors. Suppose that ${V}^{\prime}=\left({v}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{v}_{{i}_{m}}\right)$ has been selected as the control nodes with 01 value ${B}^{\prime}=\left({b}_{{i}_{1}},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{b}_{{i}_{m}}\right)$. These notions are also used in the following two problems, so detection of the attractor with the minimum score with these control nodes can be formulated as follows. Define I = {i_{1}, ⋯, i_{ m }}. The objective function can be replaced by
and the constraints replaced u_{ i } and r_{ i } by
The resulting ILP is denoted by ILPA'.
From the perspective of avoiding examination of the examined (V′, B′), we will modify ILPA. In other words, we try to find nodevalue pairs that are not the same as the previously obtained solutions. This can be tackled by some additional linear inequalities which ensures that the following nodevalue pairs differ from the previous ones. Note that these two networks have the same control, so if we can avoid obtaining the previously examined (V′, B′) for N_{1}, then we can also get different nodevalue pairs for N_{2}. Thus, we shall consider one of the networks, i.e., N_{1}. Assume that ${x}^{k}=\left({x}_{1}^{\left(k\right)},\phantom{\rule{2.77695pt}{0ex}}{x}_{2}^{\left(k\right)},\phantom{\rule{2.77695pt}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}},\phantom{\rule{2.77695pt}{0ex}}{x}_{n}^{\left(k\right)}\right)$ is the k th control previously found, where we let ${x}_{i}^{\left(k\right)}={z}_{i}^{\left(k\right)}\phantom{\rule{2.77695pt}{0ex}}\text{if}\phantom{\rule{2.77695pt}{0ex}}{w}_{i}=1$, and otherwise ${x}_{i}^{\left(k\right)}=1$. Then for each k, we add the following linear inequality:
where δ(x, y) is the delta function,
This inequality guarantees that the following must hold for at least one of the control nodes:

if ${x}_{i}^{(k)}}=1$, either ${z}_{i}^{\left(k+1\right)}=0$ or ${w}_{i}^{\left(k+1\right)}=0$ holds,

otherwise, either ${z}_{i}^{\left(k+1\right)}=1$ or ${w}_{i}^{\left(k+1\right)}=0$ holds.
We consider the case that one of the control nodes ${x}_{i}^{\left(k\right)}$ is $\left(\text{i}.\text{e}.,\phantom{\rule{2.77695pt}{0ex}}{z}_{i}^{\left(k\right)}=1\right)$ for the k th control. If ${w}_{i}^{\left(k+1\right)}=0$ holds for the (k + 1)th control, then in the next iteration, it will not be selected as a control node. Otherwise, ${w}_{i}^{\left(k+1\right)}=1$, so it is still a control node. However, the inequality guarantees that ${z}_{i}^{\left(k+1\right)}$ must be 0, which is different from the previous value $\left({z}_{i}^{\left(k\right)}=1\right)$. Thus, the above inequality ensures that the following obtained set of nodevalue pairs is different from the previous ones. Define ILPB as this modified version in the following context. We can solve the simultaneous attractor control problem for multiple networks through the following algorithm.

1.
Repeat steps 2  3.

2.
Find (V', B') yielding the maximum score of a singleton attractor using ILPB where (V', B') is not the same as any of the already examined nodes/values pairs. "NULL" should be output and halt, if the maximum score is less than θ.

3.
For (V', B'), calculate the singleton attractors with the minimum score by ILPA'. Output (V', B') and halt, if the minimum score is no less than θ.
Note that in the worst case, it may repeat this procedure exponentially many times, but we expect that the procedure will not be repeated so many times, since the expected number of singleton attractors (i.e, per (V', B')) is small, regardless of the total number of nodes (n). How to select the thresholds θ is an important issue in this program. If we know an appropriate threshold in advance, we can simply use such a θ. Here, we let θ be 1.2n, because the desired attractors almost always exist if θ ≪ 1.2n, and it often occurs in our preliminary computational experiments that the algorithm cannot find the desired attractors if θ ≫ 1.2n.
ILP for attractor control under damaging cancer cells substantially
We try to investigate if there exists a control that ensures damaging the cancer cell substantially and, under this condition, identifying singleton attractors for the normal cell. For the ILPA of ACDC (Attractor Control under Damaging Cancer cells substantially), we have to make sure that the control damages the cancer cell significantly by introducing the following inequality
A larger ξ_{1} signifies that the cancer cell is damaged substantially. Here we set the ξ_{1} as 0.6n. Then, for the normal cell, we define the objective function as ${\sum}_{i}{\gamma}_{i}\cdot \phantom{\rule{2.77695pt}{0ex}}{r}_{i}$. It should be noted that for this problem it is possible that there does not exist any singleton attractor for N_{2} since the control set must satisfy the inequality (1) ensuring significant damage to the cancer cells. In terms of ILPA', we aim at finding the maximum of the minimum score of N_{2} for the case of having multiple singleton attractors. Since we consider the multiple singleton attractors for N_{2}, we do not need to include constraints about N_{1}, so we can replace r_{ i } with
Consequently, the objective function becomes "Minimize ${\sum}_{i\notin I}\gamma i\left(\text{1}{s}_{i}\right)$". As for ILPB, to avoid examining the previously obtained (V', B') for N_{2}, we assume that ${s}^{k}=\left({s}_{1}^{\left(k\right)},{s}_{2}^{\left(k\right)},\dots ,{s}_{n}^{\left(k\right)}\right)$ is the k th control previously found, where we let ${s}_{i}^{\left(k\right)}={s}_{i}^{\left(k\right)}$ if w_{ i } = 1, otherwise ${s}_{i}^{\left(k\right)}=1$. Then for each k, we add the
following linear inequality:
We can solve this problem using the following algorithm.

1.
Repeat steps 2  3.

2.
Find (V', B') yielding the maximum score of a singleton attractor using ILPB where (V', B') is not the same as any of the already examined nodes/values pairs and that the control damages the cancer cell significantly. "NULL" should be output and halt, if the singleton attractor does not exist.

3.
For (V', B'), calculate the singleton attractors with the minimum score by ILPA'. Output (V', B') and halt, if the minimum score is greater than θ_{1}.
Here, we set θ_{1} to be 0.6n, because if θ_{1} ≪ 0.6n, then the desired attractor almost always exists, whereas if θ_{1} ≫ 0.6n, then the desired attractor seldom exists.
ILP for attractor control under keeping normal cells undamaged
We also try to investigate whether there exists a control that ensures not damaging the normal cell substantially and under this condition, we find a singleton attractor for the cancer cell. In terms of ILPA for ACKN (Attractor Control under Keeping Normal cells undamaged), considering this control should not damage the normal cell substantially, so we introduce an additional inequality
We set the ξ_{2} to be 0.6n, and we define the objective function as ${\sum}_{i}{\alpha}_{i}\phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{u}_{i}$. The singleton attractor for N_{1} is difficult to find for this problem, because the control set must satisfy the inequality (2) (not damaging the normal cell substantially). Considering that there may exist multiple singleton attractors for the cancer cell, and the worst case, it is necessary to maximize the minimum score of singleton attractors.
Thus for ILPA', we do not need to add any constraints about N_{2} for ILPA' so we replace u_{ i } by
The objective function is "Minimize ${\sum}_{i\notin I}{\alpha}_{i}{x}_{i}.$ Except these modifications, ILPB is exactly the same as in the case of SAC. We can solve this problem through the following algorithm.

1.
Repeat steps 2  3.

2.
Find (V', B') yielding the maximum score of a singleton attractor using ILPB where (V', B') is not the same as any of the already examined nodes/values pairs and that the control causes only limited damage to the normal cell. "NULL" should be output and halt, if no singleton attractor for N_{1} exist.

3.
For (V', B'), calculate the singleton attractors with the minimum score by ILPA'. Output ((V', B') and halt, if the minimum score is not less than θ_{2}.
For this problem, we set θ_{2} as 0.6n because it often occurs that the algorithm can find the desired attractor if θ_{2} ≪ 0.6n, and that it hardly finds the desired attractor if θ_{2} ≫ 0.6n.
Results
Results on simultaneous attractor control
The proposed method for SIMULTANEOUS ATTRACTOR CONTROL was applied to randomly generated BNs. The following data is according to 10 randomly generated pairs of BNs with K = 2, where for each pair of BNs, the BNs have the same nodes but different Boolean functions.

1.
time: average CPU time (seconds) per pair of BNs,

2.
#pos#rep: the number of pairs of BNs where the desired attractors can be found, and the average number of repeats for each such pair of BNs (recall that it is necessary for SAC (Simultaneous Attractor Control) to solve ILP instances repeatedly),

3.
#neg#rep: the number of pairs of BNs for which the desired attractor does not exist, and the average number of repeats for each such pair of BNs.
We set α_{ i } = 1 and γ_{ i } = 1 for all i ∈ {1, 2, ⋯, n}, and we set the maximum number of repeats to be 20. It is possible that this procedure may not decide whether or not the desired attractors exist within 20 repeats in some cases. The table shows that the number of repeats is small even though the number of nodes is large. The reason the CPU time for (140, 14) is greater than that for (160, 16) is that the latter required fewer repeats.
Results on attractor control for normal cells under damaging cancer cell substantially
For this problem, we first guarantee that this control damages the cancer cell substantially and add the additional constraint (1). It is possible that no singleton attractor exists for the normal cell under the condition that the control set guarantees significant damage to the cancer cells. As shown in Table 3, our proposed method is useful also for this problem.
Results on attractor control for cancer cells under keeping normal cells undamaged
The results also suggest that our proposed approach is efficient and effective for solving the problem of attractor control of the cancer cell with no damage to the normal cell guaranteed.
In order to verify our proposed method further, we applied ILP to some real networks (see Figure 5B) in [16]. There are two types of mice in that figure. The left figure is for BALB/c, which is more sensitive to radiation and tends to become cancer, while the right one is for C57BL/6, which is more resistant to radiation and tends not to become cancer. However, the nodes for these two networks are different and the Boolean functions are not given. Thus, we consider utilizing the right figure, which corresponds to the normal cell, and further simulating a cancer cell based on the normal cell (right figure) to guarantee that both networks have the same nodes (see Figure 2(b)). Furthermore, the red and green nodes for the normal cell represent 1 and 0, respectively. We assign most of the OR nodes to the cancer cell, and we see that the majority of nodes in the singleton attractor of the cancer cell are 1. Similarly, we assign most of the AND nodes to the normal cell, and we see that the majority of nodes in the singleton attractor of the normal cell are 0. This may mean that this is a subnetwork activated by cancer. Thus we have obtained the singleton attractor for the cancer cell and normal cell. Since our objective is to transform the state of the cancer cell into that of the normal cell, we try to find a control such that most of the 1 nodes in the cancer cell are converted into 0 nodes. Then we define a score function for both BNs. Specifically, we assume γ_{ i } = 1 and α_{ i } = 1 for all i, and we set n = 36, m = 4 and θ = 1.2·n. We can see that most of the states of the singleton attractor for cancer cells are changed into 0, indicating they are the desired attractors for the cancer cell under this control. Moreover, the CPU time is 0.03 (sec), which suggests that the proposed method might work efficiently for real mediumsize networks.
Discussions
In this paper, we formulated three novel problems, simultaneous attractor control for multiple networks, attractor control for normal cells with guaranteed damage to cancer cells, and attractor control for cancer cells guaranteed not to damage normal cells, and we applied an ILPbased approach for solving these problems in a unified manner. We further investigated the attractor control problem for multiple BNs and validated our proposed method for realistic networks. Though attractor control problems are Σphard, the experimental results have shown the efficiency of our proposed method. Furthermore, this method was seen to be useful for solving mediumsize instances of these problems. The method we proposed might not be the fastest, but it is easy and simple to implement and, furthermore, it has rooms for modifications and extensions. In particular, the use of nonlinear costs for the scoring functions is of interest for future work.
References
 1.
Belleza E, Chaos A, Kauffman S, Shmulevich I, Aldana M: Critical Dynamics in Genetic Regulatory Networks: Examples from Four Kingdoms. PLoS One. 2008, 3: e245610.1371/journal.pone.0002456.
 2.
Kauffman SA: The Origins of Order: Selforganization and Selection in Evolution. 1993, New York: Oxford Univ Press
 3.
Samuelsson B, Troein C: Superpolynomial Growth in The Number of Attractors in Kauffman Networks. Phys Rev Lett. 2003, 90 (9): 098701
 4.
Devloo V, Hansen P, Labbé : Identification of Steady States in Large Networks by Logical Analysis. Bulletin of Mathematical Biology. 2003, 65 (6): 10251051. 10.1016/S00928240(03)000612.
 5.
Garg A, Xenarios I, Mendoza L, DeMicheli G: An Efficient Method for Dynamic Analysis of Gene Regulatory Networks and in silico Gene Perturbation Experiments. proceedings of 11th Annual International Conference on Research in Computational Molecular Biology. 2007, Oakland, Galif USA, 4453: 6276.
 6.
Irons DJ: Improving the Efficiency of Attractor Cycle Identification in BNs. Physis D. 2006, 217 (1): 721.
 7.
Zhang SQ, Hayashida M, Akutsu T, Ching WK, Ng MK: Algorithms for Finding Small Attractors in Boolean Networks. EURASIP Journal on Bioinformatics and Systems Biology. 2007, 2007: 20180
 8.
Akutsu T, Kuhara S, Maruyama O, Miyano S: A System for Identifying Genetic Networks from Genetic Expression Patterns Produced by Gene Disuptions and Overexpressions. Genome Informatics. 1998, 9: 151160.
 9.
Akutsu T, Melkman AA, Tamura T, Yamamoto M: Determining a Singleton Attractor of a Boolean Network with Nested Canalyzing Functions. J Computational Biology. 2011, 18: 12751290. 10.1089/cmb.2010.0281.
 10.
Datta A, Choudhary A, Bittner ML, Dougherty ER: External Control in Markovian Genetic Regulatory Networks. Machine Learning. 2003, 52: 169191. 10.1023/A:1023909812213.
 11.
Datta A, Choudhary A, Bittner ML, Dougherty ER: External Control in Markovian Genetic Regulatory Networks: The Imperfect Information Case. Bioinformatics. 2004, 20 (6): 924930. 10.1093/bioinformatics/bth008.
 12.
Hayashida M, Tamura T, Akutsu T, Zhang SQ, Ching WK: Algorithms and Complexity analyses for Control of Singleton Attractors in Boolean Networks. EURASIP Journal on Bioinformatics and Systems Biology. 2008, 2008: 521407
 13.
Chen X, Akutsu T, Tamura T, Ching WK: Finding Optimal Control Policy in Probabilistic Boolean Networks with Hard Constraints by Using Integer Programming and Dynamic Programming. International Journal of Data Mining and Bioinformatics. 2013, 7: 322343.
 14.
Kobayashi K, Hiraishi K: An Integer Programming Approach to Optimal Control Problems in ContextSensitive Probabilistic Boolean Networks. Automatica. 2011, 47: 12601264. 10.1016/j.automatica.2011.01.035.
 15.
Akutsu T, Zhao Y, Hayashida M, Tamura T: Integer ProgrammingBased Approach to Attractor Detection and Control of Boolean Networks. IEICE TRANS INF & SYST. 2012, E95D2960.
 16.
Snijders AM, Marchetti F, Bhatnagar S, Duru N, et al.: Genetic Differences in Transcript Responses to LowDose Ionizing Radiation Identify Tissue Functions Associated with Breast Cancer Susceptibility. Plos ONE. 2012, 7 (10): e4539410.1371/journal.pone.0045394.
Declarations
Publication of this article was funded by University Grants of Kyoto University, Japan.
This article has been published as part of BMC Systems Biology Volume 8 Supplement 1, 2014: Selected articles from the Twelfth Asia Pacific Bioinformatics Conference (APBC 2014): Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/8/S1.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
The basic idea of this research was proposed by TA after some discussion with WC. YQ and TT analyzed the theoretical details. All computer experiments were conducted by YQ. YQ wrote most parts of the manuscript, and TT edited it. TA and WC supervised the work. All authors read and approved the final manuscript.
Rights and permissions
About this article
Cite this article
Qiu, Y., Tamura, T., Ching, W. et al. On control of singleton attractors in multiple Boolean networks: integer programmingbased method. BMC Syst Biol 8, S7 (2014). https://doi.org/10.1186/175205098S1S7
Published:
Keywords
 Normal Cell
 Boolean Function
 Score Function
 Boolean Network
 Minimum Score