Optimal Control of Gene Regulatory Networks with Effectiveness of Multiple Drugs: A Boolean Network Approach

Developing control theory of gene regulatory networks is one of the significant topics in the field of systems biology, and it is expected to apply the obtained results to gene therapy technologies in the future. In this paper, a control method using a Boolean network (BN) is studied. A BN is widely used as a model of gene regulatory networks, and gene expression is expressed by a binary value (0 or 1). In the control problem, we assume that the concentration level of a part of genes is arbitrarily determined as the control input. However, there are cases that no gene satisfying this assumption exists, and it is important to consider structural control via external stimuli. Furthermore, these controls are realized by multiple drugs, and it is also important to consider multiple effects such as duration of effect and side effects. In this paper, we propose a BN model with two types of the control inputs and an optimal control method with duration of drug effectiveness. First, a BN model and duration of drug effectiveness are discussed. Next, the optimal control problem is formulated and is reduced to an integer linear programming problem. Finally, numerical simulations are shown.


Introduction
In the field of systems biology, there have been a lot of studies on modeling, analysis, and control of gene regulatory networks. Especially, control of gene regulatory networks corresponds to therapeutic interventions, which are realized by radiation, chemotherapy, and so on. In order to develop gene therapy technologies (see, e.g., [1]) in the future, developing control theory of gene regulatory networks is important. Furthermore, in recent years, the important result on control of gene regulatory networks has been obtained in [2]. That is, feedback control of synthetic biological circuits has been implemented, and the experimental result in which cellular behavior is regulated by control has been obtained. This result suggests that control methods of gene regulatory networks can be realized. Motivated by the above background, we study control methods of gene regulatory networks.
Gene regulatory networks are in general expressed by ordinary/partial differential equations with high nonlinearity and high dimensionality. In order to deal with such a system, it is important to consider a simple model, and various models such as Bayesian networks, Boolean networks (BNs) [3], hybrid systems (piecewise affine models), and Petri nets have been developed so far (see, e.g., [4] for further details). In control problems, BNs and hybrid systems are frequently used [5][6][7][8]. In the hybrid systems-based approach, the class of gene regulatory networks are limited to low-dimensional systems, because the computation time to solve the control problem is too long. In a BN, dynamics such as interactions between genes are expressed by the Boolean functions [3]; that is, gene expression is expressed by a binary value (0 or 1). There is a criticism that a Boolean network is too simple as a model of gene regulatory networks (see, e.g., [9]), but this model can be relatively applied to large-scale systems. In addition, since the behavior of gene regulatory networks is stochastic by the effects of noise, it is appropriate that a Boolean function is randomly decided at each time among the candidates of the Boolean functions. From this viewpoint, a probabilistic Boolean network (PBN) has been proposed in [10]. Furthermore, a context-sensitive PBN (CS-PBN) in which the deciding time is randomly selected has been proposed as a general form of PBNs [11,12].
Furthermore, in the control theory of gene regulatory networks, the control input is given by the concentration level of a part of genes; that is, we assume that the concentration level of a part of genes can be arbitrarily determined. However, in the case where this assumption is not satisfied, it is important to consider structural control via external stimuli [13,14]. These controls are realized by multiple drugs, and it is also important to consider multiple effects such as duration of effect and side effects [15]. To our knowledge, a unified method considering these properties has not been proposed so far.
Thus, in this paper, we propose a BN model with two types of the control inputs and an optimal control method with duration of drug effectiveness. The first control input is the control input satisfying the assumption that the binary value is arbitrarily determined. The second control input is called a structural control input herein, and the dynamics, that is, the Boolean functions, are selected among the candidates of the dynamics. However, it is difficult to uniquely select one Boolean function. Hence, we suppose that one Boolean function is selected probabilistically, and the probability distribution is switched by using the structural control input. A structural control method has been discussed in [13,14], but the notion of the structural control input defined in this paper is different from that in those existing methods. Since the proposed BN model has a switch of the probability distribution, it may be regarded as a generalized version of PBNs.
In optimal control of PBNs and CS-PBNs, many results have been obtained so far (see, e.g., [11,12,[16][17][18][19][20][21]). In many existing results, state transition diagrams with 2 nodes (i.e., 2 × 2 transition probability matrices) must be computed for a PBN with states. As a result, in order to compute state transition diagrams, several issues such as memory consumption must be considered in implementation, and it is desirable to directly use a given Boolean function. The authors have proposed in [22] a control method in which state transition diagrams are not computed. In many existing results, we consider finding a control input such that the expected value of the cost function is minimized. In [22], we consider finding a control input such that the lower bound of the cost function is minimized under a certain constraint condition. Owing to this difference, state transition diagrams are not computed in the method in [22], and the optimal control problem is reduced to an integer linear programming (ILP) problem. Also in [16], ILP-based methods were proposed for other optimal control problems, and in those methods, solving multiple ILP problems is required.
In this paper, based on our previously proposed method [22], the optimal control problem with duration of drug effectiveness is reduced to an ILP problem. Since a given Boolean function is directly used, duration of drug effectiveness can be easily described as a linear inequality constraint. The proposed method provides us with a basic in control theory of gene regulatory networks. The conference paper [23] is a preliminary version of this paper. In this paper, we provide improved formulations and explanations, a discussion on duration of drug effectiveness, and a numerical simulation using the large-scale BN.
This paper is organized as follows. In Section 2.1, the Boolean networks with two kinds of the control inputs are proposed. In Section 2.2, duration of drug effectiveness is introduced. In Section 2.3, the optimal control problem is formulated. In Section 2.4, its solution method is proposed. In Section 3, two numerical examples are presented. In Section 4, we conclude this paper.
Notation. Let R denote the set of real numbers. Let {0, 1} denote the set of -dimensional vectors, which consists of elements 0 and 1. Let and 0 × denote the × identity matrix and the × zero matrix, respectively. For simplicity, we sometimes use the symbol 0 instead of 0 × and and the symbol instead of . For a matrix , ln denotes the matrix such that the ( , )th element is given as the natural logarithm of the ( , )th element in . For a matrix , denotes the transpose of .

The Boolean Networks with Control Inputs. A Boolean network (BN) with states is given by
where ∈ {0,1} is the state (e.g., the concentration of genes) and = 0, 1, 2, . . . is the discrete time. The function : {0, 1} → {0, 1} is a given Boolean function with logical operators such as AND (∧), OR (∨), and NOT (¬). If the BN (1) is deterministic, then the next state ( + 1) is uniquely determined for a given ( ). See also Example 1 for an example.
Next, the control inputs are added to a BN (1). For the BN (1) with state, consider two types of the control inputs. First, in a similar way to that of the conventional control method, the control input is added to the BN (1) as follows: where ∈ {0, 1} is the control input; that is, the value of (e.g., the concentration of genes) can be arbitrarily given, and : {0, 1} × {0, 1} → {0, 1} is a given Boolean function. The th element of the state and the th element of the control input are denoted by and , respectively. In the BN (2), ( + 1) is uniquely determined for the given ( ) and ( ). Then, consider the structural control input. Suppose that the candidates of are given by , = 1, 2, . . . , . It will be difficult to select one Boolean function uniquely. In this paper, we assume that one discrete probability distribution is selected among discrete probability distributions. Probabilistic distributions are derived from experimental results, but details are one of the future works. Then, a method for inferring a probabilistic Boolean network will be useful (see, e.g., [24]). Let , denote the probability that the Boolean function is selected in the th discrete probability distribution. Then, hold. In addition, -dimensional binary variables ∈ {0, 1} are assigned to discrete probability distributions, and let denote the th element of . The structural control input corresponds to kinds of external stimuli. Then, the equality constraint is imposed. Here, we show a simple example.
Example 1. As a simple example, consider the simplified model of an apoptosis network in Figure 1 [25]. Then, the Boolean network model expressing this apoptosis network is given by where the concentration level (high or low) of the inhibitor of apoptosis proteins (IAPs) is denoted by 1 , the concentration level of the active caspase 3 (C3a) is denoted by 2 , and the concentration level of the active caspase 8 (C8a) is denoted by 3 . The concentration level of the tumor necrosis factor (TNF, a stimulus) is denoted by and is regarded as the control input. Since the caspase C3a is responsible for cleaving or breaking many other proteins, a high level of the C3a concentration, that is, 2 = 1, implies cell near-death, otherwise, cell survival. As seen in (5), if the concentration of IAP is high ( 1 = 1) or the concentration of the caspase C8a is low ( 3 = 0), then the concentration of C3a becomes low; that is, 2 = 0. On the other hand, 1 and 3 at the next time depend on the value of 2 as well as . In this way, some dynamical interactions exist. See [25,26] for further details. Suppose that = 2 and = 2. Then, as an example of the candidates of the Boolean functions, we consider the following: We suppose that the Boolean function 1 expresses the situation that the dynamics of an apoptosis network are selected with high probability and that the Boolean function 2 expresses the situation that the state is not changed with high probability. By using 1 and 2 , one of the two discrete probability distributions { 1,1 , 1,2 } and { 2,1 , 2,2 } is selected at each time.
A BN with two types of the control inputs includes the probabilistic behavior, and we assume that the probability distribution can be controlled. From these facts, a BN studied in this paper can be regarded as a generalized form of a probabilistic Boolean network (PBN). To explain the relation between the proposed BN model and a PBN, we show a simple example.
Example 2. As a simple example, consider the PBN with three states and one control input. Suppose that the Boolean functions are given as follows: This PBN corresponds to the cases of = 4 and = 1. The candidates of the Boolean functions , = 1, 2, 3, 4, and the probabilities 1, , = 1, 2, 3, 4, are obtained as follows: , 1,4 = 0.06.
In this example, the cardinality of the finite state set {0, 1} 3 is given by 2 3 = 8, and we obtain the state transition diagram of   Figure 2. In the existing solution methods for optimal control of PBNs, the optimal control input is computed using dynamic programming with state transition diagrams.
As shown in this example, computing state transition diagrams with 2 nodes ( is the number of the state) is required in the existing solution methods for optimal control of PBNs with states, and this computation is hard for large-scale systems (see also Section 3.2). Thus, it is important to consider a new solution method. In this paper, for BNs with two types of the control inputs, a solution method using integer programming is proposed based on our previously proposed work in [22]. In the proposed method, computation of state transition diagrams such as that in Figure 2 is not needed.
Remark 3. By adding the candidates of the Boolean functions, BNs with two types of the control inputs can be transformed into BNs with only the structural control input. That is, the control input can be eliminated from (2) by fixing the value of in (2). Then, the number of the candidates of Boolean functions is 2 , and 2 combinations for must be computed in advance. To avoid this computation, we consider two types of the control inputs.

Duration of Drug Effectiveness.
The control input and the structural control input are realized by using multiple drugs. Then, we must consider the multiple effects such as duration of effect and the side effects. In this paper, we focus on the duration of drug effectiveness. In, for example, chemotherapy, therapeutic intervention is generally applied to the target cell in a cyclic manner [15]. Each therapeutic window is started by delivering the drug. The drug delivered is effective on the target cell for some period of time. This is followed by a recovery phase. However, when the drug is not delivered, the drug may be delivered in the timing that is faster than the next time in a cyclic [15]. Therefore, it is necessary to model several situations on duration of drug effectiveness. To model the duration of effect, three parameters , 1 , and 0 are defined for each input (or ). The parameters and 1 have been already defined in [15]. The parameter is the length of the drug effectiveness period. That is, if ( ) = 1, then ( is a decision variable. By using , 1 , 0 , we can consider several situations, and we show two typical examples. Example 4. First, suppose that, for the control input ∈ {0, 1} 1 , , 1 , and 0 are given as = 1, 1 = 3, and 0 = 2, respectively. Consider the case of ( ) = 1. Then, ( + 1) and ( + 2) are uniquely determined as ( + 1) = 1 and ( + 2) = 0, respectively, and ( + 3) is a decision variable. Then, ( + 2) = 0 is the recovery phase, and 1 − − 1 = 1 is its length. In the case of ( ) = 0, the relation ( + 1) = 0 holds, and ( + 2) is a decision variable.
By using the three parameters , 1 , and 0 , several situations on the duration of effect can be modeled (see also [15]). In addition, since these parameters can be given for each (or ), effectiveness of multiple drugs can be evaluated. Thus, in this paper, we consider not only two types of the control inputs but also duration of drug effectiveness.

Optimal Control Problem.
First, the following two notations are defined. Let ( ) denote the probability that some Boolean function is selected at time . In addition, the probability that some time sequence of theBoolean functions For simplicity of notation, ( 1 ), ( 1 +1), . . . , ( 2 ) are omitted in ( 1 , 2 ). Next, for the Boolean networks with states and two types of the control inputs, consider the following optimal control problem. Problem 1. Suppose that, for the Boolean network with states and two types of the control inputs, the initial state (0) = 0 , satisfying 0 ≤ ≤ 1, the control time , the parameters on duration of drug effectiveness ( ) , 1 ( ) , and 0 ( ) are given. Then, for all combinations of the Boolean functions satisfying the constraint find two control input sequences (0), (1), . . . , ( − 1) and (0), (1), . . . , ( − 1) minimizing the lower bound of the cost function subject to the constraint on duration of drug effectiveness, where , ∈ R 1× , ∈ R 1× , and ∈ R 1× are weighting vectors whose elements are nonnegative real numbers.
For simplicity of discussion, a linear function with respect to , , and is considered as a cost function. We consider that a linear cost function is appropriate from the following two reasons.
(i) For a binary variable ∈ {0, 1}, the relation 2 = holds. That is, in the cost function, the quadratic term such as 2 ( ) is not necessary.
(ii) In control of gene regulatory networks, the expression of a certain gene is frequently focused (see, e.g., [18]). For example, in the gene regulatory network related to melanoma, it important to inhibit the concentration level of the gene WNT5A [27]. In this case, it is enough to consider the cost function (12).
Furthermore, in many existing methods on optimal control of PBNs, the expected value of a nonnegative function is frequently used as a cost function (see, e.g., [11,12,[17][18][19][20][21]). However, the expected value of the state must be computed from all combinations of the Boolean functions, and this computation is hard for large-scale systems. To avoid this computation, in this paper, we evaluate the control performance by using the lower bound. If the constraint (11) is not included in Problem 1, then the behaviors are regarded as uncertain (nondeterministic) behaviors, and the best performance is derived in Problem 1. Since the combinations of the Boolean functions selected with low probability are included, performance evaluation is not appropriate. In order to exclude such combinations, we impose the constraint (11). Similar problem formulations have been considered in optimal control of stochastic hybrid systems (see, e.g., [28][29][30]). Thus, since the performance index in this paper is different from that in existing methods, it is difficult to directly compare the performance of the proposed method with those of existing methods. On the other hand, in [22], we discussed this topic from the qualitative viewpoint. In [22], the upper bound is also computed by using the control input such that the lower bound is minimized. If the lower bound and the upper bound are not improved by control, then the expected value will not be improved. Then, it is important to suitably set in the constraint (11). See [22] for further details.
We show an example for setting weighting vectors from the biological viewpoint.
Example 5. Consider the Boolean network expressing an apoptosis network in Example 1 again. For this system, we consider finding a control strategy such that a stimulus is not applied as much as possible, and cell survival is achieved; = 0 implies that a stimulus is not applied to the system, and 1 = 1 and 2 = 0 express cell survival [25]. Then, as one of the appropriate cost functions, we can consider the following cost function: By the coordinate transformation of 1 into 1 − 1 , this cost function can be rewritten as the form of (12).

Solution Method.
We propose a solution method for Problem 1. First, two lemmas are introduced as preparations. Next, Problem 1 is reduced to an integer linear programming (ILP) problem.
As preparations, two lemmas are introduced. To reduce Problem 1 to an ILP problem, it is necessary to transform a Boolean function into a polynomial on the real number field. First, the following lemma [31] is used. Lemma 6. Consider the two binary variables 1 , and 2 . Then, the following relations hold: For example, 1 ∨ ¬ 2 is equivalently transformed into 1 + (1 − 2 ) − 1 (1 − 2 ) = 1− 2 + 1 2 . Furthermore, the product of binary variables such as 1 2 can be linearized by using the following lemma [32].

Lemma 7. Suppose that the binary variables
∈ {0, 1} and ∈ J are given, where J is some index set. Then, = ∏ ∈J is equivalent to the following linear inequalities: where |J| is the cardinality of J.
Here, we define the following vector: Then, by using the natural logarithm, (0, − 1) in (11) is expressed as In this expression, one probability distribution is selected by using ( ), and the probability that a certain Boolean function is selected is determined by ( ). Then, Problem 1 is equivalent to the following problem.
By using Lemma 7, the system (15) and (∑ =1 ln ( )) ( ) can be equivalently expressed in the following linear form: where ( ) := ( ) ( ), and (22) is the linear inequality obtained by applying Lemma 7 to (15). The vector ( ) ∈ {0, 1} is an auxiliary binary variable obtained by using Lemma 7, and the dimension of ( ), that is, , is determined depending on the form of the given Boolean functions. In addition, ( ) is defined as Here in after, for simplicity of notation, ( ) is rewritten as ( ), := [ 0], and ∑ =1 ln ( ) is rewritten as ( ), := [0 ln 1 ln 2 ⋅ ⋅ ⋅ ln ]. Now, we consider transforming Problem A by using (20), (21), and (22). By using obtained from the state equation in (20), we can obtain Next, the inequality constraint where = [ ⋅ ⋅ ⋅ ]. Furthermore, from (22), we can obtain Next, consider the constraint on duration of drug effectiveness. This constraint can be expressed as a Boolean function. Then, by using Lemmas 6 and 7, it can be transformed into the following form: where V and V are binary decision variables with certain dimensions. Deriving a general form of coefficient matrices will be difficult, but for the given , 1 , and 0 , deriving coefficient matrices is easy.
We show two examples.
Example 8. Consider Example 4 again. First, consider the case of = 1, 1 = 3, and 0 = 2. Then, noting explanations in Example 4, we can obtain . . . (32) and in this case, these are equivalent to We explain (2) = V 2 , and V 2 ≤ 1−V 0 as an example. If V 0 = 1, then V 2 ≤ 0 holds. Since V 2 is binary, we can obtain V 2 = 0; that is, (2) = 0. If V 0 = 0, then V 2 ≤ 1 holds, and we can obtain (2) = V 2 . That is, (2) can take on either 0 or 1. From the previous discussion, we see that a pair of (2) = V 2 and Thus, we can obtain the forms of (30) and (31). In the case of = 5 ( is the control time in Problem 1), (30) can be obtained as where 2 = 0. We remark that, in a general case, the product such as = V 0 V 2 must be transformed into linear inequalities by using Lemma 7.
Next, consider the case of = 0, 1 = 3, and 0 = 3. Then, we can obtain the following:  and this case is one of the simplest cases. In the case of = 5 ( is the control time in Problem 1), (30) can be obtained as where 2 = 0, 1 = 0, and 2 = 0.
Problem B can be solved by a suitable solver such as IBM ILOG CPLEX [33].

Results and Discussion
In this section, we show numerical simulations. First, we consider the WNT5A network [14]. Next, in order to evaluate the proposed method from the viewpoint of the computation time, we consider an artificial example.
3.1. WNT5A Network. The gene regulatory network with the gene WNT5A is related to melanoma, and it has been extensively studied (see, e.g., [27]). The BN model ( + 1) = ( ( )) of the WNT5A network is given by the following: 1 ( + 1) = ¬ 6 ( ) , 2 ( + 1) = (¬ 2 ( ) ∧ 4 ( ) ∧ 6 ( )) ∨ {¬ 2 ( ) ∧ ( 4 ( ) ∨ 6 ( ))} , where the concentration level (high or low) of the gene WNT5A is denoted by 1 , the concentration level of the gene pirin by 2 , the concentration level of the gene S100P is denoted by 3 , the concentration level of the gene RET1 is denoted by 4 , the concentration level of the gene MART1 is denoted by 5 , the concentration level of the gene HADHB is denoted by 6 , and the concentration level of the gene STC2 is denoted by 7 . See [14] for further details. In a WNT5A network, it is important to inhibit the concentration level of the gene WNT5A [27].
The optimal control problem is formulated. For simplicity, we consider only the structural control input. Then, suppose that the number of the structural control inputs is two. If 1 ( ) = 1, then the system is given as: The case of 1 ( ) = 1 corresponds to the situation such that the dynamics of the WNT5A network, that is, ( + 1) = ( ( )), are selected with high probability. The case of 2 ( ) = 1 corresponds to the situation such that the state is not changed; that is, ( + 1) = ( ) is selected with high probability. From the previous setting, 1,1 = 0.8, 1,2 = 0.2, From the obtained inputs, we see that the system is controlled by switching two discrete probability distributions, and the obtained inputs satisfy the constraint on duration of drug effectiveness. Noting that the trivial value of * is 15, we see that in this case the effectiveness of control synthesis is clear. Finally, we discuss the computation time for solving Problem 1. The computation time of the ILP problem was less than 20 [msec], where we used IBM ILOG CPLEX 11.0 as an ILP solver on the computer with Windows Vista 32-bit, the Intel Core 2 Duo CPU 3.0 GHz, and the 4 GB memory. Since the WNT5A network considered here is small size, Problem 1 can be solved fast.

Artificial Example.
In order to evaluate the computation time for solving Problem 1, we consider one artificial example of a BN with 15 states and 3 control inputs. We stress that the existing method [11,12,[17][18][19]21] cannot be applied to such a BN. This is because it is necessary to compute the state transition diagram such as that in Figure 2, that is, the transition probability matrix with 2 × 2 . In naive implementation using MATLAB [34], matrices with 2 15 × 2 15 cannot be created due to memory consumption, where we used the computer described previously.
Next, we discuss the computation time. Consider the two cases of = 10 and = 20. Then, in the ILP problem (Problem B) obtained, the dimension of binary variables is 1420 for = 10 and 2840 for = 20, and the number of inequalities is 3381 for = 10 and 6731 for = 20. In the case of = 10, the computation time of the ILP problem was 96 [sec], where we used the computer described previously. In the case of = 20, the computation time of the ILP problem was 238 [sec]. We remark that BNs with such a size are large scale in control problems of gene regulatory networks. Thus, we conclude that Problem 1 can be solved within the practical computation time.

Conclusions
In this paper, we have proposed a Boolean network (BN) model with two types of the control inputs and an optimal control method. By using this model, several situations in control of gene regulatory networks can be modeled. To model more realistic situations, duration of drug effectiveness has also been introduced. Since duration is given for each control input, effectiveness of multiple drugs can be evaluated. Furthermore, for this BN model, the optimal control problem has been formulated, and this problem is reduced to an integer linear programming problem. Finally, numerical simulations have been shown. The proposed method provides us with a basic in the control theory of gene regulatory networks.
Recently, to simplify state transition diagrams such as that in Figure 2, a stochastic Boolean network has been proposed in [35]. The authors proposed in [36] a similar method using polynomial optimization. In addition, to simplify a given Boolean network, the Karnaugh map realization of a Boolean network has been proposed in [37]. These methods are useful for reducing the computational burden. It is one of the future works to consider the control problem with duration of drug effectiveness based on these methods.