# Verification and Optimal Control of Context-Sensitive Probabilistic Boolean Networks Using Model Checking and Polynomial Optimization

Authors: Koichi Kobayashi; Kunihiko Hiraishi
Journal: The Scientific World Journal
 (2014)

---

## Abstract

One of the significant topics in systems biology is to develop control theory of gene regulatory networks (GRNs). In typical
control of GRNs, expression of some genes is inhibited (activated) by manipulating external stimuli and expression of other genes. It is expected to apply control theory of GRNs 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 GRNs, and gene expression is expressed by a binary value (ON or OFF). In particular, a context-sensitive probabilistic Boolean network (CS-PBN), which is one of the extended models of BNs, is used. For CS-PBNs, the verification problem and the optimal control problem are considered. For the verification problem, a solution method using the probabilistic model checker PRISM is proposed. For the optimal control problem, a solution method using polynomial optimization is proposed. Finally, a numerical example on the WNT5A network, which is related to melanoma, is presented. The proposed methods provide us useful tools in control theory of GRNs.


---

## Body

## 1. Introduction

Control of gene regulatory networks (GRNs) is one of the significant topics in the field of systems biology, and is also one of the basics of therapeutic interventions (see, e.g., [1]) in the future. Furthermore, in recent years, the experimental result on control of GRNs 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 GRNs can be realized. Motivated by the above backgrounds, we study a control method of GRNs.GRNs are in general modeled 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–10]. In the hybrid systems-based approach, the classes of GRNs are limited to low-dimensional systems, because the computation time to solve the control problem is too long. In a BN, gene expression is expressed by a binary value (ON or OFF), and dynamics such as interactions between genes are expressed by Boolean functions [3]. In, for example, [11], it is pointed out that a BN is too simple as a model of GRNs. However, there is an advantage that a BN can be relatively applied to large-scale systems. Furthermore, since the behavior of GRNs is stochastic by the effects of noise, it is appropriate that a Boolean function is randomly decided at each time among the candidates of Boolean functions. Thus, a probabilistic Boolean network (PBN) has been proposed in [12], and further, a context-sensitive PBN (CS-PBN) has been proposed as a general form of PBNs [13, 14]. In a CS-PBN, the deciding time is also randomly selected.In this paper, a CS-PBN is adopted as a model of GRNs, and for CS-PBNs, the verification problem and the optimal control problem are considered. In the verification of PBNs, a solution method using the probabilistic model checker PRISM [15] has been proposed in [16]. However, this PRISM-based method for PBNs has not been extended to that for CS-PBNs. In optimal control of PBNs and CS-PBNs, many results have been obtained so far (see, e.g., [13, 14, 17–24]). In many existing results, state transition diagrams with 2n nodes (i.e., 2n×2n transition probability matrices) must be computed for a (CS-)PBN with n genes. 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, 23] control methods in which state transition diagrams are not computed. Comparing the methods in [22, 23] with other existing results [13, 14, 17–21, 24], the methods in [22, 23] can relatively handle more large-scale GRNs. The method in [22] can be applied to PBNs and CS-PBNs, but the expected value of a given nonnegative function cannot be evaluated as a cost function (objective function). In the method in [23], the expected value of a given nonnegative function can be used as a cost function, and the optimal control problem is reduced to a polynomial optimization problem. However, this method has been proposed for PBNs, and an extension to CS-PBNs has not been discussed so far. Thus, for verification of CS-PBNs, the PRISM-based method for PBNs [16] is extended to that for CS-PBNs. For optimal control of CS-PBNs, a solution method using polynomial optimization [23] is extended to that for CS-PBNs. Furthermore, the effectiveness of the proposed methods is presented by a numerical example on the WNT5A network, which is related to melanoma. The proposed methods provide us useful tools in control theory of GRNs.This paper is organized as follows. In Section2.1, a CS-PBN is explained. In Section 2.2, a solution method for the verification problem is proposed. In Section 2.3, a solution method for the optimal problem is proposed. In Section 3, a numerical example is presented. In Section 4, we conclude this paper.Notation. 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. For a matrix M, M⊤ denotes the transpose of M.


## 2. Materials and Methods

### 2.1. Context-Sensitive Probabilistic Boolean Networks

First, we introduce a probabilistic Boolean network (PBN). Consider the following PBN:(1)x1(k+1)=f(1)(k,x(k),u(k)),x2(k+1)=f(2)(k,x(k),u(k)),⋮xn(k+1)=f(n)(k,x(k),u(k)),
where x=[x1x2⋯xn]⊤∈{0,1}n is the state (e.g., the expression of genes), u=[u1u2⋯um]⊤∈{0,1}m is the control input (e.g., the expression of genes), that is, the value of u can be arbitrarily given, and k=0,1,2,… is the discrete time. For a fixed k∈{0,1,…}, f(i):{0,1,…}×{0,1}n×{0,1}m→{0,1}1 is a given Boolean function consisting of logical operators such as AND (∧), OR (∨), and NOT (¬). In deterministic Boolean networks, x(k+1) is uniquely determined for given k, x(k), and u(k). In PBNs, the candidates of f(i)(k,x(k),u(k)) are given, and for each xi, selecting one Boolean function is probabilistically independent at each time. Let fj(i)(x(k),u(k)), j=1,2,…,l(i), denote the candidates of f(i)(k,x(k),u(k)). The probability that fj(i)(x(k),u(k)) is selected is defined by
(2)cj(i):=Prob(f(i)(k,x(k),u(k))=fj(i)(x(k),u(k))).
Then, the following relation:
(3)∑j=1l(i)‍cj(i)=1
must be satisfied. 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., [25]).Example 1.
As a simple example, consider the following deterministic Boolean network of an apoptosis network [26, 27] (see also Figure 1):
(4)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. Although Boolean dynamics in the above system are synchronous, both synchronous and asynchronous dynamics will be included. From this viewpoint, we consider the following PBN induced by the above system:
(5)f(1)={f1(1)=¬x2(k)∧u(k),c1(1)=0.6,f2(1)=x1(k),c2(1)=0.4,(6)f(2)={f1(1)=¬x1(k)∧x3(k),c1(2)=0.7,f2(1)=x2(k),c2(2)=0.3,(7)f(3)={f1(3)=x2(k)∨u(k),c1(3)=0.8,f2(3)=x3(k),c2(3)=0.2,
where l(1)=l(2)=l(3)=2, and we give cj(i) satisfying ∑j=1l(i)‍cj(i)=1. In addition, all state trajectories can be expressed as the state transition diagram with 23 nodes.Figure 1
Simplified model of an apoptosis network. Activation (solid), inhibition (broken).In PBNs, we suppose that selecting one Boolean function is probabilistically independent at each time. However, it will be natural to consider the situation that switches of Boolean functions do not occur frequently. From this viewpoint, a context-sensitive PBN (CS-PBN) has been proposed in [13, 14]. In CS-PBNs, the deciding time of Boolean functions is also selected randomly. Hereafter, let q∈[0,1] denote the probability that Boolean functions are switched at time k, and a pair of the system (1) and q is called a CS-PBN.To compare CS-PBNs with PBNs, consider (7) as a simple example. In PBNs, a switch of f1(3) and f2(3) functions does not depend on the Boolean function at time k-1. In CS-PBNs, a switch of f1(3) and f2(3) is decided by the discrete-time Markov chain in Figure 2. In other words, this switch depends on the Boolean function at time k-1. Owing to this difference, a control/verification method for CS-PBNs cannot be directly derived from that for PBNs.Figure 2
Discrete-time Markov chain inf(3).


### 2.2. Verification Using Model Checking

First, the reachability problem is formulated as the verification problem studied in this paper. The reachability problem is one of the typical verification problems. For a given CS-PBN, the outputy(k)=[y1(k)y2(k)⋯yp(k)]⊤∈{0,1}p is defined, where yi=xj, i=1,2,…,p, j∈𝒥⊆{1,2,…,n}. We remark that the output does not mean the measured signal. First, the reachability problem is formulated as follows.Problem 2 (reachability problem).
Suppose that, for CS-PBN with the output, the initial statex(0)=x0, the initial Boolean function f(i)(0,x(0),u(0))=fj0(i)(i)(x(0),u(0)) (j0(i)∈{1,2,…,l(i)}), the control time N, and the target output yf are given (u(0) is not given). Then, find a maximum probability Pmax⁡ that y(k)=yf holds within time N by manipulating a control input sequence u(0),u(1),…,u(N-1).In the standard reachability problem, only terminal time is focused, and it is checked whethery(N)=yf holds or not. In this paper, we focus on not only terminal time N but also other times 0,1,…,N-1. Furthermore, since a CS-PBN has the control input, which can be regarded as a nondeterministic variable, we find a maximum probability satisfying the condition.Next, we will propose a solution method for Problem2. As a preparation, the following lemma [28] is introduced.Lemma 3.
Consider two binary variablesδ1,δ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. By using this lemma, a Boolean function can be transformed into a polynomial on the real number field.To solve Problem2, the probabilistic model checker PRISM [15] is used. PRISM supports a discrete-time Markov chain (DT-MC), a continuous-time Markov chain (CT-MC), and a Markov decision process (MDP). PRISM consists of three parts: “Model,” “Properties,” and “Simulator.” In the “Model” part, a given probabilistic system is described using the PRISM language. In the “Properties” part, the property specification language incorporates temporal logic such as PCTL (probabilistic computation tree logic) [29], and we can verify if a given PCTL formula holds. In the “Simulator,” the state trajectories can be simulated.Now, using PRISM, we propose a method for modeling a given CS-PBN. By modeling a given CS-PBN via PRISM, Problem2 can be solved. In the PRISM-based method, Boolean functions in a given PBN can be directly used. To explain the PRISM-based method, consider the PBN (5)–(7) in Example 1 and q=0.5. Suppose that the initial state and the initial Boolean function are given by x0=[111]⊤, f(1)(0)=f1(1), f(2)(0)=f1(2), and f(3)(0)=f1(3) (i.e., j0(1)=j0(2)=j0(3)=1), respectively. By using Lemma 3, each Boolean function can be transformed into some polynomial on the field of real numbers. Then, the PRISM code describing this CS-PBN is shown in Figure 3.Figure 3
PRISM code expressing the CS-PBN.In line 1, it is described that a given system is a MDP; that is, the control input (in other words, the nondeterministic variable) that must decide is included. In line 2, the probabilityq is given by q=0.5. In lines 3–7, the discrete-time Markov chain such as Figure 2 is modeled for f(1). The probabilistic variable  d1 corresponds to j∈{1,2} in fj(i). In line 4, f(1)(0) is given by f(1)(0)=f1(1). In lines 5-6, the behavior of  d1 is modeled. In line 5, it is described that if d1=1 holds, then the next state d1′ is 1 with the probability 0.6q+(1-q) and 2 with the probability 0.4q. In lines 8–12, f(1) is modeled. In line 9, it is described that x1 takes a binary value, and the initial value of x1 is given by 1. In line 10, f1(1) is modeled. In line 11, f2(1) is modeled. In a similar way, f(2) is modeled in lines 13–22, and f(3) is modeled in lines 23–32. In CS-PBNs, a discrete probabilistic distribution is given for each f(i). Hence, f(i), i=1,2,3, must be modeled separately. To associate with each module,  [CSPBN] is described. Finally, in lines 33–39, the property of the control input is described as a nondeterministic variable. Note that the initial value of the control input u(0) must be given (see line 34). Hence, PRISM must be executed for two cases of u(0)=0 and u(0)=1.The above explanation is the outline of the PRISM-based modeling method. Based on the above example, we propose a procedure for deriving the PRISM code expressing a general CS-PBN.


#### 2.2.1. Procedure for Modeling CS-PBNs

Step 1.
Transform each Boolean function into a polynomial on the real number field by using Lemma3. The obtained Boolean functions are denoted by f^j(i).Step 2.
Describe that a given system is a MDP, and giveq.Step 3.
Describe modules CSPBNmi and CSPBN i, i=1,2,…,n, as follows:module CSPBNmi
di  : [1..l(i)] init j0(i);[CSPBN]di=1→c1(i)q+(1-q):(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q:(di′=l(i));[CSPBN]di=2→c1(i)q:(di′=1)+c2(i)q+(1-q):(di′=2)+⋯+cl(i)(i)q:(di′=l(i));⋮[CSPBN]di=l(i)→c1(i)q:(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q+(1-q):(di′=l(i));endmodulemodule CSPBNi
xi:[0..1] init xi(0);[CSPBN]di=1→1.0: (xi′=f^1(i)(x,u));[CSPBN]di=2→1.0: (xi′=f^2(i)(x,u));⋮[CSPBN]di=l(i)→1.0: (xi′=f^l(i)(i)(x,u));endmoduleStep 4.
Describe the control inputui, i=1,2,…,m, as follows:module inputi
ui:[0..1] init ui(0);[PBN1]ui=0→(ui′=0)[PBN1]ui=0→(ui′=1)[PBN1]ui=1→(ui′=0)[PBN1]ui=1→(ui′=1)endmoduleFinally, consider solving Problem2. For solving this problem, we use Pmax prepared in PRISM. For example, suppose that y=[y1y2]⊤ and yf=[01]⊤. Then, in PRISM, Problem 2 can be described by
(8)Pmax=?[F<=N(y1=0)&(y2=1)].
Therefore, we see that Problem 2 can be solved using PRISM. The control input sequence u(0),u(1),…,u(N-1) is obtained simultaneously, but in PRISM 4.0.3, the obtained control input sequence cannot be displayed except for the case of N=∞. In the case of N=∞, the discrete-time Markov chain can be obtained as the closed-loop system of a given CS-PBN.


### 2.3. Optimal Control Using Polynomial Optimization

Consider the following problem.Problem 4 (optimal control problem).
Suppose that, for CS-PBN, the initial statex(0)=x0, the initial Boolean function f(i)(0,x(0),u(0))=fj0(i)(i)(x(0),u(0))(j0(i)∈{1,2,…,l(i)}), and the control time N are given (u(0) is not given). Then, find a control input sequence u(0),u(1),…,u(N-1) minimizing the cost function
(9)J=E[∑k=0N-1‍{Qx(k)+Ru(k)}+Qfx(N)∣x(0)=x0],
where Q,Qf∈ℛ1×n, R∈ℛ1×m are weighting vectors whose element is a nonnegative real number and E[·∣·] denotes a conditional expected value.According to the following two reasons, the linear cost function (9) is appropriate. (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 GRNs, expression of a certain gene is frequently focused (see, e.g., [19]). That is, in the cost function, the quadratic term such as xi(k)xj(k), i≠j, is not necessary.For a PBN, the authors have derived the following recursive representation of the expected value of the state:(10)E[xi(k+1)]=∑j=1l(i)‍cj(i)f^j(i)(E[x(k)],u(k)),
where the condition x(0)=x0 in the expected value is omitted. See [23] for further details. In this paper, this representation is extended to that in CS-PBNs.First, we present a simple example. Consider the PBN (5)–(7) in Example 1. Suppose that the initial state and the initial Boolean function are given by x0=[111]⊤, f(1)(0)=f1(1), f(2)(0)=f1(2), and f(3)(0)=f1(3), respectively.Consider deriving the expected value of the state at timek=1. Suppose that u(0)=0. Since the Boolean function at time k=0 is given, the state at time k=1 is uniquely derived as
(11)E[x(1)∣x(0)=x0,u(0)=0]=x(1)=[001]⊤.
Hereafter, the condition such as x(0)=x0, u(0)=0 in the expected value is omitted. Next, consider deriving the expected value of the state at time k=2. We remark that, for each xi, the discrete-time Markov chain such as Figure 2 can be obtained. For example, the probability that f1(1) is selected at time k=1 is 0.6q+(1-q), and the probability that f2(1) is selected at time k=1 is 0.4q. Suppose that u(1)=1. Then, we can obtain
(12)E[x1(2)]=E[f(1)(1,x(1),u(1))]=(0.6q+(1-q))E[f1(1)(x(1),u(1))]+0.4qE[f2(1)(x(1),u(1))]=(0.6q+(1-q))·(1-0)·1+0.4q·0=0.6q+(1-q),E[x2(2)]=E[f(2)(1,x(1),u(1))]=(0.7q+(1-q))E[f1(2)(x(1),u(1))]+0.3qE[f2(2)(x(1),u(1))]=(0.7q+(1-q))·(1-0)·1+0.3q·0=0.7q+(1-q),E[x3(2)]=E[f(3)(1,x(1),u(1))]=(0.8q+(1-q))E[f1(3)(x(1),u(1))]+0.2qE[f2(3)(x(1),u(1))]=(0.8q+(1-q))·(0+1-0·1)+0.2q·1=1.0.
Finally, we remark that the probability that some Boolean function fj(i) is selected is time-varying. For example, the probability that f1(1) is selected at time k=2 is (0.6q+(1-q))2+0.4q·0.6q, and the probability that f2(1) is selected at time k=2 is (0.6q+(1-q))·0.4q+0.4q·(0.2q+(1-q)).Next, consider a general case. From the observation of the above example, we can obtain the following recursive representation:(13)E[xi(k+1)]=∑j=1l(i)‍dj(i)(k)f^j(i)(E[x(k)],u(k)),d(i)(k+1)=P(i)d(i)(k),
where d(i)(k)=[d1(i)(k)d2(i)(k)⋯dl(i)(i)(k)]⊤ and
(14)P(i)=[c1(i)q+(1-q)c2(i)q⋯cl(i)(i)qc1(i)qc2(i)q+(1-q)⋯cl(i)(i)q⋮⋮⋱⋮c1(i)qc2(i)q⋯cl(i)(i)q+(1-q)].
Therefore, Problem 4 can be reduced to the following polynomial optimization problem:
(15)findE[x(k+1)]∈ℛn,u(k)∈ℛm,findE[x(k+1)]l∈k=0,1,…,N-1,min⁡Costfunction(9),subjecttoSystem(13),i=1,2,…,n,x(0)=x0,ui(k)(ui(k)-1)=0.
The constraint ui(k)(ui(k)-1)=0 guarantees that u(k) is a binary variable. A polynomial optimization problem can be solved by using a suitable solver such as SparsePOP [30].


## 2.1. Context-Sensitive Probabilistic Boolean Networks

First, we introduce a probabilistic Boolean network (PBN). Consider the following PBN:(1)x1(k+1)=f(1)(k,x(k),u(k)),x2(k+1)=f(2)(k,x(k),u(k)),⋮xn(k+1)=f(n)(k,x(k),u(k)),
where x=[x1x2⋯xn]⊤∈{0,1}n is the state (e.g., the expression of genes), u=[u1u2⋯um]⊤∈{0,1}m is the control input (e.g., the expression of genes), that is, the value of u can be arbitrarily given, and k=0,1,2,… is the discrete time. For a fixed k∈{0,1,…}, f(i):{0,1,…}×{0,1}n×{0,1}m→{0,1}1 is a given Boolean function consisting of logical operators such as AND (∧), OR (∨), and NOT (¬). In deterministic Boolean networks, x(k+1) is uniquely determined for given k, x(k), and u(k). In PBNs, the candidates of f(i)(k,x(k),u(k)) are given, and for each xi, selecting one Boolean function is probabilistically independent at each time. Let fj(i)(x(k),u(k)), j=1,2,…,l(i), denote the candidates of f(i)(k,x(k),u(k)). The probability that fj(i)(x(k),u(k)) is selected is defined by
(2)cj(i):=Prob(f(i)(k,x(k),u(k))=fj(i)(x(k),u(k))).
Then, the following relation:
(3)∑j=1l(i)‍cj(i)=1
must be satisfied. 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., [25]).Example 1.
As a simple example, consider the following deterministic Boolean network of an apoptosis network [26, 27] (see also Figure 1):
(4)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. Although Boolean dynamics in the above system are synchronous, both synchronous and asynchronous dynamics will be included. From this viewpoint, we consider the following PBN induced by the above system:
(5)f(1)={f1(1)=¬x2(k)∧u(k),c1(1)=0.6,f2(1)=x1(k),c2(1)=0.4,(6)f(2)={f1(1)=¬x1(k)∧x3(k),c1(2)=0.7,f2(1)=x2(k),c2(2)=0.3,(7)f(3)={f1(3)=x2(k)∨u(k),c1(3)=0.8,f2(3)=x3(k),c2(3)=0.2,
where l(1)=l(2)=l(3)=2, and we give cj(i) satisfying ∑j=1l(i)‍cj(i)=1. In addition, all state trajectories can be expressed as the state transition diagram with 23 nodes.Figure 1
Simplified model of an apoptosis network. Activation (solid), inhibition (broken).In PBNs, we suppose that selecting one Boolean function is probabilistically independent at each time. However, it will be natural to consider the situation that switches of Boolean functions do not occur frequently. From this viewpoint, a context-sensitive PBN (CS-PBN) has been proposed in [13, 14]. In CS-PBNs, the deciding time of Boolean functions is also selected randomly. Hereafter, let q∈[0,1] denote the probability that Boolean functions are switched at time k, and a pair of the system (1) and q is called a CS-PBN.To compare CS-PBNs with PBNs, consider (7) as a simple example. In PBNs, a switch of f1(3) and f2(3) functions does not depend on the Boolean function at time k-1. In CS-PBNs, a switch of f1(3) and f2(3) is decided by the discrete-time Markov chain in Figure 2. In other words, this switch depends on the Boolean function at time k-1. Owing to this difference, a control/verification method for CS-PBNs cannot be directly derived from that for PBNs.Figure 2
Discrete-time Markov chain inf(3).


## 2.2. Verification Using Model Checking

First, the reachability problem is formulated as the verification problem studied in this paper. The reachability problem is one of the typical verification problems. For a given CS-PBN, the outputy(k)=[y1(k)y2(k)⋯yp(k)]⊤∈{0,1}p is defined, where yi=xj, i=1,2,…,p, j∈𝒥⊆{1,2,…,n}. We remark that the output does not mean the measured signal. First, the reachability problem is formulated as follows.Problem 2 (reachability problem).
Suppose that, for CS-PBN with the output, the initial statex(0)=x0, the initial Boolean function f(i)(0,x(0),u(0))=fj0(i)(i)(x(0),u(0)) (j0(i)∈{1,2,…,l(i)}), the control time N, and the target output yf are given (u(0) is not given). Then, find a maximum probability Pmax⁡ that y(k)=yf holds within time N by manipulating a control input sequence u(0),u(1),…,u(N-1).In the standard reachability problem, only terminal time is focused, and it is checked whethery(N)=yf holds or not. In this paper, we focus on not only terminal time N but also other times 0,1,…,N-1. Furthermore, since a CS-PBN has the control input, which can be regarded as a nondeterministic variable, we find a maximum probability satisfying the condition.Next, we will propose a solution method for Problem2. As a preparation, the following lemma [28] is introduced.Lemma 3.
Consider two binary variablesδ1,δ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. By using this lemma, a Boolean function can be transformed into a polynomial on the real number field.To solve Problem2, the probabilistic model checker PRISM [15] is used. PRISM supports a discrete-time Markov chain (DT-MC), a continuous-time Markov chain (CT-MC), and a Markov decision process (MDP). PRISM consists of three parts: “Model,” “Properties,” and “Simulator.” In the “Model” part, a given probabilistic system is described using the PRISM language. In the “Properties” part, the property specification language incorporates temporal logic such as PCTL (probabilistic computation tree logic) [29], and we can verify if a given PCTL formula holds. In the “Simulator,” the state trajectories can be simulated.Now, using PRISM, we propose a method for modeling a given CS-PBN. By modeling a given CS-PBN via PRISM, Problem2 can be solved. In the PRISM-based method, Boolean functions in a given PBN can be directly used. To explain the PRISM-based method, consider the PBN (5)–(7) in Example 1 and q=0.5. Suppose that the initial state and the initial Boolean function are given by x0=[111]⊤, f(1)(0)=f1(1), f(2)(0)=f1(2), and f(3)(0)=f1(3) (i.e., j0(1)=j0(2)=j0(3)=1), respectively. By using Lemma 3, each Boolean function can be transformed into some polynomial on the field of real numbers. Then, the PRISM code describing this CS-PBN is shown in Figure 3.Figure 3
PRISM code expressing the CS-PBN.In line 1, it is described that a given system is a MDP; that is, the control input (in other words, the nondeterministic variable) that must decide is included. In line 2, the probabilityq is given by q=0.5. In lines 3–7, the discrete-time Markov chain such as Figure 2 is modeled for f(1). The probabilistic variable  d1 corresponds to j∈{1,2} in fj(i). In line 4, f(1)(0) is given by f(1)(0)=f1(1). In lines 5-6, the behavior of  d1 is modeled. In line 5, it is described that if d1=1 holds, then the next state d1′ is 1 with the probability 0.6q+(1-q) and 2 with the probability 0.4q. In lines 8–12, f(1) is modeled. In line 9, it is described that x1 takes a binary value, and the initial value of x1 is given by 1. In line 10, f1(1) is modeled. In line 11, f2(1) is modeled. In a similar way, f(2) is modeled in lines 13–22, and f(3) is modeled in lines 23–32. In CS-PBNs, a discrete probabilistic distribution is given for each f(i). Hence, f(i), i=1,2,3, must be modeled separately. To associate with each module,  [CSPBN] is described. Finally, in lines 33–39, the property of the control input is described as a nondeterministic variable. Note that the initial value of the control input u(0) must be given (see line 34). Hence, PRISM must be executed for two cases of u(0)=0 and u(0)=1.The above explanation is the outline of the PRISM-based modeling method. Based on the above example, we propose a procedure for deriving the PRISM code expressing a general CS-PBN.


### 2.2.1. Procedure for Modeling CS-PBNs

Step 1.
Transform each Boolean function into a polynomial on the real number field by using Lemma3. The obtained Boolean functions are denoted by f^j(i).Step 2.
Describe that a given system is a MDP, and giveq.Step 3.
Describe modules CSPBNmi and CSPBN i, i=1,2,…,n, as follows:module CSPBNmi
di  : [1..l(i)] init j0(i);[CSPBN]di=1→c1(i)q+(1-q):(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q:(di′=l(i));[CSPBN]di=2→c1(i)q:(di′=1)+c2(i)q+(1-q):(di′=2)+⋯+cl(i)(i)q:(di′=l(i));⋮[CSPBN]di=l(i)→c1(i)q:(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q+(1-q):(di′=l(i));endmodulemodule CSPBNi
xi:[0..1] init xi(0);[CSPBN]di=1→1.0: (xi′=f^1(i)(x,u));[CSPBN]di=2→1.0: (xi′=f^2(i)(x,u));⋮[CSPBN]di=l(i)→1.0: (xi′=f^l(i)(i)(x,u));endmoduleStep 4.
Describe the control inputui, i=1,2,…,m, as follows:module inputi
ui:[0..1] init ui(0);[PBN1]ui=0→(ui′=0)[PBN1]ui=0→(ui′=1)[PBN1]ui=1→(ui′=0)[PBN1]ui=1→(ui′=1)endmoduleFinally, consider solving Problem2. For solving this problem, we use Pmax prepared in PRISM. For example, suppose that y=[y1y2]⊤ and yf=[01]⊤. Then, in PRISM, Problem 2 can be described by
(8)Pmax=?[F<=N(y1=0)&(y2=1)].
Therefore, we see that Problem 2 can be solved using PRISM. The control input sequence u(0),u(1),…,u(N-1) is obtained simultaneously, but in PRISM 4.0.3, the obtained control input sequence cannot be displayed except for the case of N=∞. In the case of N=∞, the discrete-time Markov chain can be obtained as the closed-loop system of a given CS-PBN.


## 2.2.1. Procedure for Modeling CS-PBNs

Step 1.
Transform each Boolean function into a polynomial on the real number field by using Lemma3. The obtained Boolean functions are denoted by f^j(i).Step 2.
Describe that a given system is a MDP, and giveq.Step 3.
Describe modules CSPBNmi and CSPBN i, i=1,2,…,n, as follows:module CSPBNmi
di  : [1..l(i)] init j0(i);[CSPBN]di=1→c1(i)q+(1-q):(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q:(di′=l(i));[CSPBN]di=2→c1(i)q:(di′=1)+c2(i)q+(1-q):(di′=2)+⋯+cl(i)(i)q:(di′=l(i));⋮[CSPBN]di=l(i)→c1(i)q:(di′=1)+c2(i)q:(di′=2)+⋯+cl(i)(i)q+(1-q):(di′=l(i));endmodulemodule CSPBNi
xi:[0..1] init xi(0);[CSPBN]di=1→1.0: (xi′=f^1(i)(x,u));[CSPBN]di=2→1.0: (xi′=f^2(i)(x,u));⋮[CSPBN]di=l(i)→1.0: (xi′=f^l(i)(i)(x,u));endmoduleStep 4.
Describe the control inputui, i=1,2,…,m, as follows:module inputi
ui:[0..1] init ui(0);[PBN1]ui=0→(ui′=0)[PBN1]ui=0→(ui′=1)[PBN1]ui=1→(ui′=0)[PBN1]ui=1→(ui′=1)endmoduleFinally, consider solving Problem2. For solving this problem, we use Pmax prepared in PRISM. For example, suppose that y=[y1y2]⊤ and yf=[01]⊤. Then, in PRISM, Problem 2 can be described by
(8)Pmax=?[F<=N(y1=0)&(y2=1)].
Therefore, we see that Problem 2 can be solved using PRISM. The control input sequence u(0),u(1),…,u(N-1) is obtained simultaneously, but in PRISM 4.0.3, the obtained control input sequence cannot be displayed except for the case of N=∞. In the case of N=∞, the discrete-time Markov chain can be obtained as the closed-loop system of a given CS-PBN.


## 2.3. Optimal Control Using Polynomial Optimization

Consider the following problem.Problem 4 (optimal control problem).
Suppose that, for CS-PBN, the initial statex(0)=x0, the initial Boolean function f(i)(0,x(0),u(0))=fj0(i)(i)(x(0),u(0))(j0(i)∈{1,2,…,l(i)}), and the control time N are given (u(0) is not given). Then, find a control input sequence u(0),u(1),…,u(N-1) minimizing the cost function
(9)J=E[∑k=0N-1‍{Qx(k)+Ru(k)}+Qfx(N)∣x(0)=x0],
where Q,Qf∈ℛ1×n, R∈ℛ1×m are weighting vectors whose element is a nonnegative real number and E[·∣·] denotes a conditional expected value.According to the following two reasons, the linear cost function (9) is appropriate. (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 GRNs, expression of a certain gene is frequently focused (see, e.g., [19]). That is, in the cost function, the quadratic term such as xi(k)xj(k), i≠j, is not necessary.For a PBN, the authors have derived the following recursive representation of the expected value of the state:(10)E[xi(k+1)]=∑j=1l(i)‍cj(i)f^j(i)(E[x(k)],u(k)),
where the condition x(0)=x0 in the expected value is omitted. See [23] for further details. In this paper, this representation is extended to that in CS-PBNs.First, we present a simple example. Consider the PBN (5)–(7) in Example 1. Suppose that the initial state and the initial Boolean function are given by x0=[111]⊤, f(1)(0)=f1(1), f(2)(0)=f1(2), and f(3)(0)=f1(3), respectively.Consider deriving the expected value of the state at timek=1. Suppose that u(0)=0. Since the Boolean function at time k=0 is given, the state at time k=1 is uniquely derived as
(11)E[x(1)∣x(0)=x0,u(0)=0]=x(1)=[001]⊤.
Hereafter, the condition such as x(0)=x0, u(0)=0 in the expected value is omitted. Next, consider deriving the expected value of the state at time k=2. We remark that, for each xi, the discrete-time Markov chain such as Figure 2 can be obtained. For example, the probability that f1(1) is selected at time k=1 is 0.6q+(1-q), and the probability that f2(1) is selected at time k=1 is 0.4q. Suppose that u(1)=1. Then, we can obtain
(12)E[x1(2)]=E[f(1)(1,x(1),u(1))]=(0.6q+(1-q))E[f1(1)(x(1),u(1))]+0.4qE[f2(1)(x(1),u(1))]=(0.6q+(1-q))·(1-0)·1+0.4q·0=0.6q+(1-q),E[x2(2)]=E[f(2)(1,x(1),u(1))]=(0.7q+(1-q))E[f1(2)(x(1),u(1))]+0.3qE[f2(2)(x(1),u(1))]=(0.7q+(1-q))·(1-0)·1+0.3q·0=0.7q+(1-q),E[x3(2)]=E[f(3)(1,x(1),u(1))]=(0.8q+(1-q))E[f1(3)(x(1),u(1))]+0.2qE[f2(3)(x(1),u(1))]=(0.8q+(1-q))·(0+1-0·1)+0.2q·1=1.0.
Finally, we remark that the probability that some Boolean function fj(i) is selected is time-varying. For example, the probability that f1(1) is selected at time k=2 is (0.6q+(1-q))2+0.4q·0.6q, and the probability that f2(1) is selected at time k=2 is (0.6q+(1-q))·0.4q+0.4q·(0.2q+(1-q)).Next, consider a general case. From the observation of the above example, we can obtain the following recursive representation:(13)E[xi(k+1)]=∑j=1l(i)‍dj(i)(k)f^j(i)(E[x(k)],u(k)),d(i)(k+1)=P(i)d(i)(k),
where d(i)(k)=[d1(i)(k)d2(i)(k)⋯dl(i)(i)(k)]⊤ and
(14)P(i)=[c1(i)q+(1-q)c2(i)q⋯cl(i)(i)qc1(i)qc2(i)q+(1-q)⋯cl(i)(i)q⋮⋮⋱⋮c1(i)qc2(i)q⋯cl(i)(i)q+(1-q)].
Therefore, Problem 4 can be reduced to the following polynomial optimization problem:
(15)findE[x(k+1)]∈ℛn,u(k)∈ℛm,findE[x(k+1)]l∈k=0,1,…,N-1,min⁡Costfunction(9),subjecttoSystem(13),i=1,2,…,n,x(0)=x0,ui(k)(ui(k)-1)=0.
The constraint ui(k)(ui(k)-1)=0 guarantees that u(k) is a binary variable. A polynomial optimization problem can be solved by using a suitable solver such as SparsePOP [30].


## 3. Results and Discussion

In this section, we present a numerical example on the WNT5A network. First, the WNT5A network is explained. Next, computational results are presented.


### 3.1. WNT5A Network

Consider the GRN with the gene WNT5A, which is related to melanoma. A Boolean network model is given by(16)x1(k+1)=¬x6(k),x2(k+1)=(¬x2(k)∧x4(k)∧x6(k))∨{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 by x1, the concentration level of the gene pirin by x2, the concentration level of the gene S100P by x3, the concentration level of the gene RET1 by x4, the concentration level of the gene MART1 by x5, the concentration level of the gene HADHB by x6, and the concentration level of the gene STC2 by x7. See [31] for further details.Next, suppose that the control inputu is given by x2 (the concentration level of the gene pirin), according to the discussion in [19]. By replacing x2 and x3,x4,…,x7 with u and x2,x3,…,x6, respectively, we can obtain the following model:
(17)x1(k+1)=fd(1)(x(k),u(k))=¬x5(k),x2(k+1)=fd(2)(x(k),u(k))=¬x6(k),x3(k+1)=fd(3)(x(k),u(k))=x3(k),x4(k+1)=fd(4)(x(k),u(k))=¬x6(k)∨u(k),x5(k+1)=fd(5)(x(k),u(k))=x2(k)∨x3(k),x6(k+1)=fd(6)(x(k),u(k))=x6(k)∨¬u(k).Furthermore, we add the probabilistic behavior as follows:(18)xi(k+1)={fd(i)(x(k),u(k)),withtheprobabilityc1,xi(k),withtheprobabilityc2,
where l(i)=2 holds. Thus, we can obtain the PBN model expressing a WNT5A network.


### 3.2. Computational Result on Verification

Consider solving Problem4. For the PBN (18), we assume that c1 and c2 are given by c1=0.5 and c2=0.5, respectively. In the WNT5A network, it is important to inhibit the concentration level x1 of the gene WNT5A [32]. From this fact, we set y=x1 and yf=0. The initial state is given by x0=[100100]⊤. The initial Boolean function is given by fd(i). In addition, we set N=5.Next, we show the computation result. Then, we can obtainPmax⁡=0.7215 for q=0.3, Pmax⁡=0.6587 for q=0.5, and Pmax⁡=0.6489 for q=0.7. It is desirable that Pmax⁡ is close to 1. Hence, we see that the performance is degraded for a larger q.


### 3.3. Computational Result on Optimal Control

Consider solving Problem4. For the PBN (18), we assume that c1 and c2 are given by c1=0.8 and c2=0.2, respectively. Since the concentration level x1 must be inhibited, the weights Q,Qf, and R in Problem 4 are given by
(19)Q=[100000],R=1,Qf=[1000000],
respectively. The initial state is given by x0=[110100]⊤. The initial Boolean function is given by fd(i). In addition, we set N=5 and q=0.3.Next, we show the computation result. By solving Problem4, we can obtain u(0)=u(1)=1, u(2)=u(3)=u(4)=0. The expected value of the state at each time is obtained as
(20)E[x(1)]=[110110]⊤,E[x(2)]=[0.0610110]⊤,E[x(3)]=[0.006110110]⊤,E[x(4)]=[0.00080.2200.2210.9866]⊤,E[x(5)]=[0.00010.044800.04480.33850.998]⊤.
Hence, we see that the concentration level x1 of the gene WNT5A is inhibited with time.In addition, the optimal valueJ* of the cost function was 5.23. For q=0.5 and q=0.7, we can obtain J*=5.41 and J*=5.57, respectively. From these values, we see that the performance is degraded for a larger q.


## 3.1. WNT5A Network

Consider the GRN with the gene WNT5A, which is related to melanoma. A Boolean network model is given by(16)x1(k+1)=¬x6(k),x2(k+1)=(¬x2(k)∧x4(k)∧x6(k))∨{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 by x1, the concentration level of the gene pirin by x2, the concentration level of the gene S100P by x3, the concentration level of the gene RET1 by x4, the concentration level of the gene MART1 by x5, the concentration level of the gene HADHB by x6, and the concentration level of the gene STC2 by x7. See [31] for further details.Next, suppose that the control inputu is given by x2 (the concentration level of the gene pirin), according to the discussion in [19]. By replacing x2 and x3,x4,…,x7 with u and x2,x3,…,x6, respectively, we can obtain the following model:
(17)x1(k+1)=fd(1)(x(k),u(k))=¬x5(k),x2(k+1)=fd(2)(x(k),u(k))=¬x6(k),x3(k+1)=fd(3)(x(k),u(k))=x3(k),x4(k+1)=fd(4)(x(k),u(k))=¬x6(k)∨u(k),x5(k+1)=fd(5)(x(k),u(k))=x2(k)∨x3(k),x6(k+1)=fd(6)(x(k),u(k))=x6(k)∨¬u(k).Furthermore, we add the probabilistic behavior as follows:(18)xi(k+1)={fd(i)(x(k),u(k)),withtheprobabilityc1,xi(k),withtheprobabilityc2,
where l(i)=2 holds. Thus, we can obtain the PBN model expressing a WNT5A network.


## 3.2. Computational Result on Verification

Consider solving Problem4. For the PBN (18), we assume that c1 and c2 are given by c1=0.5 and c2=0.5, respectively. In the WNT5A network, it is important to inhibit the concentration level x1 of the gene WNT5A [32]. From this fact, we set y=x1 and yf=0. The initial state is given by x0=[100100]⊤. The initial Boolean function is given by fd(i). In addition, we set N=5.Next, we show the computation result. Then, we can obtainPmax⁡=0.7215 for q=0.3, Pmax⁡=0.6587 for q=0.5, and Pmax⁡=0.6489 for q=0.7. It is desirable that Pmax⁡ is close to 1. Hence, we see that the performance is degraded for a larger q.


## 3.3. Computational Result on Optimal Control

Consider solving Problem4. For the PBN (18), we assume that c1 and c2 are given by c1=0.8 and c2=0.2, respectively. Since the concentration level x1 must be inhibited, the weights Q,Qf, and R in Problem 4 are given by
(19)Q=[100000],R=1,Qf=[1000000],
respectively. The initial state is given by x0=[110100]⊤. The initial Boolean function is given by fd(i). In addition, we set N=5 and q=0.3.Next, we show the computation result. By solving Problem4, we can obtain u(0)=u(1)=1, u(2)=u(3)=u(4)=0. The expected value of the state at each time is obtained as
(20)E[x(1)]=[110110]⊤,E[x(2)]=[0.0610110]⊤,E[x(3)]=[0.006110110]⊤,E[x(4)]=[0.00080.2200.2210.9866]⊤,E[x(5)]=[0.00010.044800.04480.33850.998]⊤.
Hence, we see that the concentration level x1 of the gene WNT5A is inhibited with time.In addition, the optimal valueJ* of the cost function was 5.23. For q=0.5 and q=0.7, we can obtain J*=5.41 and J*=5.57, respectively. From these values, we see that the performance is degraded for a larger q.


## 4. Conclusions

In this paper, we discussed verification and optimal control for a context-sensitive probabilistic Boolean network (CS-PBN), which is one of the models for gene regulatory networks (GRNs). In verification, the PRISM-based method for PBNs [16] was extended to that for CS-PBNs. In optimal control, the optimal control method for PBNs [23] was extended to that for CS-PBNs. A CS-PBN is a generalized version of a PBN, and it enables us to consider several situations. Furthermore, as a numerical example, we considered the WNT5A network, which is related to melanoma. The proposed methods provide us useful tools in control theory of GRNs.In recent years, a stochastic Boolean network [33] has been proposed as a new representation of PBNs. In addition, to simplify a given Boolean network, the Karnaugh map realization of a Boolean network has been proposed in [34]. These modeling methods will be useful for reducing the computational burden. Future efforts will focus on applying these modeling methods to the control problem and the verification for CS-PBNs.



---

*Source: 968341-2014-01-23.xml*