Verification and optimal control of context-sensitive probabilistic Boolean networks using model checking and polynomial optimization.

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.

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 2 nodes (i.e., 2 × 2 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 Section 2.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} 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: where is the state (e.g., the expression of genes), 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 ( : {0,1,…}×{0,1} × {0,1} → {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 ((k, x(k), u(k)) are given, and for each x , selecting one Boolean function is probabilistically independent at each time. Let f ((x(k), u(k)), j = 1,2,…, l(i), denote the candidates of f ((k, x(k), u(k)). The probability that f ((x(k), u(k)) is selected is defined by Then, the following relation: 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): where the concentration level (high or low) of the inhibitor of apoptosis proteins (IAP) is denoted by x 1, the concentration level of the active caspase 3 (C3a) by x 2, and the concentration level of the active caspase 8 (C8a) by x 3. 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: where l(1) = l(2) = l(3) = 2, and we give c ( satisfying ∑ c ( = 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 f 1 (3) and f 2 (3) functions does not depend on the Boolean function at time k − 1. In CS-PBNs, a switch of f 1 (3) and f 2 (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 in f (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 output is defined, where y = x , 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 state x(0) = x 0, the initial Boolean function f ((0, x(0), u(0)) = f ((x(0), u(0)) (j 0(i)∈{1,2,…, l(i)}), the control time N, and the target output y are given (u(0) is not given). Then, find a maximum probability P max⁡ that y(k) = y 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 whether y(N) = y 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 Problem 2. As a preparation, the following lemma [28] is introduced.

Lemma 3

Consider two binary variables δ 1, δ 2. Then, the following relations hold.

¬δ 1 is equivalent to 1 − δ 1.

δ 1∨δ 2 is equivalent to δ 1 + δ 2 − δ 1 δ 2.

δ 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 Problem 2, 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, Problem 2 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 , f (1)(0) = f 1 (1), f (2)(0) = f 1 (2), and f (3)(0) = f 1 (3) (i.e., j 0(1) = j 0(2) = j 0(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 probability q 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 f (. In line 4, f (1)(0) is given by f (1)(0) = f 1 (1). In lines 5-6, the behavior of d1 is modeled. In line 5, it is described that if d 1 = 1 holds, then the next state d 1′ 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 x 1 takes a binary value, and the initial value of x 1 is given by 1. In line 10, f 1 (1) is modeled. In line 11, f 2 (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 (. Hence, f (, 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 Lemma 3. The obtained Boolean functions are denoted by .

Step 2

Describe that a given system is a MDP, and give q.

Step 3

Describe modules CSPBNm i and CSPBN i, i = 1,2,…, n, as follows:

module CSPBNm i

d   : [1..l(i)] init j 0(i);

[CSPBN] d = 1→c 1 ( q + (1 − q):(d ′ = 1) + c 2 ( q : (d ′ = 2)+⋯+c ( q : (d ′ = l(i));

[CSPBN] d = 2 → c 1 ( q : (d ′ = 1) + c 2 ( q + (1 − q):(d ′ = 2)+⋯+c ( q : (d ′ = l(i));

⋮

[CSPBN] d = l(i) → c 1 ( q : (d ′ = 1) + c 2 ( q : (d ′ = 2)+⋯+c ( q + (1 − q):(d ′ = l(i));

endmodule

module CSPBN i

x :[0..1] init x (0);

[CSPBN] d = 1→1.0: ();

[CSPBN] d = 2→1.0: ();

⋮

[CSPBN] d = l(i)→1.0: ();

endmodule

Step 4

Describe the control input u , i = 1,2,…, m, as follows:

module input i

u :[0..1] init u (0);

[PBN1] u = 0→(u ′ = 0)

[PBN1] u = 0→(u ′ = 1)

[PBN1] u = 1→(u ′ = 0)

[PBN1] u = 1→(u ′ = 1)

endmodule

Finally, consider solving Problem 2. For solving this problem, we use prepared in PRISM. For example, suppose that and . Then, in PRISM, Problem 2 can be described by 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 state x(0) = x 0, the initial Boolean function f ((0, x(0), u(0)) = f ((x(0), u(0))(j 0(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 where Q, Q ∈ ℛ 1×, R ∈ ℛ 1× 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 x 2(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 x (k)x (k), i≠j, is not necessary.

For a PBN, the authors have derived the following recursive representation of the expected value of the state: where the condition x(0) = x 0 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 , f (1)(0) = f 1 (1), f (2)(0) = f 1 (2), and f (3)(0) = f 1 (3), respectively.

Consider deriving the expected value of the state at time k = 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 Hereafter, the condition such as x(0) = x 0, 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 x , the discrete-time Markov chain such as Figure 2 can be obtained. For example, the probability that f 1 (1) is selected at time k = 1 is 0.6q + (1 − q), and the probability that f 2 (1) is selected at time k = 1 is 0.4q. Suppose that u(1) = 1. Then, we can obtain Finally, we remark that the probability that some Boolean function f ( is selected is time-varying. For example, the probability that f 1 (1) is selected at time k = 2 is (0.6q + (1 − q))2 + 0.4q · 0.6q, and the probability that f 2 (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: where and Therefore, Problem 4 can be reduced to the following polynomial optimization problem: The constraint u (k)(u (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 where the concentration level (high or low) of the gene WNT5A is denoted by x 1, the concentration level of the gene pirin by x 2, the concentration level of the gene S100P by x 3, the concentration level of the gene RET1 by x 4, the concentration level of the gene MART1 by x 5, the concentration level of the gene HADHB by x 6, and the concentration level of the gene STC2 by x 7. See [31] for further details.

Next, suppose that the control input u is given by x 2 (the concentration level of the gene pirin), according to the discussion in [19]. By replacing x 2 and x 3, x 4,…, x 7 with u and x 2, x 3,…, x 6, respectively, we can obtain the following model:

Furthermore, we add the probabilistic behavior as follows: where l(i) = 2 holds. Thus, we can obtain the PBN model expressing a WNT5A network.

3.2. Computational Result on Verification

Consider solving Problem 4. For the PBN (18), we assume that c 1 and c 2 are given by c 1 = 0.5 and c 2 = 0.5, respectively. In the WNT5A network, it is important to inhibit the concentration level x 1 of the gene WNT5A [32]. From this fact, we set y = x 1 and y = 0. The initial state is given by . The initial Boolean function is given by f (. In addition, we set N = 5.

Next, we show the computation result. Then, we can obtain P max⁡ = 0.7215 for q = 0.3, P max⁡ = 0.6587 for q = 0.5, and P max⁡ = 0.6489 for q = 0.7. It is desirable that P max⁡ 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 Problem 4. For the PBN (18), we assume that c 1 and c 2 are given by c 1 = 0.8 and c 2 = 0.2, respectively. Since the concentration level x 1 must be inhibited, the weights Q, Q , and R in Problem 4 are given by respectively. The initial state is given by . The initial Boolean function is given by f (. In addition, we set N = 5 and q = 0.3.

Next, we show the computation result. By solving Problem 4, 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 Hence, we see that the concentration level x 1 of the gene WNT5A is inhibited with time.

In addition, the optimal value J* 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.