On control of singleton attractors in multiple Boolean networks: integer programming-based method

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 anti-cancer 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 programming-based 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 multiple-BN 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 post-genome 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 NP-hard 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 average-case or worst-case 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ⁿ × 2ⁿ matrices, which makes it difficult to apply them to large-scale 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 programming-based approach for constructing the attractor control problem for moderate size BNs. The integer linear programming-based (ILP) approach is effective to determine the optimal solutions subject to a series of constraints. Though [15] have shown that attractor control is ${\text{Σ}}_2^p\text{-hard}$, the ILP-based method can be applied to medium-size 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 anti-cancer 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₁, v₂, ⋯, v _n } and a list of Boolean functions F = {f₁, f₂, ⋯, f _n } where f _i : {0, 1}ⁿ → {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

$$v_i\left(t+1\right)=f_i\left(v\left(t\right)\right),i=1,2,\cdots ,n$$

where

$$v\left(t\right)=\left[v_1\left(t\right),\cdots ,v_n\left(t\right)\right]$$

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.

$$\begin{array}{c}{v}_1\left(t+1\right)=v_3\left(t\right),\hfill \\ v_2\left(t+1\right)=\overline{v_1\left(t\right)},\hfill \\ v_3\left(t+1\right)=v_1\left(t\right)\wedge \overline{v_2\left(t\right)}.\hfill \end{array}$$

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ⁿ possible states. Thus, the state transition table and the state transition diagram have 2ⁿ rows and 2ⁿ vertices respectively.

Table 1 The truth table

Table 2 Results on Simultaneous Attractor Control

Figure 1

Example of a Boolean network.

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.,

$$v\left(p\right)=f(f\left(\cdots f\left(v\left(0\right)\right)\cdots \right)=v\left(0\right),$$

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₁(V, F₁) and N₂(V, F₂), representing a cancer cell and a normal cell, respectively. Here, V is a set of genes (nodes) while F₁ and F₂ are sets of regulation functions (i.e., sets of Boolean functions along with input genes) for N₁ and N₂, 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 anti-cancer 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},\cdots ,v_{i_m}$ and these nodes are assigned by Boolean values of $b_{i_1},\cdots ,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₁, ⋯, 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}ⁿ to the set of real numbers. For simplicity, we assume these functions are given in linear form as follows:

$$h\left(v\right)={\displaystyle \sum _i}{\alpha }_i\cdot \left(1-w_i\right)\cdot v_i\text{and}g\left(v\right)={\displaystyle \sum _i}{\beta }_i\cdot \left(1-w_i\right)\cdot v_i$$

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

$$G\left(v\right)=h\left(v\right)+g\left(v\right)$$

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 0-1 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, α₂ = 2, α₃ = 3 and m = 1. If v₁ 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₂, 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₃, 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₁) 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₁) and a singleton attractor for normal cells (N₂). To ensure that the control can cause significant damage to the cancer cells, we introduce a threshold ξ₁ for N₁ to be given later. The score function G(v) becomes g(v). It is possible that there exist multiple singleton attractors for N₂, so we maximize the minimum score of the singleton attractors and give the threshold θ₁ 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 θ₁ and ξ₁,

Problem: Find 0-1 assignment to m control nodes where the minimum score of singleton attractors of N₂ is no less than θ₁, and the score of the singleton attractor for N₁ is greater than ξ₁, 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₂), and whether there still exists a singleton attractor for cancer cells (N₁). We give a threshold ξ₂ 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₁, we consider maximizing the minimum score function integrating the threshold θ₂ for the minimum score of N₁.

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 θ₂ and ξ₂,

Problem: Find m nodes and a 0-1 assignment of these control nodes for which the minimum score of singleton attractors for N₁ is no less than θ₂ and the score of singleton attractors for N₂ is greater than ξ₂, 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 0-1 variable that corresponds to v _i for N₁ and N₂, respectively.

ILP for ATTRACTOR DETECTION

Here we review how ATTRACTOR DETECTION is formalized as ILP in [15]. Define δ _b (x) by

$${\delta }_b\left(x\right)=\left\{\begin{array}{c}x\text{if}b=1\hfill \\ \overline{x}\text{otherwise}.\hfill \end{array}\right.$$

Then if Boolean function has k inputs, we represent it by

$$f_i\left(x_{i_1},\cdots ,x_{i_k}\right)=\underset{b_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{\vee }{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)\wedge {\delta }_{b_1}\left(x_{i_1}\right)\wedge \cdots \wedge {\delta }_{b_k}\left(x_{i_k}\right).$$

Since Boolean constraints cannot be used in ILP directly, we define τ _b (x) as

$${\tau }_b\left(x\right)=\left\{\begin{array}{cc}x\hfill & \text{if}b=1\hfill \\ 1-x\hfill & \text{otherwise}.\hfill \end{array}\right.$$

We introduce $x_{i,b_{i_1}\cdots b_{i_k}}$ in order to represent whether each 0-1 assignment for $\left(b_{i_1},\cdots ,b_{i_k}\right)$ with $f_i\left(b_{i_1},\cdots ,b_{i_k}\right)=1$ is satisfied. If $f_i\left(b_{i_1},\cdots ,b_{i_k}\right)=1$, we give constraints

$$\begin{array}{c}{x}_{i,b_{i_1}\cdots b_{i_k}}\ge \left({\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right)\right)-\left(k-1\right),\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}\le \frac{1}{k}{\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right),\hfill \end{array}$$

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},\cdots ,b_{i_k}\right)=0$, a constraint $x_{i,b_{i_1}\cdots b_{i_k}}=0$ is simply added. These constraints guarantee that

$$x_{i,b_{i_1}\cdots b_{i_k}}=1$$

if and only if

$$f_i\left(b_{i_1},\cdots ,b_{i_k}\right)\wedge {\delta }_{b_1}\left(x_{i_1}\right)\wedge \cdots \wedge {\delta }_{b_k}\left(x_{i_k}\right)=1.$$

Finally, for every x _i , we add the following constraints

$$x_i\le {\displaystyle \sum _{b_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k}}{x}_{i,b_{i_1}\cdots b_{i_k}}$$

and

$$x_i\ge \frac{1}{2^k}{\displaystyle \sum _{b_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k}}{x}_{i,b_{i_1}\cdots b_{i_k}}.$$

The above constraints are included so as to ensure that $x_i=f\left(x_{i_1},\cdots ,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₁ is determined by $x_1=f_1\left(x_2,x_3\right)=x_2\vee {\overline{x}}_3$. Then f₁(x₂, x₃) can be represented as

$$\begin{array}{cc}{f}_1\left(x_2,x_3\right)\hfill & =\left(f_1\left(0,0\right)\wedge {\overline{x}}_2\wedge {\overline{x}}_3\right)\vee \left(f_1\left(0,1\right)\wedge {\overline{x}}_2\wedge x_3\right)\hfill \\ \left(f_1\left(1,0\right)\wedge x_2\wedge {\overline{x}}_3\right)\vee \left(f_1\left(1,1\right)\wedge x_2\wedge x_3\right)\hfill \\ =\left({\overline{x}}_2\wedge {\overline{x}}_3\right)\vee \left(x_2\wedge {\overline{x}}_3\right)\vee \left(x_2\wedge x_3\right).\hfill \end{array}$$

Thus, this Boolean function can be converted into the following inequalities:

$$\left\{\begin{array}{ccc}{x}_{1,00}\hfill & \ge \hfill & \left(1-x_2\right)+\left(1-x_3\right)-1=-x_2-x_3+1,\hfill \\ x_{1,00}\hfill & \le \hfill & \frac{1}{2}\left(1-x_2+1-x_3\right),\hfill \\ x_{1,01}\hfill & =\hfill & 0,\hfill \\ x_{1,10}\hfill & \ge \hfill & x_2+\left(1-x_3\right)-1=x_2-x_3,\hfill \\ x_{1,10}\hfill & \le \hfill & \frac{1}{2}\left(x_2+1-x_3\right),\hfill \\ x_{1,11}\hfill & \ge \hfill & x_2+x_3-1,\hfill \\ x_{1,11}\hfill & \le \hfill & \frac{1}{2}\left(x_2+x_3\right),\hfill \\ x_1\hfill & \le \hfill & x_{1,00}+x_{1,01}+x_{1,10}+x_{1,11},\hfill \\ x_1\hfill & \ge \hfill & \frac{1}{4}\left(x_{1,00}+x_{1,01}+x_{1,10}+x_{1,11}\right).\hfill \end{array}\right.$$

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 ₂

subject to

$$\left\{\begin{array}{c}{x}_{i,b_{i_1}\cdots b_{i_k}}\ge \left({\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right)\right)-\left(k-1\right),\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}\le \frac{1}{k}{\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right),\hfill \\ \text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k,\text{such}\text{that}{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)=1,\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}=0,\text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k,\text{such}\text{that}{f}_i(b_{i_1},\cdots ,b_{i_{\text{k}}}\text{)}\text{ = }\text{0,}\hfill \\ x_i\le {\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{x}_{i,b_{i_1}\cdots b_{i_k}},\text{for}\text{all}i\in \left[1\cdots n\right],\hfill \\ x_i\ge \frac{1}{2^k}{\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{x}_{i,b_{i_1}\cdots b_{i_k}},\hfill \\ x_i\in \left\{0,1\right\},\text{for}\text{all}i\in \left[1\cdots n\right]\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}\in \left\{0,1\right\},\text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k.\hfill \end{array}\right.$$

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₂". 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},x_{i,t,b_{i_1},\cdots ,b_{i_k}}\text{for}t\in \left\{0,1,\cdots ,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₁.

Let x _i be

$$x_i=\left\{\begin{array}{c}{p}_i\text{if}{w}_i=0;\hfill \\ z_i\text{if}{w}_i=1.\hfill \end{array}\right.$$

This function can be represented by

$$\begin{array}{c}{p}_i-w_i\le x_i\le p_i+w_i,\\ z_i-\left(1-w_i\right)\le x_i\le z_i+\left(1-w_i\right).\end{array}$$

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 0-1 assignment. Similarly, for N₂ of the normal cell, another variable s _i is defined as

$$s_i=\left\{\begin{array}{c}{q}_i\text{if}{w}_i=0;\hfill \\ z_i\text{if}{w}_i=1.\hfill \end{array}\right.$$

This function can also be represented by

$$\begin{array}{c}{q}_i-w_i\le s_i\le q_i+w_i\\ z_i-\left(1-w_i\right)\le s_i\le z_i+\left(1-w_i\right).\end{array}$$

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 0-1 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₁ and N₂ take the following forms, respectively:

$$h\left(x\right)={\displaystyle \sum _i}{\alpha }_i\cdot \left(1-w_i\right)\cdot x_i$$

and

$$g\left(s\right)={\displaystyle \sum _i}{\beta }_i\cdot \left(1-w_i\right)\cdot s_i.$$

We can assume without loss of generality that α _i ≥ 0 and β _i ≤ 0. Then, for the first BN, N₁, 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{Σ}}_i{\alpha }_i\cdot 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₂, we introduce the additional binary variable γ _i such that γ _i = -β _i . Thus, the score function for g(s) becomes ${\text{Σ}}_i{\gamma }_i\cdot \left(1-w_i\right)\cdot \left(1-s_i\right)$. Furthermore, the score function can be converted into ${\text{Σ}}_i{\gamma }_i\cdot 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{Σ}}_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\cdot u_i+{\gamma }_i\cdot r_i$

subject to

$$\left\{\begin{array}{c}{x}_{i,b_{i_1}\cdots b_{i_k}}\ge \left({\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right)\right)-\left(k-1\right),\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}\le \frac{1}{k}{\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(x_{i_j}\right),\hfill \\ \text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\text{such}\text{that}{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)=1\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}=0,\hfill \\ \text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\text{such}\text{that}{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)=0\hfill \\ p_i\le {\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{x}_{i,b_{i_1}\cdots b_{i_k}},\text{for}\text{all}i\in \left[1\cdots n\right]\hfill \\ p_i\ge \frac{1}{2^k}{\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{x}_{i,b_{i_1}\cdots b_{i_k}},\hfill \\ p_i-w_i\le x_i\le p_i+w_i,\hfill \\ z_i+w_i-1\le x_i\le z_i-w_i+1,\hfill \\ u_i\le x_i,u_i\le 1-w_i,\text{for}\text{each}i\in \left[1\cdots n\right]\hfill \\ s_{i,b_{i_1}\cdots b_{i_k}}\ge \left({\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(s_{i_j}\right)\right)-\left(k-1\right),\hfill \\ s_{i,b_{i_1}\cdots b_{i_k}}\le \frac{1}{k}\left({\sum }_{j=1,\cdots ,k}{\tau }_{b_{i_j}}\left(s_{i_j}\right)\right),\hfill \\ \text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\text{such}\text{that}{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)=1\hfill \\ s_{i,b_{i_1}\cdots b_{i_k}}=0,\text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\text{such}\text{that}{f}_i\left(b_{i_1},\cdots ,b_{i_k}\right)=0\hfill \\ q_i\le {\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{s}_{i,b_{i_1}\cdots b_{i_k}},\text{for}\text{all}i\in \left[1\cdots n\right]\hfill \\ q_i\ge \frac{1}{2^k}{\sum }_{b_{i_1},\cdots b_{i_k}\in {\left\{0,1\right\}}^k}{s}_{i,b_{i_1}\cdots b_{i_k}},\hfill \\ q_i-w_i\le s_i\le q_i+w_i,\hfill \\ z_i+w_i-1\le s_i\le z_i-w_i+1,\hfill \\ r_i\le 1-s_i,r_i\le 1-w_i,\text{for}\text{each}i\in \left[1\cdots n\right]\hfill \\ x_i,s_i,p_i,q_i,z_i,w_i,u_i,r_i\in \left\{0,1\right\},\text{for}\text{all}i\in \left[1\cdots n\right],\hfill \\ x_{i,b_{i_1}\cdots b_{i_k}}\in \left\{0,1\right\},\text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\hfill \\ s_{i,b_{i_1}\cdots b_{i_k}}\in \left\{0,1\right\},\text{for}\text{all}i\in \left[1\cdots n\right]\text{and}{b}_{i_1}\cdots b_{i_k}\in {\left\{0,1\right\}}^k\hfill \\ {\sum }_{i=1,\cdots ,n}{w}_i=m.\hfill \end{array}\right.$$

We denote this ILP formulation as ILP-A. 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},\cdots ,v_{i_m}\right)$ has been selected as the control nodes with 0-1 value $B^{\prime }=\left(b_{i_1},\cdots ,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₁, ⋯, i _m }. The objective function can be replaced by

$$\text{“Minimize }{\displaystyle \sum _{i\notin I}}{\alpha }_i{x}_i+{\gamma }_i\left(1-s_i\right)”$$

and the constraints replaced u _i and r _i by

$$\left\{\begin{array}{cc}{x}_i=z_i,s_i=z_i,w_i=1,z_i=b_i\hfill & \text{for}\text{all}i\in I\hfill \\ w_i=0\hfill & \text{for}\text{all}i\notin I.\hfill \end{array}\right.$$

The resulting ILP is denoted by ILP-A'.

From the perspective of avoiding examination of the examined (V′, B′), we will modify ILP-A. In other words, we try to find node-value 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 node-value 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₁, then we can also get different node-value pairs for N₂. Thus, we shall consider one of the networks, i.e., N₁. Assume that $x^k=\left(x_1^{\left(k\right)},x_2^{\left(k\right)},\cdots ,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)}\text{if}{w}_i=1$, and otherwise $x_i^{\left(k\right)}=-1$. Then for each k, we add the following linear inequality:

$${\displaystyle \sum _{x_i^{\left(k\right)}\ne -1}}\left(\delta \left(x_i^{\left(k\right)},1\right)\left(-z_i^{\left(k+1\right)}-w_i^{\left(k+1\right)}\right)+\delta \left(x_i^{\left(k\right)},0\right)\left(z_i^{\left(k+1\right)}-w_i^{\left(k+1\right)}\right)\right)\ge 1-{\displaystyle \sum _{x_i^{\left(k\right)}\ne -1}}\left(1+x_i^{\left(k\right)}\right),$$

where δ(x, y) is the delta function,

$$\delta \left(x,y\right)=\left\{\begin{array}{ccc}1\hfill & \text{if}\hfill & x=y\hfill \\ 0\hfill & \text{if}\hfill & x\ne y.\hfill \end{array}\right.$$

This inequality guarantees that the following must hold for at least one of the control nodes:

• if ${\displaystyle 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}.,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 node-value pairs is different from the previous ones. Define ILP-B 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 ILP-B 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 ILP-A'. 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 ILP-A 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

$${\displaystyle \sum _i}{\gamma }_i{u}_i\ge {\xi }_1.$$

(1)

A larger ξ₁ signifies that the cancer cell is damaged substantially. Here we set the ξ₁ as 0.6n. Then, for the normal cell, we define the objective function as ${\sum }_i{\gamma }_i\cdot r_i$. It should be noted that for this problem it is possible that there does not exist any singleton attractor for N₂ since the control set must satisfy the inequality (1) ensuring significant damage to the cancer cells. In terms of ILP-A', we aim at finding the maximum of the minimum score of N₂ for the case of having multiple singleton attractors. Since we consider the multiple singleton attractors for N₂, we do not need to include constraints about N₁, so we can replace r _i with

$$\begin{array}{ccc}{s}_i=z_i,w_i=1,z_i=b_i\hfill & \text{for}\text{all}\hfill & i\in I\hfill \\ w_i=0\hfill & \text{for all}\hfill & i\notin I\hfill \end{array}$$

Consequently, the objective function becomes "Minimize ${\sum }_{i\notin I}\gamma i\left(\text{1}-s_i\right)$". As for ILP-B, to avoid examining the previously obtained (V', B') for N₂, 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:

$${\displaystyle \sum _{s_i^{\left(k\right)}\ne -1}}\left(\delta \left(s_i^{\left(k\right)},1\right)\left(-z_i^{\left(k+1\right)}-w_i^{\left(k+1\right)}\right)+\delta \left(s_i^{\left(k\right)},0\right)\left(z_i^{\left(k+1\right)}-w_i^{\left(k+1\right)}\right)\right)\ge 1-{\displaystyle \sum _{s_i^{\left(k\right)}\ne -1}}\left(1+s_i^{\left(k\right)}\right),$$

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 ILP-B 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 ILP-A'. Output (V', B') and halt, if the minimum score is greater than θ₁.

Here, we set θ₁ to be 0.6n, because if θ₁ ≪ 0.6n, then the desired attractor almost always exists, whereas if θ₁ ≫ 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 ILP-A 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

$${\displaystyle \sum _i}{\gamma }_i{r}_i\ge {\xi }_{\text{2}}.$$

(2)

We set the ξ₂ to be 0.6n, and we define the objective function as ${\sum }_i{\alpha }_i\cdot u_i$. The singleton attractor for N₁ 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 ILP-A', we do not need to add any constraints about N₂ for ILP-A' so we replace u _i by

$$\begin{array}{ccc}{x}_i=z_i,w_i=1,z_i=b_i\hfill & \text{for}\text{all}\hfill & i\in I\hfill \\ w_i=0\hfill & \text{for all}\hfill & i\notin I\hfill \end{array}$$

The objective function is "Minimize ${\sum }_{i\notin I}{\alpha }_i{x}_i.$ Except these modifications, ILP-B 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 ILP-B 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₁ exist.

• 3. For (V', B'), calculate the singleton attractors with the minimum score by ILP-A'. Output ((V', B') and halt, if the minimum score is not less than θ₂.

For this problem, we set θ₂ as 0.6n because it often occurs that the algorithm can find the desired attractor if θ₂ ≪ 0.6n, and that it hardly finds the desired attractor if θ₂ ≫ 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.

Table 3 Results on Attractor Control for the Normal Cells Given that the Control Damage the Cancer Cells Significantly

Table 4 Results on Attractor Control for the Cancer Cells Given that the control has limited Damage to the Normal Cells

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 medium-size networks.

Figure 2

Normal cell vs cancer cell.

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 ILP-based 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 Σp-hard, the experimental results have shown the efficiency of our proposed method. Furthermore, this method was seen to be useful for solving medium-size 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 non-linear costs for the scoring functions is of interest for future work.