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

Authors: Koichi Kobayashi; Kunihiko Hiraishi
Journal: BioMed Research International
 (2013)

---

## Abstract

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.


---

## Body

## 1. 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–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–21]). In many existing results, state transition diagrams with 2n nodes (i.e., 2n×2n transition probability matrices) must be computed for a PBN with n 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 Section2.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 1.
Letℛ denote the set of real numbers. Let {0,1}n denote the set of n-dimensional vectors, which consists of elements 0 and 1. Let In and 0m×n denote the n×n identity matrix and  the m×n zero matrix, respectively. For simplicity, we sometimes use the symbol 0 instead of 0m×n and  and the symbol I instead of In. For a matrix M, lnM denotes the matrix such that the (i,j)th element is given as the natural logarithm of the (i,j)th element in M. For a matrix M, MT denotes the transpose of M.


## 2. Materials and Methods

### 2.1. The Boolean Networks with Control Inputs

A Boolean network (BN) withn states is given by
(1)x(k+1)=fa(x(k)),wherex∈{0,1}n is the state (e.g., the concentration of genes) and  k=0,1,2,… is the discrete time. The function fa:{0,1}n→{0,1}n is a given Boolean function with logical operators such as AND (∧), OR (∨), and NOT (¬). If the BN (1) is deterministic, then the next state x(k+1) is uniquely determined for a given x(k). See also Example 1 for an example.Next, the control inputs are added to a BN (1). For the BN (1) with n 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:
(2)x(k+1)=f(x(k),u(k)),whereu∈{0,1}m is the control input; that is, the value of u (e.g., the concentration of genes) can be arbitrarily given, and f:{0,1}n×{0,1}m→{0,1}n is a given Boolean function. The ith element of the state x and the ith element of the control input u are denoted by xi and ui, respectively. In the BN (2), x(k+1) is uniquely determined for the given x(k) and u(k).Then, consider the structural control input. Suppose that the candidates off are given by fi, i=1,2,…,l. It will be difficult to select one Boolean function uniquely. In this paper, we assume that one discrete probability distribution is selected among ms 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 ri,j denote the probability that the Boolean function fj is selected in the ith discrete probability distribution. Then,
(3)∑j=1l‍ri,j=1,i=1,2,…,ms,hold. In addition,ms-dimensional binary variables us∈{0,1}ms are assigned to ms discrete probability distributions, and let uis denote the ith element of us. The structural control input us corresponds to ms kinds of external stimuli. Then, the equality constraint
(4)∑i=1ms‍uis(k)=1is imposed. Here, we show a simple example.Example 1.
As a simple example, consider the simplified model of an apoptosis network in Figure1 [25]. Then, the Boolean network model expressing this apoptosis network is given by
(5)x1(k+1)=¬x2(k)∧u(k),x2(k+1)=¬x1(k)∧x3(k),x3(k+1)=x2(k)∨u(k),
where the concentration level (high or low) of the inhibitor of apoptosis proteins (IAPs) is denoted byx1, the concentration level of the active caspase 3 (C3a) is denoted by x2, and the concentration level of the active caspase 8 (C8a) is denoted  by x3. The concentration level of the tumor necrosis factor (TNF, a stimulus) is denoted by u 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, x2=1, implies cell near-death, otherwise, cell survival. As seen in (5), if the concentration of IAP is high (x1=1) or the concentration of the caspase C8a is low (x3=0), then the concentration of C3a becomes low; that is, x2=0. On the other hand, x1 and x3 at the next time depend on the value of x2 as well as u. In this way, some dynamical interactions exist. See [25, 26] for further details.
Suppose thatl=2 and ms=2. Then, as an example of the candidates of the Boolean functions, we consider the following:
(6)f1=[¬x2(k)∧u(k)¬x1(k)∧x3(k)x2(k)∨u(k)],r1,1=0.8,r2,1=0.1,f2=[x1(k)x2(k)x3(k)],r1,2=0.2,r2,2=0.9.
We suppose that the Boolean functionf1 expresses the situation that the dynamics of an apoptosis network are selected with high probability and that the Boolean function f2 expresses the situation that the state is not changed with high probability. By using u1s and u2s, one of the two discrete probability distributions {r1,1,r1,2} and {r2,1,r2,2} is selected at each time.Figure 1
A simplified model of an apoptosis network: activation (solid) and inhibition (broken).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:(7)x1(k+1)={x3(k)∨u(k),withtheprobability0.8,¬x3(k),withtheprobability0.2,x2(k+1)=x1(k)∧¬x3(k),withtheprobability1.0,x3(k+1)={x1(k)∧¬x2(k),withtheprobability0.7,x2(k)∨u(k),withtheprobability0.3.
This PBN corresponds to the cases ofl=4 and ms=1. The candidates of the Boolean functions fi, i=1,2,3,4, and the probabilities r1,j, j=1,2,3,4, are obtained as follows:
(8)f1=[x3(k)∨u(k)x1(k)∧¬x3(k)x1(k)∧¬x2(k)],r1,1=0.56,f2=[x3(k)∨u(k)x1(k)∧¬x3(k)x2(k)∨u(k)],r1,2=0.24,f3=[¬x3(k)x1(k)∧¬x3(k)x1(k)∧¬x2(k)],r1,3=0.14,f4=[¬x3(k)x1(k)∧¬x3(k)x2(k)∨u(k)],r1,4=0.06.
Next, consider the state orbit of this PBN. In PBNs, one Boolean function is probabilistically selected at each time. Then, forx(0)=[000]T and u(0)=0, we obtain the following:
(9)Prob(x(1)=[000]T∣x(0)=[000]T)=0.8,Prob(x(1)=[100]T∣x(0)=[000]T)=0.2.
In this example, the cardinality of the finite state set{0,1}3 is given by 23=8, and we obtain the state transition diagram of Figure 2 by computing the transition from each value of the state. In Figure 2, the number assigned to each node denotes x1, x2, and  x3 (each element of the state), and the number assigned to each arc denotes the transition probability from some state to another state. For simplicity of illustration, the state transitions from x(k)=[000]T,[001]T,[010]T,[110]T are illustrated in 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.Figure 2
The state transition diagram withu(k)=0.As shown in this example, computing state transition diagrams with2n nodes (n is the number of the state) is required in the existing solution methods for optimal control of PBNs with n 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 inputu can be eliminated from (2) by fixing the value of u in (2). Then, the number of the candidates of Boolean functions is 2ml, and 2m combinations for u must be computed in advance. To avoid this computation, we consider two types of the control inputs.


### 2.2. Duration of Drug Effectiveness

The control inputu and the structural control input us 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 Lui,Wui1, and Wui0 are defined for each input ui (or uis). The parameters Lui and Wui1 have been already defined in [15].The parameterLui is the length of the drug effectiveness period. That is, if ui(k)=1, then ui(k+1)=ui(k+2)=⋯=ui(k+Lui)=1 holds. Next, Wui1(>Lui) is explained. If ui(k)=1, then ui(k+1),ui(k+2),…,ui(k+Wui1-1) is uniquely determined depending on Lui, and ui(k+Wui1) is a decision variable. Then, Wui1-Lui-1 corresponds to the length of a recovery phase. Finally, Wui0 is explained. If ui(k)=0, then ui(k+1)=ui(k+2)=⋯=ui(k+Wui0-1)=0 holds, and ui(k+Wui0) is a decision variable. By using Lui,Wui1,Wui0, we can consider several situations, and we show two typical examples.Example 4.
First, suppose that, for the control inputu∈{0,1}1, Lu, Wu1, and Wu0 are given as Lu=1, Wu1=3, and Wu0=2, respectively. Consider the case of u(k)=1. Then, u(k+1) and u(k+2) are uniquely determined as u(k+1)=1 and u(k+2)=0, respectively, and u(k+3) is a decision variable. Then, u(k+2)=0 is the recovery phase, and Wu1-Lu-1=1 is its length. In the case of u(k)=0, the relation u(k+1)=0 holds, and u(k+2) is a decision variable.
Another example is shown. Suppose that, for the control inputu∈{0,1}1, Lu, Wu1, and Wu0 are given as Lu=0, Wu1=3, and Wu0=3, respectively. In both, the case of u(k)=1 and the case of u(k)=0, u(k+1)=u(k+2)=0 holds, and u(k+3) is a decision variable. In this case, u(k+1)=u(k+2)=0 is the recovery phase, and Wu1-Lu-1=2 is its length.By using the  three parametersLui,Wui1, and Wui0, several situations on the duration of effect can be modeled (see also [15]). In addition, since these parameters can be given for each ui (or uis), 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.


### 2.3. Optimal Control Problem

First, the following two notations are defined. Letπi(k) denote the probability that some Boolean function fi is selected at time k. In addition, the probability that some time sequence of theBoolean functions fi(k1),fi(k1+1),…,fi(k2) is selected at time interval [k1,k2] is denoted by
(10)π(k1,k2):=∏k=k1k2‍πi(k)(k).For simplicity of notation,i(k1),i(k1+1),…,i(k2) are omitted in π(k1,k2).Next, for the Boolean networks withn states and two types of the control inputs, consider the following optimal control problem.Problem 1.
Suppose that, for the Boolean network withn states and two types of the control inputs, the initial state x(0)=x0, ρ satisfying 0≤ρ≤1, the control time N, the parameters on duration of drug effectiveness Lui(uis),Wui(uis)1, and Wui(uis)0 are given. Then, for all combinations of the Boolean functions satisfying the  constraint
(11)π(0,N-1)≥ρ,
find two control input sequencesu(0),u(1),…,u(N-1) and us(0),us(1),…,us(N-1) minimizing the lower bound of the cost function
(12)J=∑k=0N-1‍{Qx(k)+Ru(k)+Rsus(k)}+Qfx(N)
subject to the constraint on duration of drug effectiveness, whereQ,Qf∈ℛ1×n, R∈ℛ1×m, and R∈ℛ1×ms are weighting vectors whose elements are nonnegative real numbers.For simplicity of discussion, a linear function with respect tox, u, and us 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 xi2(k) 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–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–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 Example1 again. For this system, we consider finding a control strategy such that a stimulus u is not applied as much as possible, and cell survival is achieved; u=0 implies that a stimulus is not applied to the system, and x1=1 and x2=0 express cell survival [25]. Then, as one of the appropriate cost functions, we can consider the following cost function:
(13)J=∑i=0N-1‍{10|x1(i)-1|+10|x2(i)-0|+u(i)}J=+100|x1(N)-1|+100|x2(N)-0|.
By the coordinate transformation ofx1 into 1-x1, this cost function can be rewritten as the form of (12).


### 2.4. Solution Method

We propose a solution method for Problem1. 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 Problem1 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:(i)
¬δ1 is equivalent to 1-δ1,(ii)
δ1∨δ2 is equivalent to δ1+δ2-δ1δ2,(iii)
δ1∧δ2 is equivalent to δ1δ2.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δj∈{0,1} and j∈𝒥 are given, where 𝒥 is some index set. Then, z=∏j∈𝒥‍δj is equivalent to the following linear inequalities:
(14)∑j∈𝒥‍δj-z≤|𝒥|-1,-∑j∈𝒥‍δj+|𝒥|z≤0,
where|𝒥| is the cardinality of 𝒥.From Lemmas6 and 7, we see that any Boolean function can be equivalently transformed into a pair of some linear function and some linear inequality. See [31, 32] for further details. For example, δ1∨¬δ2 is equivalent to a pair of 1-δ2+z and z=δ1δ2. By using Lemma 7, z=δ1δ2 can be expressed as a set of linear inequalities.Now, we consider reducing Problem1 to an ILP problem.By using Lemma6, the candidates of the Boolean functions fi(x(k),u(k)), i=1,2,…,l, are transformed into a polynomial on the real number field. Let f^i(x(k),u(k)) denote the polynomial obtained. Then, consider the following system using f^i(x(k),u(k)):
(15)x(k+1)=∑i=1l‍{δi(k)f^i(x(k),u(k))},whereδ1(k),δ2(k),…,δl(k) are binary variables satisfying
(16)∑i=1l‍δi(k)=1.The binary vectorδ(k):=[δ1(k)δ2(k)⋯δl(k)]T is used to select the polynomial f^i and to express (11) as a linear form. Here, we define the following vector:
(17)Si:=[ri,1ri,2⋯ri,l].Then, by using the natural logarithm,π(0,N-1) in (11) is expressed as
(18)lnπ(0,N-1)=∑k=0N-1‍(∑i=1ms‍lnSiuis(k))δ(k).In this expression, one probability distribution is selected by usinguis(k), and the probability that a certain Boolean function is selected is determined by δ(k). Then, Problem 1 is equivalent to the following problem.Problem A
(19)
find
u
(
k
)
,
u
s
(
k
)
,
δ
(
k
)
,
k
=
0,1
,
…
,
N
-
1
,
min
Cost function
(
12
)
,
subject  to
System
(
15
)
,
x
(
0
)
=
x
0
,
Inequality constraint:
∑
k
=
0
N
-
1
‍
(
∑
i
=
1
m
s
‍
ln
S
i
δ
i
s
(
k
)
)
δ
(
k
)
≥
ln
ρ
,
Equality
constraint
(
4
)
,
(
16
)
,
Constraint  on  duration  of  drug  effectiveness.By using Lemma7, the system (15) and (∑i=1ms‍lnSiδis(k))δ(k) can be equivalently expressed in the following linear form:
(20)x(k+1)=Ax(k)+Buu(k)+Bsus(k)+Bbzb(k),(21)(∑i=1ms‍lnSiδis(k))δ(k)=∑i=1ms‍lnSizis(k),(22)Ex(k)+Fuu(k)+Fsus(k)+Fzz(k)≤G,wherezis(k):=δis(k)δ(k), and (22) is the linear inequality obtained by applying Lemma 7 to (15). The vector zb(k)∈{0,1}b is an auxiliary binary variable obtained by using Lemma 7, and the dimension of zb(k), that is, b, is determined depending on the form of the given Boolean functions. In addition, z(k) is defined as
(23)z(k)≔[(zb(k))T(z1s(k))T(z2s(k))T⋯(zmss(k))T]Tz(k)∈{0,1}b+msl.Here in after, for simplicity of notation,Bbzb(k) is rewritten as Bzz(k), Bz:=[Bb0], and ∑i=1ms‍lnSizis(k) is rewritten as Cz(k), C:=[0lnS1lnS2⋯lnSms].Now, we consider transforming Problem A by using (20), (21), and (22). By using
(24)x(k)=Akx0+∑i=1k‍Ai-1(Buu(k-i)+Bsus(k-i)+Bzz(k))obtained from the state equation in (20), we can obtain
(25)x¯=A¯x0+B¯uu¯+B¯su¯s+B¯zz¯,where(26)x¯:=[(x(0))T(x(1))T⋯(x(N))T]T,u¯:=[(u(0))T(u(1))T⋯(u(N-1))T]T,u¯s:=[(us(0))T(us(1))T⋯(us(N-1))T]T,z¯:=[(z(0))T(z(1))T⋯(z(N-1))T]T,A¯=[IAA2⋮AN],B¯=[00⋯0I0⋯0A⋱⋱⋮⋮⋱⋱0AN-1⋯AI],B¯u=B¯[Bu0⋱0Bu],B¯s=B¯[Bs0⋱0Bs],B¯z=B¯[Bz0⋱0Bz].Next, the inequality constraint∑k=0N-1‍(∑i=1ms‍lnSiδis(k))δ(k)≥lnρ in Problem A is equivalent to
(27)-C¯z¯≤-lnρ,whereC¯=[CC⋯C]. Furthermore, from (22), we can obtain
(28)E¯x¯+F¯uu¯+F¯su¯s+F¯zz¯≤G¯,where(29)E¯=[E00⋱⋮0E0],F¯u=[Fu0⋱0Fu],F¯s=[Fs0⋱0Fs],F¯z=[Fz0⋱0Fz],G¯=[GG⋮G].Next, consider the constraint on duration of drug effectiveness. This constraint can be expressed as a Boolean function. Then, by using Lemmas6 and 7, it can be transformed into the following form:
(30)u¯=V1v¯+V2,W1v¯≤W2,(31)u¯s=Vs1v¯s+Vs2,Ws1v¯s≤Ws2,wherev¯ and v¯s are binary decision variables with certain dimensions. Deriving a general form of coefficient matrices will be difficult, but for the given Lui, Wui1, and Wui0, deriving coefficient matrices is easy.We show two examples.Example 8.
Consider Example4 again. First, consider the case of Lu=1, Wu1=3, and Wu0=2. Then, noting explanations in Example 4, we can obtain
(32)u(0)=v0,u(1)=v0,u(2)=(1-v0)v2,u(3)=(1-v0)v2+v0v3,u(4)=v0v3+(1-v2)v4,u(5)=(1-v2)v4+v2v5,u(5)⋮
and in this case, these are equivalent to(33)u(0)=v0,u(1)=v0,u(2)=v2,v2≤1-v0,u(3)=v2+v3,v2≤1-v0,v3≤v0,u(4)=v3+v4,v3≤v0,v4≤1-v2,0≤v3+v4≤1,u(5)=v4+v5,v4≤1-v2,v2≤v5,0≤v4+v5≤1,⋮
We explainu(2)=v2, and v2≤1-v0 as an example. If v0=1, then v2≤0 holds. Since v2 is binary, we can obtain v2=0; that is, u(2)=0. If v0=0, then v2≤1 holds, and we can obtain u(2)=v2. That is, u(2) can take on either 0 or 1. From the previous discussion, we see that a pair of u(2)=v2 and v2≤1-v0 is equivalent to u(2)=(1-v0)v2. Thus, we can obtain the forms of (30) and (31). In the case of N=5 (N is the control time in Problem 1), (30) can be obtained as
(34)[u(0)u(1)u(2)u(3)u(4)]︸u¯=[10001000010001100011]︸V1[v0v2v3v4]︸v¯,[1100-1010010100-1-10011]︸W1[v0v2v3v4]︸v¯≤[10101]︸W2,
whereV2=0. We remark that, in a general case, the product such as z=v0v2 must be transformed into linear inequalities by using Lemma 7.
Next, consider the case ofLu=0, Wu1=3, and Wu0=3. Then, we can obtain the following:
(35)u(0)=v0,u(1)=0,u(2)=0,u(3)=v3,u(4)=0,u(5)=0,u(5)⋮
and this case is one of the simplest cases. In the case ofN=5 (N is the control time in Problem 1), (30) can be obtained as
(36)[u(0)u(1)u(2)u(3)u(4)]︸u¯=[1000000100]︸V1[v0v3]︸v¯,
whereV2=0, W1=0, and W2=0.Finally, the cost function (12) is rewritten as
(37)J=Q¯x¯+R¯u¯+R¯su¯s,whereQ¯=[Q⋯QQf], R¯=[R⋯R], and R¯s=[Rs⋯Rs]. By substituting (25), (30), and (31) into (28) and (37), Problem A is equivalent to the following ILP problem.Problem B
(38)
find
v
¯
,
v
¯
s
,
z
¯
,
min
(
R
¯
+
Q
¯
B
¯
u
)
V
1
v
¯
+
(
R
¯
s
+
Q
¯
B
¯
s
)
V
s
1
v
¯
s
+
Q
¯
B
¯
z
z
¯
+
Q
¯
A
¯
x
0
,
subject  to
-
C
¯
z
¯
≤
-
ln
ρ
(
F
¯
u

+

E
¯
B
¯
u
)
V
1
v
¯

+

(
F
¯
s

+

E
¯
B
¯
s
)
V
s
1
v
¯
s

+

(
F
¯
z
+
E
¯
B
¯
z
)
z
¯
≤
G
¯
-
F
¯
u
V
1
-
F
¯
s
V
s
1
-
E
¯
(
A
¯
x
0
+
B
¯
u
V
2
+
B
¯
s
V
s
2
)
,
W
1
v
¯
≤
W
2
,
W
s
1
v
¯
s
≤
W
s
2
.Problem B can be solved by a suitable solver such as IBM ILOG CPLEX [33].


## 2.1. The Boolean Networks with Control Inputs

A Boolean network (BN) withn states is given by
(1)x(k+1)=fa(x(k)),wherex∈{0,1}n is the state (e.g., the concentration of genes) and  k=0,1,2,… is the discrete time. The function fa:{0,1}n→{0,1}n is a given Boolean function with logical operators such as AND (∧), OR (∨), and NOT (¬). If the BN (1) is deterministic, then the next state x(k+1) is uniquely determined for a given x(k). See also Example 1 for an example.Next, the control inputs are added to a BN (1). For the BN (1) with n 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:
(2)x(k+1)=f(x(k),u(k)),whereu∈{0,1}m is the control input; that is, the value of u (e.g., the concentration of genes) can be arbitrarily given, and f:{0,1}n×{0,1}m→{0,1}n is a given Boolean function. The ith element of the state x and the ith element of the control input u are denoted by xi and ui, respectively. In the BN (2), x(k+1) is uniquely determined for the given x(k) and u(k).Then, consider the structural control input. Suppose that the candidates off are given by fi, i=1,2,…,l. It will be difficult to select one Boolean function uniquely. In this paper, we assume that one discrete probability distribution is selected among ms 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 ri,j denote the probability that the Boolean function fj is selected in the ith discrete probability distribution. Then,
(3)∑j=1l‍ri,j=1,i=1,2,…,ms,hold. In addition,ms-dimensional binary variables us∈{0,1}ms are assigned to ms discrete probability distributions, and let uis denote the ith element of us. The structural control input us corresponds to ms kinds of external stimuli. Then, the equality constraint
(4)∑i=1ms‍uis(k)=1is imposed. Here, we show a simple example.Example 1.
As a simple example, consider the simplified model of an apoptosis network in Figure1 [25]. Then, the Boolean network model expressing this apoptosis network is given by
(5)x1(k+1)=¬x2(k)∧u(k),x2(k+1)=¬x1(k)∧x3(k),x3(k+1)=x2(k)∨u(k),
where the concentration level (high or low) of the inhibitor of apoptosis proteins (IAPs) is denoted byx1, the concentration level of the active caspase 3 (C3a) is denoted by x2, and the concentration level of the active caspase 8 (C8a) is denoted  by x3. The concentration level of the tumor necrosis factor (TNF, a stimulus) is denoted by u 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, x2=1, implies cell near-death, otherwise, cell survival. As seen in (5), if the concentration of IAP is high (x1=1) or the concentration of the caspase C8a is low (x3=0), then the concentration of C3a becomes low; that is, x2=0. On the other hand, x1 and x3 at the next time depend on the value of x2 as well as u. In this way, some dynamical interactions exist. See [25, 26] for further details.
Suppose thatl=2 and ms=2. Then, as an example of the candidates of the Boolean functions, we consider the following:
(6)f1=[¬x2(k)∧u(k)¬x1(k)∧x3(k)x2(k)∨u(k)],r1,1=0.8,r2,1=0.1,f2=[x1(k)x2(k)x3(k)],r1,2=0.2,r2,2=0.9.
We suppose that the Boolean functionf1 expresses the situation that the dynamics of an apoptosis network are selected with high probability and that the Boolean function f2 expresses the situation that the state is not changed with high probability. By using u1s and u2s, one of the two discrete probability distributions {r1,1,r1,2} and {r2,1,r2,2} is selected at each time.Figure 1
A simplified model of an apoptosis network: activation (solid) and inhibition (broken).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:(7)x1(k+1)={x3(k)∨u(k),withtheprobability0.8,¬x3(k),withtheprobability0.2,x2(k+1)=x1(k)∧¬x3(k),withtheprobability1.0,x3(k+1)={x1(k)∧¬x2(k),withtheprobability0.7,x2(k)∨u(k),withtheprobability0.3.
This PBN corresponds to the cases ofl=4 and ms=1. The candidates of the Boolean functions fi, i=1,2,3,4, and the probabilities r1,j, j=1,2,3,4, are obtained as follows:
(8)f1=[x3(k)∨u(k)x1(k)∧¬x3(k)x1(k)∧¬x2(k)],r1,1=0.56,f2=[x3(k)∨u(k)x1(k)∧¬x3(k)x2(k)∨u(k)],r1,2=0.24,f3=[¬x3(k)x1(k)∧¬x3(k)x1(k)∧¬x2(k)],r1,3=0.14,f4=[¬x3(k)x1(k)∧¬x3(k)x2(k)∨u(k)],r1,4=0.06.
Next, consider the state orbit of this PBN. In PBNs, one Boolean function is probabilistically selected at each time. Then, forx(0)=[000]T and u(0)=0, we obtain the following:
(9)Prob(x(1)=[000]T∣x(0)=[000]T)=0.8,Prob(x(1)=[100]T∣x(0)=[000]T)=0.2.
In this example, the cardinality of the finite state set{0,1}3 is given by 23=8, and we obtain the state transition diagram of Figure 2 by computing the transition from each value of the state. In Figure 2, the number assigned to each node denotes x1, x2, and  x3 (each element of the state), and the number assigned to each arc denotes the transition probability from some state to another state. For simplicity of illustration, the state transitions from x(k)=[000]T,[001]T,[010]T,[110]T are illustrated in 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.Figure 2
The state transition diagram withu(k)=0.As shown in this example, computing state transition diagrams with2n nodes (n is the number of the state) is required in the existing solution methods for optimal control of PBNs with n 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 inputu can be eliminated from (2) by fixing the value of u in (2). Then, the number of the candidates of Boolean functions is 2ml, and 2m combinations for u must be computed in advance. To avoid this computation, we consider two types of the control inputs.


## 2.2. Duration of Drug Effectiveness

The control inputu and the structural control input us 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 Lui,Wui1, and Wui0 are defined for each input ui (or uis). The parameters Lui and Wui1 have been already defined in [15].The parameterLui is the length of the drug effectiveness period. That is, if ui(k)=1, then ui(k+1)=ui(k+2)=⋯=ui(k+Lui)=1 holds. Next, Wui1(>Lui) is explained. If ui(k)=1, then ui(k+1),ui(k+2),…,ui(k+Wui1-1) is uniquely determined depending on Lui, and ui(k+Wui1) is a decision variable. Then, Wui1-Lui-1 corresponds to the length of a recovery phase. Finally, Wui0 is explained. If ui(k)=0, then ui(k+1)=ui(k+2)=⋯=ui(k+Wui0-1)=0 holds, and ui(k+Wui0) is a decision variable. By using Lui,Wui1,Wui0, we can consider several situations, and we show two typical examples.Example 4.
First, suppose that, for the control inputu∈{0,1}1, Lu, Wu1, and Wu0 are given as Lu=1, Wu1=3, and Wu0=2, respectively. Consider the case of u(k)=1. Then, u(k+1) and u(k+2) are uniquely determined as u(k+1)=1 and u(k+2)=0, respectively, and u(k+3) is a decision variable. Then, u(k+2)=0 is the recovery phase, and Wu1-Lu-1=1 is its length. In the case of u(k)=0, the relation u(k+1)=0 holds, and u(k+2) is a decision variable.
Another example is shown. Suppose that, for the control inputu∈{0,1}1, Lu, Wu1, and Wu0 are given as Lu=0, Wu1=3, and Wu0=3, respectively. In both, the case of u(k)=1 and the case of u(k)=0, u(k+1)=u(k+2)=0 holds, and u(k+3) is a decision variable. In this case, u(k+1)=u(k+2)=0 is the recovery phase, and Wu1-Lu-1=2 is its length.By using the  three parametersLui,Wui1, and Wui0, several situations on the duration of effect can be modeled (see also [15]). In addition, since these parameters can be given for each ui (or uis), 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.


## 2.3. Optimal Control Problem

First, the following two notations are defined. Letπi(k) denote the probability that some Boolean function fi is selected at time k. In addition, the probability that some time sequence of theBoolean functions fi(k1),fi(k1+1),…,fi(k2) is selected at time interval [k1,k2] is denoted by
(10)π(k1,k2):=∏k=k1k2‍πi(k)(k).For simplicity of notation,i(k1),i(k1+1),…,i(k2) are omitted in π(k1,k2).Next, for the Boolean networks withn states and two types of the control inputs, consider the following optimal control problem.Problem 1.
Suppose that, for the Boolean network withn states and two types of the control inputs, the initial state x(0)=x0, ρ satisfying 0≤ρ≤1, the control time N, the parameters on duration of drug effectiveness Lui(uis),Wui(uis)1, and Wui(uis)0 are given. Then, for all combinations of the Boolean functions satisfying the  constraint
(11)π(0,N-1)≥ρ,
find two control input sequencesu(0),u(1),…,u(N-1) and us(0),us(1),…,us(N-1) minimizing the lower bound of the cost function
(12)J=∑k=0N-1‍{Qx(k)+Ru(k)+Rsus(k)}+Qfx(N)
subject to the constraint on duration of drug effectiveness, whereQ,Qf∈ℛ1×n, R∈ℛ1×m, and R∈ℛ1×ms are weighting vectors whose elements are nonnegative real numbers.For simplicity of discussion, a linear function with respect tox, u, and us 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 xi2(k) 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–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–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 Example1 again. For this system, we consider finding a control strategy such that a stimulus u is not applied as much as possible, and cell survival is achieved; u=0 implies that a stimulus is not applied to the system, and x1=1 and x2=0 express cell survival [25]. Then, as one of the appropriate cost functions, we can consider the following cost function:
(13)J=∑i=0N-1‍{10|x1(i)-1|+10|x2(i)-0|+u(i)}J=+100|x1(N)-1|+100|x2(N)-0|.
By the coordinate transformation ofx1 into 1-x1, this cost function can be rewritten as the form of (12).


## 2.4. Solution Method

We propose a solution method for Problem1. 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 Problem1 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:(i)
¬δ1 is equivalent to 1-δ1,(ii)
δ1∨δ2 is equivalent to δ1+δ2-δ1δ2,(iii)
δ1∧δ2 is equivalent to δ1δ2.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δj∈{0,1} and j∈𝒥 are given, where 𝒥 is some index set. Then, z=∏j∈𝒥‍δj is equivalent to the following linear inequalities:
(14)∑j∈𝒥‍δj-z≤|𝒥|-1,-∑j∈𝒥‍δj+|𝒥|z≤0,
where|𝒥| is the cardinality of 𝒥.From Lemmas6 and 7, we see that any Boolean function can be equivalently transformed into a pair of some linear function and some linear inequality. See [31, 32] for further details. For example, δ1∨¬δ2 is equivalent to a pair of 1-δ2+z and z=δ1δ2. By using Lemma 7, z=δ1δ2 can be expressed as a set of linear inequalities.Now, we consider reducing Problem1 to an ILP problem.By using Lemma6, the candidates of the Boolean functions fi(x(k),u(k)), i=1,2,…,l, are transformed into a polynomial on the real number field. Let f^i(x(k),u(k)) denote the polynomial obtained. Then, consider the following system using f^i(x(k),u(k)):
(15)x(k+1)=∑i=1l‍{δi(k)f^i(x(k),u(k))},whereδ1(k),δ2(k),…,δl(k) are binary variables satisfying
(16)∑i=1l‍δi(k)=1.The binary vectorδ(k):=[δ1(k)δ2(k)⋯δl(k)]T is used to select the polynomial f^i and to express (11) as a linear form. Here, we define the following vector:
(17)Si:=[ri,1ri,2⋯ri,l].Then, by using the natural logarithm,π(0,N-1) in (11) is expressed as
(18)lnπ(0,N-1)=∑k=0N-1‍(∑i=1ms‍lnSiuis(k))δ(k).In this expression, one probability distribution is selected by usinguis(k), and the probability that a certain Boolean function is selected is determined by δ(k). Then, Problem 1 is equivalent to the following problem.Problem A
(19)
find
u
(
k
)
,
u
s
(
k
)
,
δ
(
k
)
,
k
=
0,1
,
…
,
N
-
1
,
min
Cost function
(
12
)
,
subject  to
System
(
15
)
,
x
(
0
)
=
x
0
,
Inequality constraint:
∑
k
=
0
N
-
1
‍
(
∑
i
=
1
m
s
‍
ln
S
i
δ
i
s
(
k
)
)
δ
(
k
)
≥
ln
ρ
,
Equality
constraint
(
4
)
,
(
16
)
,
Constraint  on  duration  of  drug  effectiveness.By using Lemma7, the system (15) and (∑i=1ms‍lnSiδis(k))δ(k) can be equivalently expressed in the following linear form:
(20)x(k+1)=Ax(k)+Buu(k)+Bsus(k)+Bbzb(k),(21)(∑i=1ms‍lnSiδis(k))δ(k)=∑i=1ms‍lnSizis(k),(22)Ex(k)+Fuu(k)+Fsus(k)+Fzz(k)≤G,wherezis(k):=δis(k)δ(k), and (22) is the linear inequality obtained by applying Lemma 7 to (15). The vector zb(k)∈{0,1}b is an auxiliary binary variable obtained by using Lemma 7, and the dimension of zb(k), that is, b, is determined depending on the form of the given Boolean functions. In addition, z(k) is defined as
(23)z(k)≔[(zb(k))T(z1s(k))T(z2s(k))T⋯(zmss(k))T]Tz(k)∈{0,1}b+msl.Here in after, for simplicity of notation,Bbzb(k) is rewritten as Bzz(k), Bz:=[Bb0], and ∑i=1ms‍lnSizis(k) is rewritten as Cz(k), C:=[0lnS1lnS2⋯lnSms].Now, we consider transforming Problem A by using (20), (21), and (22). By using
(24)x(k)=Akx0+∑i=1k‍Ai-1(Buu(k-i)+Bsus(k-i)+Bzz(k))obtained from the state equation in (20), we can obtain
(25)x¯=A¯x0+B¯uu¯+B¯su¯s+B¯zz¯,where(26)x¯:=[(x(0))T(x(1))T⋯(x(N))T]T,u¯:=[(u(0))T(u(1))T⋯(u(N-1))T]T,u¯s:=[(us(0))T(us(1))T⋯(us(N-1))T]T,z¯:=[(z(0))T(z(1))T⋯(z(N-1))T]T,A¯=[IAA2⋮AN],B¯=[00⋯0I0⋯0A⋱⋱⋮⋮⋱⋱0AN-1⋯AI],B¯u=B¯[Bu0⋱0Bu],B¯s=B¯[Bs0⋱0Bs],B¯z=B¯[Bz0⋱0Bz].Next, the inequality constraint∑k=0N-1‍(∑i=1ms‍lnSiδis(k))δ(k)≥lnρ in Problem A is equivalent to
(27)-C¯z¯≤-lnρ,whereC¯=[CC⋯C]. Furthermore, from (22), we can obtain
(28)E¯x¯+F¯uu¯+F¯su¯s+F¯zz¯≤G¯,where(29)E¯=[E00⋱⋮0E0],F¯u=[Fu0⋱0Fu],F¯s=[Fs0⋱0Fs],F¯z=[Fz0⋱0Fz],G¯=[GG⋮G].Next, consider the constraint on duration of drug effectiveness. This constraint can be expressed as a Boolean function. Then, by using Lemmas6 and 7, it can be transformed into the following form:
(30)u¯=V1v¯+V2,W1v¯≤W2,(31)u¯s=Vs1v¯s+Vs2,Ws1v¯s≤Ws2,wherev¯ and v¯s are binary decision variables with certain dimensions. Deriving a general form of coefficient matrices will be difficult, but for the given Lui, Wui1, and Wui0, deriving coefficient matrices is easy.We show two examples.Example 8.
Consider Example4 again. First, consider the case of Lu=1, Wu1=3, and Wu0=2. Then, noting explanations in Example 4, we can obtain
(32)u(0)=v0,u(1)=v0,u(2)=(1-v0)v2,u(3)=(1-v0)v2+v0v3,u(4)=v0v3+(1-v2)v4,u(5)=(1-v2)v4+v2v5,u(5)⋮
and in this case, these are equivalent to(33)u(0)=v0,u(1)=v0,u(2)=v2,v2≤1-v0,u(3)=v2+v3,v2≤1-v0,v3≤v0,u(4)=v3+v4,v3≤v0,v4≤1-v2,0≤v3+v4≤1,u(5)=v4+v5,v4≤1-v2,v2≤v5,0≤v4+v5≤1,⋮
We explainu(2)=v2, and v2≤1-v0 as an example. If v0=1, then v2≤0 holds. Since v2 is binary, we can obtain v2=0; that is, u(2)=0. If v0=0, then v2≤1 holds, and we can obtain u(2)=v2. That is, u(2) can take on either 0 or 1. From the previous discussion, we see that a pair of u(2)=v2 and v2≤1-v0 is equivalent to u(2)=(1-v0)v2. Thus, we can obtain the forms of (30) and (31). In the case of N=5 (N is the control time in Problem 1), (30) can be obtained as
(34)[u(0)u(1)u(2)u(3)u(4)]︸u¯=[10001000010001100011]︸V1[v0v2v3v4]︸v¯,[1100-1010010100-1-10011]︸W1[v0v2v3v4]︸v¯≤[10101]︸W2,
whereV2=0. We remark that, in a general case, the product such as z=v0v2 must be transformed into linear inequalities by using Lemma 7.
Next, consider the case ofLu=0, Wu1=3, and Wu0=3. Then, we can obtain the following:
(35)u(0)=v0,u(1)=0,u(2)=0,u(3)=v3,u(4)=0,u(5)=0,u(5)⋮
and this case is one of the simplest cases. In the case ofN=5 (N is the control time in Problem 1), (30) can be obtained as
(36)[u(0)u(1)u(2)u(3)u(4)]︸u¯=[1000000100]︸V1[v0v3]︸v¯,
whereV2=0, W1=0, and W2=0.Finally, the cost function (12) is rewritten as
(37)J=Q¯x¯+R¯u¯+R¯su¯s,whereQ¯=[Q⋯QQf], R¯=[R⋯R], and R¯s=[Rs⋯Rs]. By substituting (25), (30), and (31) into (28) and (37), Problem A is equivalent to the following ILP problem.Problem B
(38)
find
v
¯
,
v
¯
s
,
z
¯
,
min
(
R
¯
+
Q
¯
B
¯
u
)
V
1
v
¯
+
(
R
¯
s
+
Q
¯
B
¯
s
)
V
s
1
v
¯
s
+
Q
¯
B
¯
z
z
¯
+
Q
¯
A
¯
x
0
,
subject  to
-
C
¯
z
¯
≤
-
ln
ρ
(
F
¯
u

+

E
¯
B
¯
u
)
V
1
v
¯

+

(
F
¯
s

+

E
¯
B
¯
s
)
V
s
1
v
¯
s

+

(
F
¯
z
+
E
¯
B
¯
z
)
z
¯
≤
G
¯
-
F
¯
u
V
1
-
F
¯
s
V
s
1
-
E
¯
(
A
¯
x
0
+
B
¯
u
V
2
+
B
¯
s
V
s
2
)
,
W
1
v
¯
≤
W
2
,
W
s
1
v
¯
s
≤
W
s
2
.Problem B can be solved by a suitable solver such as IBM ILOG CPLEX [33].


## 3. 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 x(k+1)=fa(x(k)) of the WNT5A network is given by the following:
(39)x1(k+1)=¬x6(k),x2(k+1)=(¬x2(k)∧x4(k)∧x6(k))x2(k+1)∨{¬x2(k)∧(x4(k)∨x6(k))},x3(k+1)=¬x7(k),x4(k+1)=x4(k),x5(k+1)=x2(k)∨¬x7(k),x6(k+1)=x3(k)∨x4(k),x7(k+1)=¬x2(k)∨x7(k),where the concentration level (high or low) of the gene WNT5A is denoted byx1, the concentration level of the gene pirin by x2, the concentration level of the gene S100P is denoted by x3, the concentration level of the gene RET1 is denoted  by x4, the concentration level of the gene MART1 is denoted by x5, the concentration level of the gene HADHB is denoted by x6, and the concentration level of the gene STC2 is denoted by x7. 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. Ifu1s(k)=1, then the system is given as:
(40)x(k+1)={fa(x(k))withtheprobability0.8,x(k)withtheprobability0.2.Ifu2s(k)=1, then the system is given as:
(41)x(k+1)={fa(x(k))withtheprobability0.1,x(k)withtheprobability0.9.The case ofu1s(k)=1 corresponds to the situation such that the dynamics of the WNT5A network, that is, x(k+1)=fa(x(k)), are selected with high probability. The case of u2s(k)=1 corresponds to the situation such that the state is not changed; that is, x(k+1)=x(k) is selected with high probability. From the previous setting, r1,1=0.8, r1,2=0.2, r2,1=0.1, and r2,2=0.9. For this WNT5A network with structural control inputs, consider solving Problem 1. Q,Qf, and Rs in Problem 1 are given as Q=[1000000], Qf=[10000000], and Rs=[00], respectively. The initial state is given as x0=[1010100]T. In addition, the control time N in Problem 1 is given by N=5. Finally, the constraint on duration of drug effectiveness is imposed for only u2s(k). The parameters Lu2s, Wu2s1, and Wu2s0 are given as Lu2s=1, Wu2s1=3, and Wu2s0=2, respectively. Thus, we can obtain the ILP problem (Problem B), where the dimension of binary variables is 130 and the number of inequalities is 264.We show the computational result. LetJ_* denote the optimal value of the lower bound of a given cost function in Problem 1. Let J¯* denote the upper bound of the cost function derived by using the optimal control input. First, consider the case of ρ=10-5. Then, we can obtain J_*=2 and J¯*=15, and us(k) is obtained as
(42)us(0)=⋯=us(4)=[10].Noting thatr2,1=0.1 and ρ=10-5(=0.15), all combinations of the Boolean functions are considered, and the value of ρ is not appropriate. In particular, J¯*=15 implies that  x1(k)=1, k=0,1,…,5, and is the trivial upper bound.Next, consider the case ofρ=0.2. Then, we can obtain J_*=J¯*=4, and us(k) is obtained as
(43)us(0)=us(1)=[01],us(2)=us(3)=us(4)=[10].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 ofJ¯* is 15, we see that in this case the effectiveness of control synthesis is clear.Finally, we discuss the computation time for solving Problem1. 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.


### 3.2. Artificial Example

In order to evaluate the computation time for solving Problem1, we consider one artificial example of a BN with 15 states and 3 control inputs. We stress that the existing method [11, 12, 17–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 2n×2n. In naive implementation using MATLAB [34], matrices with 215×215 cannot be created due to memory consumption, where we used the computer described previously.The optimal control problem is formulated. In this example, we consider 3 control inputs and 2 structural control inputs. Ifu1s(k)=1, then the system is given as follows:
(44)x(k+1)={f1(x(k),u(k))withtheprobability0.8,f2(x(k),u(k))withtheprobability0.2.Ifu2s(k)=1, then the system is given as follows:
(45)x(k+1)={f1(x(k),u(k))withtheprobability0.2,f2(x(k),u(k))withtheprobability0.8.The Boolean functionf1 is given by the following:
(46)x1(k+1)=x1(k)∧¬x6(k)∨u3(k),x2(k+1)=¬x4(k)∧u1(k)∨u3(k),x3(k+1)=x5(k)∧u1(k)∨¬x10(k)∧x12(k)∧u3(k),x4(k+1)=x2(k)∧x5(k)∧¬u1(k)∨¬x14(k),x5(k+1)=¬u1(k)∧x6(k)∧x7(k)∨x12(k)∧x14(k),x6(k+1)=x1(k)∧x6(k)∧x10(k)∨¬x15(k),x7(k+1)=x6(k)∧x7(k)∧x8(k)∨u2(k)∧¬u3(k),x8(k+1)=x5(k)∧¬u1(k)∨x10(k)∧u2(k)∧x13(k),x9(k+1)=x3(k)∧u1(k)∨¬x8(k)∧x11(k),x10(k+1)=x6(k),x11(k+1)=x6(k)∧x10(k)∨¬u2(k)∧u3(k),x12(k+1)=x12(k)∧¬x15(k),x13(k+1)=¬u3(k),x14(k+1)=¬x14(k)∧u3(k),x15(k+1)=x14(k)∧x15(k).The Boolean functionf2 is given by the following:
(47)x1(k+1)=x2(k)∧x4(k)∧¬x8(k),x2(k+1)=¬x2(k)∧x3(k)∨u3(k),x3(k+1)=x1(k)∨¬x2(k)∧x3(k)∧x4(k),x4(k+1)=¬x1(k)∧x2(k)∧u1(k)∨x14(k),x5(k+1)=¬u2(k)∧x13(k)∧x14(k)∧x15(k),x6(k+1)=¬x2(k)∨x5(k)∧u1(k)∧¬u3(k),x7(k+1)=u1(k)∧x13(k),x8(k+1)=x5(k)∧x13(k),x9(k+1)=¬x6(k),x10(k+1)=x2(k)∧u2(k)∧¬x12(k)∧u3(k),x11(k+1)=¬x5(k)∧u1(k)∨¬x15(k)∧u3(k),x12(k+1)=x7(k)∧x3(k),x13(k+1)=x7(k)∧u2(k)∧u3(k),x14(k+1)=x12(k)∨x14(k)∧u3(k),x15(k+1)=x8(k).From the previous setting,r1,1=0.8, r1,2=0.2, r2,1=0.2, r2,2=0.8 hold. In Problem 1, Q,Qf,R, and Rs are given as Q=[1⋯1], Qf=[10⋯10], R=[101], and  Rs=[00], respectively. The initial state and the parameter ρ are given as x0=[1⋯1]T and ρ=10-4, respectively. The constraint on duration of drug effectiveness is imposed for u1(k) and u2s(k). For u1, the parameters Lu1, Wu11, and Wu10 are given as Lu1=0, Wu11=3, and Wu10=3, respectively. For u2s, the parameters Lu2s, Wu2s1, and Wu2s0 are given as Lu2s=1, Wu2s1=3, and Wu2s0=2, respectively.Next, we discuss the computation time. Consider the two cases ofN=10 and N=20. Then, in the ILP problem (Problem B) obtained, the dimension of binary variables is 1420 for N=10 and 2840 for N=20, and the number of inequalities is 3381 for N=10 and 6731 for N=20. In the case of N=10, the computation time of the ILP problem was 96 [sec], where we used the computer described previously. In the case of N=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.


## 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 x(k+1)=fa(x(k)) of the WNT5A network is given by the following:
(39)x1(k+1)=¬x6(k),x2(k+1)=(¬x2(k)∧x4(k)∧x6(k))x2(k+1)∨{¬x2(k)∧(x4(k)∨x6(k))},x3(k+1)=¬x7(k),x4(k+1)=x4(k),x5(k+1)=x2(k)∨¬x7(k),x6(k+1)=x3(k)∨x4(k),x7(k+1)=¬x2(k)∨x7(k),where the concentration level (high or low) of the gene WNT5A is denoted byx1, the concentration level of the gene pirin by x2, the concentration level of the gene S100P is denoted by x3, the concentration level of the gene RET1 is denoted  by x4, the concentration level of the gene MART1 is denoted by x5, the concentration level of the gene HADHB is denoted by x6, and the concentration level of the gene STC2 is denoted by x7. 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. Ifu1s(k)=1, then the system is given as:
(40)x(k+1)={fa(x(k))withtheprobability0.8,x(k)withtheprobability0.2.Ifu2s(k)=1, then the system is given as:
(41)x(k+1)={fa(x(k))withtheprobability0.1,x(k)withtheprobability0.9.The case ofu1s(k)=1 corresponds to the situation such that the dynamics of the WNT5A network, that is, x(k+1)=fa(x(k)), are selected with high probability. The case of u2s(k)=1 corresponds to the situation such that the state is not changed; that is, x(k+1)=x(k) is selected with high probability. From the previous setting, r1,1=0.8, r1,2=0.2, r2,1=0.1, and r2,2=0.9. For this WNT5A network with structural control inputs, consider solving Problem 1. Q,Qf, and Rs in Problem 1 are given as Q=[1000000], Qf=[10000000], and Rs=[00], respectively. The initial state is given as x0=[1010100]T. In addition, the control time N in Problem 1 is given by N=5. Finally, the constraint on duration of drug effectiveness is imposed for only u2s(k). The parameters Lu2s, Wu2s1, and Wu2s0 are given as Lu2s=1, Wu2s1=3, and Wu2s0=2, respectively. Thus, we can obtain the ILP problem (Problem B), where the dimension of binary variables is 130 and the number of inequalities is 264.We show the computational result. LetJ_* denote the optimal value of the lower bound of a given cost function in Problem 1. Let J¯* denote the upper bound of the cost function derived by using the optimal control input. First, consider the case of ρ=10-5. Then, we can obtain J_*=2 and J¯*=15, and us(k) is obtained as
(42)us(0)=⋯=us(4)=[10].Noting thatr2,1=0.1 and ρ=10-5(=0.15), all combinations of the Boolean functions are considered, and the value of ρ is not appropriate. In particular, J¯*=15 implies that  x1(k)=1, k=0,1,…,5, and is the trivial upper bound.Next, consider the case ofρ=0.2. Then, we can obtain J_*=J¯*=4, and us(k) is obtained as
(43)us(0)=us(1)=[01],us(2)=us(3)=us(4)=[10].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 ofJ¯* is 15, we see that in this case the effectiveness of control synthesis is clear.Finally, we discuss the computation time for solving Problem1. 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.


## 3.2. Artificial Example

In order to evaluate the computation time for solving Problem1, we consider one artificial example of a BN with 15 states and 3 control inputs. We stress that the existing method [11, 12, 17–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 2n×2n. In naive implementation using MATLAB [34], matrices with 215×215 cannot be created due to memory consumption, where we used the computer described previously.The optimal control problem is formulated. In this example, we consider 3 control inputs and 2 structural control inputs. Ifu1s(k)=1, then the system is given as follows:
(44)x(k+1)={f1(x(k),u(k))withtheprobability0.8,f2(x(k),u(k))withtheprobability0.2.Ifu2s(k)=1, then the system is given as follows:
(45)x(k+1)={f1(x(k),u(k))withtheprobability0.2,f2(x(k),u(k))withtheprobability0.8.The Boolean functionf1 is given by the following:
(46)x1(k+1)=x1(k)∧¬x6(k)∨u3(k),x2(k+1)=¬x4(k)∧u1(k)∨u3(k),x3(k+1)=x5(k)∧u1(k)∨¬x10(k)∧x12(k)∧u3(k),x4(k+1)=x2(k)∧x5(k)∧¬u1(k)∨¬x14(k),x5(k+1)=¬u1(k)∧x6(k)∧x7(k)∨x12(k)∧x14(k),x6(k+1)=x1(k)∧x6(k)∧x10(k)∨¬x15(k),x7(k+1)=x6(k)∧x7(k)∧x8(k)∨u2(k)∧¬u3(k),x8(k+1)=x5(k)∧¬u1(k)∨x10(k)∧u2(k)∧x13(k),x9(k+1)=x3(k)∧u1(k)∨¬x8(k)∧x11(k),x10(k+1)=x6(k),x11(k+1)=x6(k)∧x10(k)∨¬u2(k)∧u3(k),x12(k+1)=x12(k)∧¬x15(k),x13(k+1)=¬u3(k),x14(k+1)=¬x14(k)∧u3(k),x15(k+1)=x14(k)∧x15(k).The Boolean functionf2 is given by the following:
(47)x1(k+1)=x2(k)∧x4(k)∧¬x8(k),x2(k+1)=¬x2(k)∧x3(k)∨u3(k),x3(k+1)=x1(k)∨¬x2(k)∧x3(k)∧x4(k),x4(k+1)=¬x1(k)∧x2(k)∧u1(k)∨x14(k),x5(k+1)=¬u2(k)∧x13(k)∧x14(k)∧x15(k),x6(k+1)=¬x2(k)∨x5(k)∧u1(k)∧¬u3(k),x7(k+1)=u1(k)∧x13(k),x8(k+1)=x5(k)∧x13(k),x9(k+1)=¬x6(k),x10(k+1)=x2(k)∧u2(k)∧¬x12(k)∧u3(k),x11(k+1)=¬x5(k)∧u1(k)∨¬x15(k)∧u3(k),x12(k+1)=x7(k)∧x3(k),x13(k+1)=x7(k)∧u2(k)∧u3(k),x14(k+1)=x12(k)∨x14(k)∧u3(k),x15(k+1)=x8(k).From the previous setting,r1,1=0.8, r1,2=0.2, r2,1=0.2, r2,2=0.8 hold. In Problem 1, Q,Qf,R, and Rs are given as Q=[1⋯1], Qf=[10⋯10], R=[101], and  Rs=[00], respectively. The initial state and the parameter ρ are given as x0=[1⋯1]T and ρ=10-4, respectively. The constraint on duration of drug effectiveness is imposed for u1(k) and u2s(k). For u1, the parameters Lu1, Wu11, and Wu10 are given as Lu1=0, Wu11=3, and Wu10=3, respectively. For u2s, the parameters Lu2s, Wu2s1, and Wu2s0 are given as Lu2s=1, Wu2s1=3, and Wu2s0=2, respectively.Next, we discuss the computation time. Consider the two cases ofN=10 and N=20. Then, in the ILP problem (Problem B) obtained, the dimension of binary variables is 1420 for N=10 and 2840 for N=20, and the number of inequalities is 3381 for N=10 and 6731 for N=20. In the case of N=10, the computation time of the ILP problem was 96 [sec], where we used the computer described previously. In the case of N=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.


## 4. 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 Figure2, 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.



---

*Source: 246761-2013-08-21.xml*