Optimal Control of Asynchronous Boolean Networks Modeled by Petri Nets

Abstract

A Boolean network model is one of the models of gene regulatory networks, and is widely used in analysis and control. Although a Boolean network is a class of discrete-time nonlinear systems and expresses the synchronous behavior, it is important to consider the asynchronous behavior. In this paper, using a Petri net, a new modeling method of asynchronous Boolean networks with control inputs is proposed. Furthermore, the optimal control problem of Petri nets expressing asynchronous Boolean networks is formulated, and a solution method is proposed. The proposed approach provides us a new control method of gene regulatory networks.

Introduction

In recent years, there have been a lot of studies on modeling, analysis, and control of gene regulatory networks in both the control community and the theoretical biology community. 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, hybrid systems (piecewise affine models), and Petri nets have been developed so far (see e.g., [14]). In control problems, Boolean networks and hybrid systems are frequently used [1, 3, 4, 17, 18]. However, in the hybrid systems-based approach, a class of gene regulatory networks are limited to low-dimensional systems, because the computation time to solve the control problem is too long. In Boolean networks, dynamics such as interactions between genes are expressed by Boolean functions [15]. Although there is a criticism that a Boolean network is too simple as a model of gene regulatory networks, this model can be relatively applied to large-scale systems. Furthermore, since the behavior of gene regulatory networks is probabilistic by the effects of noise, a probabilistic Boolean network (PBN) has been proposed in [20].

Although a Boolean network is a class of discrete-time nonlinear systems and expresses the synchronous behavior, it is important to consider the asynchronous behavior. Asynchronous Boolean networks have been proposed in [12, 22]. In [12], we assume that the updating time of the concentration level of each gene is given in advance. In [22], asynchronous Boolean networks are modeled by non-deterministic dynamical systems. Furthermore, the asynchronous behavior is expressed as the probabilistic behavior. However, in these two methods, the asynchronous behavior is not directly modeled.

On the other hand, a Petri net is well known as a model expressing the asynchronous behavior [25]. A Petri net is a class of directed bipartite graphs, in which the nodes represent transitions and places. The methods to express asynchronous Boolean networks as Petri nets have been proposed in [8, 21]. However, in these methods, the control input is not considered, and these methods cannot be directly applied to the control problem. The control input in biological networks has the following significance. For example, the value of the control input expresses whether a stimulus is given to a cell. Then the control input is designed to obtain the state trajectory that transits from the initial state to the desired one. So the control input can represent the current status of therapeutic interventions, which are realized by radiation, chemotherapy, and so on. In order to develop gene therapy technologies (see e.g., [16, 19]) in future, it is important to consider control methods of Boolean networks.

Thus in this paper, the optimal control problem of asynchronous Boolean networks modeled by Petri nets is discussed. First, based on the method proposed in [8] and the notation of external input places [13, 23], we propose a new method to transform asynchronous Boolean networks with control inputs into Petri nets with external input places. Furthermore, the obtained Petri net is transformed into a logical dynamical system, which is a class of linear systems with binary states and binary inputs. Next, the optimal control problem is formulated, and the biological significance is also discussed by using a simple example. Finally, the solution method of the optimal control problem is proposed. In the proposed solution method, the optimal control problem is reduced to an integer linear programming (ILP) problem. The proposed approach provides us a new control method of gene regulatory networks.

This paper is organized as follows. In Section 2, synchronous Boolean networks and asynchronous Boolean networks are introduced. In Section 3, Petri nets expressing asynchronous Boolean networks are derived. In Section 4, we explain the method to transform the obtained Petri net into a logical dynamical system. In Section 5, the optimal control problem is formulated, and in Section 6, a solution method is proposed. In Section 7, we conclude this paper. Notation: Let R denote the set of real numbers. Let {0, 1}m×n denote the set of m × n matrices, which consists of elements 0 and 1. For the finite set M , let |M | denote the number of elements. Let In and 0m×n denote the n × n identity matrix and the m × n zero matrix, respectively. For simplicity of notation, we sometimes use the symbol 0 instead of 0m×n, and the symbol I instead of In. For a matrix M , let M T denote the transpose of M . 2

Synchronous/Asynchronous Boolean Networks First, we explain synchronous Boolean networks (SBNs). where x ∈ {0, 1}n is the state (e.g., the concentration of genes), k = 0, 1, 2, . . . is the discrete time. fa : {0, 1}n → {0, 1}n is a given Boolean function with logical operators such as AND (∧), OR (∨), and NOT (¬). Since the SBN ( 1 ) is deterministic, x(k + 1) is uniquely determined for a given x(k).

To consider the control problems, we add the control input to the SBN ( 1 ) as follows:

x(k + 1) = f (x(k), u(k)) where u ∈ {0, 1}m is the control input, i.e., 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 i-th element of the state x, the i-th element of the control input u and the i-th element of the Boolean function f are denoted by xi, ui and fi, respectively. Also in the SBN ( 2 ), x(k + 1) is uniquely determined for given x(k) and u(k).

Next, we explain asynchronous Boolean networks (ABNs). Some methods for expressing ABNs have proposed so far [12, 22]. In this section, we explain a method in [22] for expressing ABNs as nondeterministic systems. A Petri net-based approach will be explained in Section 3. In the case that ABNs are expressed as nondeterministic systems, the behavior of ABNs is obtained by the union of the behaviors of the following n SBNs: Σi : xi(k + 1) = fi(x(k), u(k)), xj (k + 1) = xj (k), ∀j ∈ {1, 2, . . . , n} \ {i}, where i = 1, 2, . . . , n.

We show an example of SBNs and ABNs.

Example 1. As a simple example, consider the following SBN of an apoptosis network [9]: ⎧ 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 (IAP) is denoted by x1, the concentration level of the active caspase 3 (C3a) by x2, and the concentration level of the active caspase 8 (C8a) by x3. The concentration level of the tumor necrosis factor (TNF, a stimulus) is denoted by u, and is regarded as the control input.

In the case of synchronous Boolean dynamics, state transitions can be computed by directly using ( 4 ). For example, for x(0) = [ 1 1 1 ]T and u(k) = 0, we obtain x( 1 ) = [ 0 0 1 ]T . By computing the transition from each state, we obtain the state transition diagram in Fig. 1 (left). In Fig. 1 (left), the number assigned to each node denotes x1, x2, x3 (elements of the state), Consider the following SBN: x(k + 1) = fa(x(k)) ( 1 ) ( 2 ) ( 3 ) ( 4 ) State transitions can be computed by using ( 5 ), ( 6 ), ( 7 ). For example, for x(0) = [ 1 1 1 ]T and u(k) = 0, we obtain x( 1 ) = {[ 0 1 1 ]T , [ 1 0 1 ]T , [ 1 1 1 ]T }. In a similar way, by computing the transition from each state, we obtain the state transition diagram in Fig. 1 (right).

Comparing the left figure with the right figure in Fig. 1, we see that a part of behaviors is clearly different.

In this paper, ABNs are modeled by Petri nets, not multiple SBNs. By using Petri nets, we can consider several situations. For example, although each xi is independently activated in ( 3 ), activation of combinations of xi can be considered. 3

Transformation of Boolean Networks into Petri Nets Using each Boolean function fi in the SBN ( 2 ), consider to express an ABN as a Petri net. Based on a complementary-place transformation, a Petri net expressing an ABN has been proposed in [8], but the control input has not been considered. In this paper, as an extension of the method in [8], we propose a modeling method of a Petri net expressing an ABN with the control input.

First, some notations are prepared. By I(j), j = 1, 2, . . . , n, denote the state and the control input included in the Boolean function fj. In the example of ( 4 ), we obtain I( 1 ) = {x2, u}, I( 2 ) = {x1, x3}, and I( 3 ) = {x2, u}. Next, we define a logical parameter Kj(X) ∈ {0, 1}, X ⊆ I(j), j = 1, 2, . . . , n. If the value of each element included in X is ‘1’, and the value of the other element is ‘0’, then either Kj(X) = 1 or Kj(X) = 0 is determined. In the example of ( 4 ), we obtain K1(∅) = 0, K1({x2}) = 0, K1({u}) = 1, K1({x2, u}) = 0, K2(∅) = 0, K2({x1}) = 0, K2({x3}) = 1, K2({x1, x3}) = 0,

K3(∅) = 0, K3({x2}) = 1, K3({u}) = 1, K3({x2, u}) = 1.

Next, consider to derive a Petri net expressing Boolean networks. In the derived Petri net, the number of places is given as 2(n + m), that is, for each xi in ( 2 ), two places xi and xi are prepared. In a similar way, for each ui in ( 2 ), two places ui and ui are prepared. xi and ui are called complementary places [8]. The number of transitions is given as n

2 + 22 i+=122|I(i)|. In the example of ( 4 ), the number of transitions is given as 2 2 = 12. From the property of Boolean networks, the following assumptions are made.

Assumption 1 The maximum number of tokens in each place is equal to 1. Assumption 2 A sum of the number of tokens in xi (ui) and that in xi (ui) is equal to 1.

In addition, suppose that ui and ui are given as an external input place [13, 23]. In ui and ui, a token is arbitrary generated, but the above two assumptions must be satisfied.

Under the above preparations, we define a Petri net expressing an ABN. In [8], the Petri net expressing an ABN without the control input is defined. The following definition gives the Petri net expressing an ABN with the control input, and is an extension of the definition in [8].

Definition 1. For a given SBN ( 2 ), the Petri net expressing an ABN is defined as follows:

Nc = (P ∪ Pc, T, P re, P ost) ( 8 ) where • P = {x1, x1, x2, x2, . . . , xn, xn} is the set of places, • Pc = {u1, u1, u2, u2, . . . , um, um} is the set of external input places, • T = {txi,X , i = 1, 2, . . . , n, X ⊆ I(i)} is the set of transitions, • P re : (P ∪ Pc) × T → {0, 1} is the mapping defining arcs between places and transitions, • P ost : T × P → {0, 1} is the mapping defining arcs between transitions and transitions.

The functions P re and P ost are defined as follows:

(i) Case of xi ∈ I(i) (xi is not a self-regulator): For a given transition txi,X , the following terms are defined (all the other terms are equal to zero): P re(xi, txi,X ) = P ost(txi,X , xi) = 1 − Ki(X), P re(xi, txi,X ) = P ost(txi,X , xi) = 1 − Ki(X), P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ X, P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ I(i) − X, P re(uj , txi,X ) = 1, ∀uj ∈ X,

P re(uj , txi,X ) = 1, ∀uj ∈ I(i) − X. (ii) Case of xi ∈ I(i) (xi is a self-regulator): Consider a given transition txi,X . if xi ∈ X, then only the case of Ki(X) = 0 is considered. Therefore, the following terms are defined:

P re(xi, txi,X ) = P ost(txi,X , xj ) = 1, P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ X, xj = xi, P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ I(i) − X, P re(uj , txi,X ) = 1, ∀uj ∈ X,

P re(uj , txi,X ) = 1, ∀uj ∈ I(i) − X. if xi ∈ X, then only the case of Ki(X) = 1 is considered. Therefore, the following terms are defined:

P re(xi, txi,X ) = P ost(txi,X , xj ) = 1, P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ X, P re(xj , txi,X ) = P ost(txi,X , xj ) = 1, ∀xj ∈ I(i) − X, xj = xi, P re(uj , txi,X ) = 1, ∀uj ∈ I(i) − X,

P re(uj , txi,X ) = 1, ∀uj ∈ X.

In the above definition, a sum of the number of tokens in ui and that in u¯i becomes zero by firing some transition. In this case, to satisfy Assumption 1 and Assumption 2, a token is generated in either ui or u¯i.

We show a simple example.

Example 2. Consider the following simple SBN: x1(k + 1) = x2(k), x2(k + 1) = u(k) ( 9 ) From x1(k + 1) = x2(k), we obtain K1(∅) = 0 and K1({x2}) = 1. In a similar way, from x1(k + 1) = u(k), we obtain K2(∅) = 0 and K1({u}) = 1. Then we consider four transitions tx1,∅, tx1,{x2}, tx2,∅, and tx2,{u}. We denote these transitions by tx1 , tx1,x2 , tx2 , and tx2,u, respectively. Then we obtain the Petri net in Fig. 2. In addition, P re and P ost in ( 8 ) are obtained as follows: Suppose that one token is included in place x¯1, x¯2, and u. Then the transition tx2,u may fire. If the transition tx2,u fire, then one token is moved from x¯2 and u to x2. A pair of x2 and x¯2 satisfies Assumption 1 and Assumption 2, but a pair of u and u¯ does not satisfy Assumption 2. So one token must be added in either u and u¯ with fire. 4

Transformation of Petri Nets into Logical Dynamical Systems To consider the optimal control problem, it is desirable to transform a Petri net ( 8 ) into some linear form. In this section, based on a framework on modeling of hybrid dynamical systems [5], logical dynamical systems expressing Petri nets are derived.

Logical dynamical systems are given as a pair of linear state equations and linear inequality constraints with binary state variables and binary control input variables. A general form of logical dynamical systems is defined as follows. Definition 2. A logical dynamical system is given as x(k + 1) = Ax(k) + Bv(k), Cx(k) + Dv(k) ≤ E ( 10 ) ( 11 ) where x(k) ∈ {0, 1}nd is the state variable, and v(k) ∈ {0, 1}md is the input variable including auxiliary variables,

In modeling of hybrid dynamical systems, a mixed logical dynamical (MLD) system [5] is well known. The MLD system can be derived by replacing x(k) ∈ {0, 1}nd and v(k) ∈ {0, 1}md in ( 10 ), ( 11 ) with x(k) ∈ Rnc × {0, 1}nd and v(k) ∈ Rmc × {0, 1}md . Since a Petri net is a class of discrete event systems, the MLD system is not used, and the logical dynamical system ( 10 ), ( 11 ) is used.

Let us consider to transform the Petri net ( 8 ) into the logical dynamical system ( 10 ), ( 11 ). Some notations are prepared. By xi(k), x¯i(k), ui(k), u¯i(k) ∈ {0, 1}, denote existence or non-existence of a token in place xi, x¯i, ui, u¯i ∈ {0, 1} at time k. k may be regarded as the k-th firing in a firing sequence. Since xi(k), x¯i(k), ui(k), u¯i(k) are binary variables, Assumption 1 satisfies. To satisfy Assumption 2, xi(k) + x¯i(k) = 1 and ui(k) + u¯i(k) = 1 are imposed. Next, by txi,Xj (k) ∈ {0, 1}, i = 1, 2, . . . , n, j = 1, 2, . . . , 2|I(i)|, denote fire in the transition txi,Xj . If txi,Xj (k) = 1, then the transition txi,Xj fires at time k. Otherwise, txi,Xj does not fire. By using these notations, Petri nets ( 8 ) is transformed into logical dynamical systems.

First, a simple example is shown.

Example 3. Consider the Petri net in Fig. 2. From the property of fire, we obtain the following system expressing the Petri net in Fig. 2: x1(k + 1) = tx1,x2 (k)x¯1(k)x2(k) + (1 − tx1,x2 (k))x1(k) −tx1 (k)x1(k)x¯2(k), −tx1,x2 (k)x¯1(k)x2(k), x¯1(k + 1) = tx1 (k)x1(k)x¯2(k) + (1 − tx1 (k))x¯1(k) x2(k + 1) = tx1,x2 (k)x¯1(k)x2(k) + tx2,u(k)x¯2(k)u(k)

+(1 − tx1,x2 (k) − tx2,u(k))x2(k) − tx2 (k)x2(k)u¯(k), x¯2(k + 1) = tx1 (k)x1(k)x¯2(k) + tx2 (k)x2(k)u¯(k) +(1 − tx1 (k) − tx2 (k))x¯2(k) −tx1 (k)x1(k)x¯2(k) − tx2,u(k)x¯2(k)u(k).

If the number of firing transitions at each time is limited to 1, then the inequality condition

tx1 (k) + tx1,x2 (k) + tx2 (k) + tx2,u(k) ≤ 1 is imposed. Then x1(k), x¯1(k), x2(k), x¯2(k) ∈ {0, 1}, x1(k) + x¯1(k) = 1, and x2(k) + x¯2(k) = 1 hold thanks to the condition ( 16 ), u(k), u¯(k) ∈ {0, 1}, u(k) + u¯(k) = 1, and the initial condition x1(0), x¯1(0), x2(0), x¯2(0) ∈ {0, 1}, x1(0) + x¯1(0) = 1, x2(0) + x¯2(0) = 1. The system ( 12 )–( 15 ) is a nonlinear system, and can be linearized by using Lemma 1 in Appendix A. For example, z1(k) = tx1,x2 (k)x¯1(k)x2(k) is equivalent to ( 12 ) ( 13 ) ( 14 ) ( 15 ) ( 16 ) tx1,x2 (k) + x¯1(k) + x2(k) − z1(k) ≤ 2, −tx1,x2 (k) − x¯1(k) − x2(k) + 3z1(k) ≤ 0.

By applying Lemma 1 to the other terms, we can obtain the logical dynamical system expressing the Petri net in Fig. 2.

From this example, we see that the Petri net ( 8 ) can be expressed by the logical dynamical system ( 10 ), ( 11 ). Furthermore, when the Petri net ( 8 ) is expressed by the logical dynamical system ( 10 ), ( 11 ), variables x, v are given as x(k) = [ x1(k) x¯1(k) x2(k) x¯2(k) · · · xn(k) x¯n(k) ]T , v(k) = [ U (k) T (k) Z(k) ]T , U (k) = [ u1(k) u¯1(k) u2(k) u¯2(k) · · · um(k) u¯m(k) ]T ,

T (k) = tx1,X1 (k) · · · tx1,X2|I( 1 )| (k) · · · txn,X1 (k) · · · txn,X2|I(n)| (k) where Z(k) is a auxiliary binary variable obtained by applying Lemma 1. Matrices/vectors A, B, C, D, E in ( 10 ), ( 11 ) can be derived from P re, P ost in the Petri net ( 8 ) and Lemma 1.

Remark 1. One of the simple methods for modeling of ABNs is to express ABNs as switched systems with 2n subsystems. 2n subsystems are derived by all combinations of Boolean functions f1, f2, . . . , fn in ( 2 ). In the case of ( 9 ) in Example 2, the following 2n = 4 subsystems: T Σ1 : x1(k + 1) = x1(k), x2(k + 1) = x2(k), Σ2 : x1(k + 1) = x2(k), x2(k + 1) = x2(k), Σ3 : x1(k + 1) = x1(k), x2(k + 1) = u(k), Σ4 : x1(k + 1) = x2(k), x2(k + 1) = u(k) are obtained. From these subsystems, we can obtain the following system expressing an ABN: x1(k + 1) = (δ1(k) + δ3(k))x1(k) + (1 − δ1(k) − δ3(k))x2(k), x2(k + 1) = (δ1(k) + δ2(k))x2(k) + (1 − δ1(k) − δ2(k))u(k) ( 17 ) ( 18 ) where δ1(k), δ2(k), δ3(k) are binary variables satisfying δ1(k) + δ2(k) + δ3(k) ≤ 1, and correspond to Σ1, Σ2, Σ3, respectively. Furthermore, z1 := δ1x1, z2 := δ3x1, z3 := δ1x2, z4 := δ3x2, z5 := δ2x2, z6 := δ1u, and z7 := δ2u are defined, and Lemma 1 is applied to z1, z2, . . . , z7. Thus we can obtain the logical dynamical system ( 10 ), ( 11 ). This method is called here a direct approach.

Comparing ( 17 ), ( 18 ) with ( 12 )–( 15 ), we see that the system ( 17 ), ( 18 ) is simpler than the system ( 12 )–( 15 ). However, for general cases, this fact does not hold. Here, we focus on the dimension of binary variables to switch Boolean functions. In ( 12 )–( 15 ), this dimension corresponds to the dimension of binary variables assigned to transitions, and is given as 4. In ( 17 ), ( 18 ), this dimension is given as 3. In general, in the direct approach, the dimension of binary variables to switch Boolean functions is given as 2n − 1. In the proposed Petri net-based approach, this dimension, i.e., |T (k)| is given as in=1 2|I(i)|. In real gene regulatory networks, it is well known that |I(i)| is relatively smaller than n (see e.g., [2]). For example, in the case of n = 10 and |I(i)| = 3, 2n − 1 = 1023 and in=1 2|I(i)| = 80 are obtained. Thus, in modeling of real gene regulatory networks, it is not appropriate to use the direct approach, and the proposed Petri net-based method provides us a simpler modeling method of ABNs.

N−1 For the logical dynamical system ( 10 ), ( 11 ) expressing the Petri net ( 8 ), consider the following optimal control problem.

Problem 1. For the logical dynamical system ( 10 ), ( 11 ) expressing the Petri net ( 8 ), suppose that the initial state x(0) = x0 satisfying Assumption 1 is given. Then find an input sequence v(0), v( 1 ), . . . , v(N − 1) minimizing the linear cost function

J = {Qx(k) + Rv(k)} + Qf x(N ) ( 19 ) k=0 where Q, Qf ∈ R1×nd , R ∈ R1×md are weighting vectors whose element is a non-negative real number.

For simplicity of discussion, a linear function with respect to x and u is considered as a cost function, but a quadratic cost function may be used. In addition, using the offset vector xd ∈ {0, 1}nd and vd ∈ {0, 1}md , x(k) and v(k) may be replaced to xˆ(k) := x(k) − xd and vˆ(k) := v(k) − vd. Then it is necessary that the cost function ( 19 ) is also replaced to J = iN=−01 {Q|xˆ(i)| + R|vˆ(i)|} + Qf |xˆ(N )|. In addition, N must be determined according to a given biological network. Although a longer N is desirable, the computation time to solve Problem 1 must be also considered. For a small N , Problem 1 may be repeatedly solved at each time. This policy is well known as model predictive control [6].

In Problem 1, we assume that an input sequence v(0), v( 1 ), . . . , v(N − 1) is arbitrarily determined. However, there is a possibility that a given biological system does not satisfy this assumption. Then suppose that some candidates of input sequences are given. In Problem 1, the optimal input sequence minimizing the cost function ( 19 ) is selected among the set of the candidates B ⊆ {0, 1}mdN . This extension is easy. In this sense, Problem 1 can be applied to optimal control of asynchronous Boolean networks such that the updating time of each state is given in advance [12]. Of course, this problem can also be applied to optimal control of SBNs. Thus Problem 1 includes several situations.

Next, we show an example for setting weighting vectors from the biological viewpoint.

Example 4. Consider the Boolean network expressing an apoptosis network in Example 1 again. From ( 4 ), we obtain the Petri net ( 8 ) with 6 places, 2 external input places, and 12 transitions. In addition, from the obtained Petri net, we obtain the logical dynamical system ( 10 ), ( 11 ). For the obtained logical dynamical system, we consider to find a control strategy such that a stimulus is not applied as much as possible, and cell survival is achieved. u(k) = 0 implies that a stimulus is not applied to the system, and x1(k) = 1, x2(k) = 0 express cell survival [9]. Then as one of appropriate cost functions, we can consider the following cost function

J =

N−1 k=0 {10|x1(k) − 1| + 10|x2(k) − 0| + u(k)}

By the appropriate coordinate transformation, this cost function can be rewritten as the form of ( 19 ). See also [10–12] for biological examples on the optimal control problems. 6

Reduction to an Integer Linear Programming Problem Finally, let us consider to reduce Problem 1 to an integer linear programming (ILP) problem.

Problem 1 can be rewritten by using ( 10 ), ( 11 ). First, by using x(k) = Akx0 +

Ai−1Bv(k − i) k i=1 obtained from the state equation ( 10 ), we obtain x¯ = A¯x0 + B¯v¯ ( 20 ) where x¯ := [ xT (0) xT ( 1 ) · · · xT (N ) ]T , v¯ := [ vT (0) vT ( 1 ) · · · vT (N − 1) ]T and Furthermore, from the linear inequality ( 11 ), we obtain where

⎢ A¯ := ⎢⎢ ⎢⎢ ⎣ ⎡ Ind ⎤

AA2 ⎥⎥⎥ , B¯ := ⎢⎢⎢⎢ ... ⎥⎥⎦ ⎢⎢

⎣ AN ⎡ 0nd×md 0 · · · 0 ⎤

ABB .0. . ·. ·. .· 0... ⎥⎥⎥⎥ .

... . . . . . . 0 ⎥⎦⎥

AN−1B · · · AB B C¯x¯ + D¯ v¯ ≤ E

¯ C¯ := [ block-diag(C, C, . . . , C) 0 ] , D¯ := block-diag(D, D, . . . , D), E¯ :=

ET

ET

· · · ET T .

J = Q¯x¯ + R¯v¯ Next, the cost function ( 19 ) can also be rewritten as ( 21 ) ( 22 ) where Q¯ := [ Q · · · Q Qf ] and R¯ := [ R · · · R ]. By substituting ( 20 ) into ( 21 ) and ( 22 ), we obtain the following theorem. Theorem 1. Problem 1 is equivalent to the following ILP problem.

Problem A:

subject to find v¯ ∈ {0, 1}mdN , min (R¯ + Q¯B¯)v¯ + Q¯A¯x0,

C¯B¯ + C¯ E¯ − C¯A¯x0 .

−L¯W¯ v¯ ≤ −ln ρ

Problem A can be solved by using a suitable ILP solver such as IBM ILOG CPLEX Optimizer [24]. 7

Conclusion

In this paper, we have discussed optimal control of asynchronous Boolean networks with control inputs. First, we have proposed a method to transform Boolean networks with control inputs into Petri net with external input places. Next, after the obtained Petri net is transformed into a logical dynamical system, the optimal control problem has been formulated. This problem is a general formulation including several biological situations. Finally, the optimal control problem has been reduced to an integer linear programming problem. The proposed approach will be effective for control of several biological systems modeled by Boolean networks.

One of the future works is to apply the proposed approach to biological Boolean networks. From the practical viewpoint, an extension to probabilistic Boolean networks is also important. In addition, for large-scale Boolean networks, the computation time to solve the problem will be long. So it is significant to consider to reduce the computation time to solve the problem. Then one of methods to overcome this difficulty is to use a SAT (satisfiability problem) solver such as zChaff [26]. It is also important to consider approximate solution methods.

This work was supported by Grant-in-Aid for Young Scientists (B) 23760387 and Scientific Research (C) 21500009.

A

Linearization of The Product of Binary Variables The product of binary variables can be linearized by using the following lemma [7].

Lemma 1. Suppose that binary variables δj ∈ {0, 1}, j ∈ J are given, where J is some index set. Then z = j∈J δj is equivalent to the following linear inequalities δj − z ≤ |J | − 1, −

δj + |J |z ≤ 0 j∈J j∈J where |J | is the cardinality of J .

References

1. T. Akutsu , M. Hayashida , W.-K. Ching, and M. K. Ng , Control of Boolean networks: Hardness results and algorithms for tree structured networks , Journal of Theoretical Biology , vol. 244 , pp. 670 - 679 , 2007 .

2. M. Andrecut , S. A. Kauffman , and A. M. Madni , Evidence of Scale-Free Topology in Gene Regulatory Network of Human Tissues , Int'l Journal of Modern Physics C , vol. 19 , no. 2 , pp. 283 - 290 , 2008 .

3. S. Azuma , E. Yanagisawa, and J. Imura , Controllability Analysis of Biosystems Based on Piecewise Affine Systems Approach, IEEE Trans. on Automatic Control , vol. 53 , no. 1 , pp. 139 - 152 , 2008 .

4. C. Belta , J. Schug , T. Dang , V. Kumar , G. J. Pappas , H. Rubin , and P. Dunlap , Stability and reachability analysis of a hybrid model of luminescence in the marine bacterium Vibrio fischeri , Proc. 40th IEEE Conf. on Decision and Control , pp. 869 - 874 , 2001 .

5. A. Bemporad and M. Morari , Control of systems integrating logic, dynamics, and constraints . Automatica , Vol. 35 , pp. 407 - 427 , 1999 .

6. E. F. Camacho , C. Bordons , Model Predictive Control, 2nd Edition , Springer, 2004 .

7. T. M. Cavalier , P. M. Pardalos , and A. L. Soyster , Modeling and integer programming techniques applied to propositional calculus , Computer & Operations Research , vol. 17 , no. 6 , pp. 561 - 570 , 1990 .

8. C. Chaouiya , E. Remy, P. Ruet , and D. Thieffry , Qualitative Modelling of Genetic Networks: from Logical Regulatory Graphs to Standard Petri Nets, 25th Int'l Conf. on Application and Theory of Petri Nets, LNCS 3099 , pp. 137 - 156 , 2004 .

9. M. Chaves , Methods for Qualitative Analysis of Genetic Networks , Proc. European Control Conference 2009 , pp. 671 - 676 , 2009 .

10. A. Datta , A. Choudhary , M. L. Bittner , and E. R. Dougherty , External Control in Markovian Genetic Regulatory Networks, Machine Learning , vol. 52 , pp. 169 - 191 , 2003 .

11. A. Datta , A. Choudhary , M. L. Bittner , and E. R. Dougherty , External control in Markovian genetic regulatory networks: the imperfect information case , Bioinfomatics , vol. 20 , pp. 924 - 930 , 2004 .

12. B. Faryabi , J.-F. Chamberland , G. Vahedi , A. Datta , and E. R. Dougherty , Optimal Intervention in Asynchronous Genetic Regulatory Networks , IEEE Journal of Selected Topics in Signal Processing , Vol. 2 , No. 3 , pp. 412 - 423 , 2008 .

13. A. Ichikawa and K. Hiraishi , Analysis and Control of Discrete Event Systems Represented by Petri Nets, Discrete Event Systems: Models and Applications , LNCIS 103, pp. 115 - 134 , 1988 .

14. H. D. Jong, Modeling and Simulation of Genetic Regulatory Systems: A Literature Review , Journal of Computational Biology , vol. 9 , pp. 67 - 103 , 2002 .

15. S. A. Kauffman , Metabolic stability and epigenesis in randomly constructed genetic nets , Journal of Theoretical Biology , vol. 22 , pp. 437 - 467 , 1969 .

16. M.-S. Kim et al., DNA demethylation in hormone-induced transcriptional derepression , Nature , vol. 461 , pp. 1007 - 1012 , 2009 .

17. K. Kobayashi , J. Imura , and K. Hiraishi , Polynomial-Time Algorithm for Controllability Test of a Class of Boolean Biological Networks , EURASIP Journal on Bioinformatics and Systems Biology , vol. 2010 , Article ID 210685, 12 pages, 2010 .

18. C. J. Langmead and S. K. Jha , Symbolic Approaches to Finding Control Strategies in Boolean Networks , Journal of Bioinformatics and Computational Biology , vol. 7 , no. 2 , pp. 323 - 338 , 2009 .

19. H. Santos-Rosa and C. Caldas , Chromatin modifier enzymes, the histone code and cancer , European Journal of Cancer , vol. 41 , pp. 2381 - 2402 , 2005 .

20. I. Shmulevich , E. R. Dougherty , S. Kim , and W. Zhang , Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks , Bioinformatics , vol. 18 , pp. 261 - 274 , 2002 .

21. L.J. Steggles , R. Banks , and A. Wipat , Modelling and Analysing Genetic Networks: From Boolean Networks to Petri Nets , Int'l Conf. on Computational Methods in Systems Biology, LNBI 4210 , pp. 127 - 141 , 2006 .

22. L. Tournier and M. Chaves , Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics , Journal of Theoretical Biology , vol. 260 , pp. 196 - 209 , 2009 .

23. T. Ushio, Maximally Permissive Feedback and Modular Control Synthesis in Petri Nets with External Input Places , IEEE Trans on Automatic Control , Vol. 35 , No. 7 , pp. 844 - 848 , 1990 .

24. http://www-01.ibm.com/software/integration/optimization/cplex-optimizer/

25. http://www.informatik.uni-hamburg.de/TGI/PetriNets/

26. http://www.princeton.edu/˜chaff/zchaff.html