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

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.

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 𝑛 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 𝑛-dimensional vectors, which consists of elements 0 and 1. For a matrix 𝑀, 𝑀^⊤ denotes the transpose of 𝑀.

2. Materials and Methods

2.1. Context-Sensitive Probabilistic Boolean Networks

First, we introduce a probabilistic Boolean network (PBN). Consider the following PBN: 𝑥₁(𝑘+1)=𝑓⁽¹⁾((𝑘)(𝑘,𝑥𝑘,𝑥,𝑢))₂(𝑘+1)=𝑓⁽²⁾((𝑘)(𝑘,⋮𝑥𝑘,𝑥,𝑢))_𝑛(𝑘+1)=𝑓^(𝑛)(𝑘,𝑥(𝑘),𝑢(𝑘)), (1) where 𝑥=𝑥₁𝑥₂⋯𝑥_𝑛^⊤∈{0,1}^𝑛 is the state (e.g., the expression of genes), 𝑢=𝑢₁𝑢₂⋯𝑢_𝑚^⊤∈{0,1}^𝑚 is the control input (e.g., the expression of genes), that is, the value of 𝑢 can be arbitrarily given, and 𝑘=0,1,2,… is the discrete time. For a fixed 𝑘∈{0,1,…}, 𝑓^(𝑖)∶{0,1,…}×{0,1}^𝑛×{0,1}^𝑚→{0,1}¹ is a given Boolean function consisting of logical operators such as AND (∧), OR (∨), and NOT (¬). In deterministic Boolean networks, 𝑥(𝑘+1) is uniquely determined for given 𝑘, 𝑥(𝑘), and 𝑢(𝑘). In PBNs, the candidates of 𝑓^(𝑖)(𝑘,𝑥(𝑘),𝑢(𝑘)) are given, and for each 𝑥_𝑖, selecting one Boolean function is probabilistically independent at each time. Let 𝑓_𝑗^(𝑖)(𝑥(𝑘),𝑢(𝑘)), 𝑗=1,2,…,𝑙(𝑖), denote the candidates of 𝑓^(𝑖)(𝑘,𝑥(𝑘),𝑢(𝑘)). The probability that 𝑓_𝑗^(𝑖)(𝑥(𝑘),𝑢(𝑘)) is selected is defined by 𝑐_𝑗^(𝑖)∶=Prob𝑓^(𝑖)(𝑘,𝑥(𝑘),𝑢(𝑘))=𝑓_𝑗^(𝑖)(𝑥(𝑘),𝑢(𝑘)). (2) Then, the following relation: ^𝑙(𝑖)_𝑗=1𝑐_𝑗^(𝑖)=1 (3) 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): 𝑥₁(𝑘+1)=¬𝑥₂(𝑘)(𝑘),𝑥∧𝑢₂(𝑘+1)=¬𝑥₁(𝑘)∧𝑥₃(𝑘),𝑥₃(𝑘+1)=𝑥₂(𝑘)∨𝑢(𝑘), (4) where the concentration level (high or low) of the inhibitor of apoptosis proteins (IAP) is denoted by 𝑥₁, the concentration level of the active caspase 3 (C3a) by 𝑥₂, and the concentration level of the active caspase 8 (C8a) by 𝑥₃. The concentration level of the tumor necrosis factor (TNF, a stimulus) is denoted by 𝑢 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: 𝑓⁽¹⁾=𝑓₁⁽¹⁾=¬𝑥₂(𝑘)∧𝑢(𝑘),𝑐₁⁽¹⁾𝑓=0.6,₂₍₁₎=𝑥₁(𝑘),𝑐₂⁽¹⁾𝑓=0.4,₍₂₎=𝑓₁⁽¹⁾=¬𝑥₁(𝑘)∧𝑥₃(𝑘),𝑐₁⁽²⁾𝑓=0.7,₂₍₁₎=𝑥₂(𝑘),𝑐₂⁽²⁾𝑓=0.3,₍₃₎=𝑓₁⁽³⁾=𝑥₂(𝑘)∨𝑢(𝑘),𝑐₁⁽³⁾𝑓=0.8,₂₍₃₎=𝑥₃(𝑘),𝑐₂⁽³⁾=0.2, (5)(6)(7) where 𝑙(1)=𝑙(2)=𝑙(3)=2, and we give 𝑐_𝑗^(𝑖) satisfying ∑^𝑙(𝑖)_𝑗=1𝑐_𝑗^(𝑖)=1. In addition, all state trajectories can be expressed as the state transition diagram with 2³ 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 𝑞∈[0,1] denote the probability that Boolean functions are switched at time 𝑘, and a pair of the system (1) and 𝑞 is called a CS-PBN.

To compare CS-PBNs with PBNs, consider (7) as a simple example. In PBNs, a switch of 𝑓₁⁽³⁾ and 𝑓₂⁽³⁾ functions does not depend on the Boolean function at time 𝑘−1. In CS-PBNs, a switch of 𝑓₁⁽³⁾ and 𝑓₂⁽³⁾ is decided by the discrete-time Markov chain in Figure 2. In other words, this switch depends on the Boolean function at time 𝑘−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 𝑓(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 𝑦(𝑘)=𝑦₁(𝑘)𝑦₂(𝑘)⋯𝑦_𝑝(𝑘)^⊤∈{0,1}^𝑝 is defined, where 𝑦_𝑖=𝑥_𝑗, 𝑖=1,2,…,𝑝, 𝑗∈𝒥⊆{1,2,…,𝑛}. 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 𝑥(0)=𝑥₀, the initial Boolean function 𝑓^(𝑖)(0,𝑥(0),𝑢(0))=𝑓_𝑗^(𝑖)_0(𝑖)(𝑥(0),𝑢(0)) (𝑗₀(𝑖)∈{1,2,…,𝑙(𝑖)}), the control time 𝑁, and the target output 𝑦_𝑓 are given (𝑢(0) is not given). Then, find a maximum probability 𝑃_max that 𝑦(𝑘)=𝑦_𝑓 holds within time 𝑁 by manipulating a control input sequence 𝑢(0),𝑢(1),…,𝑢(𝑁−1).

In the standard reachability problem, only terminal time is focused, and it is checked whether 𝑦(𝑁)=𝑦_𝑓 holds or not. In this paper, we focus on not only terminal time 𝑁 but also other times 0,1,…,𝑁−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 𝛿₁,𝛿₂. Then, the following relations hold.(i)¬𝛿₁ is equivalent to 1−𝛿₁.(ii)𝛿₁∨𝛿₂ is equivalent to 𝛿₁+𝛿₂−𝛿₁𝛿₂.(iii)𝛿₁∧𝛿₂ is equivalent to 𝛿₁𝛿₂.

For example, 𝛿₁∨¬𝛿₂ is equivalently transformed into 𝛿₁+(1−𝛿₂)−𝛿₁(1−𝛿₂)=1−𝛿₂+𝛿₁𝛿₂. 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 𝑞=0.5. Suppose that the initial state and the initial Boolean function are given by 𝑥₀=111^⊤, 𝑓⁽¹⁾(0)=𝑓₁⁽¹⁾, 𝑓⁽²⁾(0)=𝑓₁⁽²⁾, and 𝑓⁽³⁾(0)=𝑓₁⁽³⁾ (i.e., 𝑗₀(1)=𝑗₀(2)=𝑗₀(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 𝑞 is given by 𝑞=0.5. In lines 3–7, the discrete-time Markov chain such as Figure 2 is modeled for 𝑓⁽¹⁾. The probabilistic variable d1 corresponds to 𝑗∈{1,2} in 𝑓_𝑗^(𝑖). In line 4, 𝑓⁽¹⁾(0) is given by 𝑓⁽¹⁾(0)=𝑓₁⁽¹⁾. In lines 5-6, the behavior of d1 is modeled. In line 5, it is described that if 𝑑₁=1 holds, then the next state 𝑑^′₁ is 1 with the probability 0.6𝑞+(1−𝑞) and 2 with the probability 0.4𝑞. In lines 8–12, 𝑓⁽¹⁾ is modeled. In line 9, it is described that 𝑥₁ takes a binary value, and the initial value of 𝑥₁ is given by 1. In line 10, 𝑓₁⁽¹⁾ is modeled. In line 11, 𝑓₂⁽¹⁾ is modeled. In a similar way, 𝑓⁽²⁾ is modeled in lines 13–22, and 𝑓⁽³⁾ is modeled in lines 23–32. In CS-PBNs, a discrete probabilistic distribution is given for each 𝑓^(𝑖). Hence, 𝑓^(𝑖), 𝑖=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 𝑢(0) must be given (see line 34). Hence, PRISM must be executed for two cases of 𝑢(0)=0 and 𝑢(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 𝑞.

Step 3. Describe modules CSPBNm 𝑖 and CSPBN 𝑖, 𝑖=1,2,…,𝑛, as follows: module CSPBNm 𝑖 𝑑_𝑖  : [1..𝑙(𝑖)] init 𝑗₀(𝑖); [CSPBN] 𝑑_𝑖=1→𝑐₁^(𝑖)𝑞+(1−𝑞)∶(𝑑^′_𝑖=1)+𝑐₂^(𝑖)𝑞∶(𝑑^′_𝑖=2)+⋯+𝑐^(𝑖)_𝑙(𝑖)𝑞∶(𝑑^′_𝑖=𝑙(𝑖)); [CSPBN] 𝑑_𝑖=2→𝑐₁^(𝑖)𝑞∶(𝑑^′_𝑖=1)+𝑐₂^(𝑖)𝑞+(1−𝑞)∶(𝑑^′_𝑖=2)+⋯+𝑐^(𝑖)_𝑙(𝑖)𝑞∶(𝑑^′_𝑖=𝑙(𝑖)); ⋮ [CSPBN] 𝑑_𝑖=𝑙(𝑖)→𝑐₁^(𝑖)𝑞∶(𝑑^′_𝑖=1)+𝑐₂^(𝑖)𝑞∶(𝑑^′_𝑖=2)+⋯+𝑐^(𝑖)_𝑙(𝑖)𝑞+(1−𝑞)∶(𝑑^′_𝑖=𝑙(𝑖)); endmodule module CSPBN 𝑖 𝑥_𝑖∶[0..1] init 𝑥_𝑖(0); [CSPBN] 𝑑_𝑖=1→1.0∶ (𝑥^′_𝑖=𝑓₁^(𝑖)(𝑥,𝑢)); [CSPBN] 𝑑_𝑖=2→1.0∶ (𝑥^′_𝑖=𝑓₂^(𝑖)(𝑥,𝑢)); ⋮ [CSPBN] 𝑑_𝑖=𝑙(𝑖)→1.0∶ (𝑥^′_𝑖=𝑓^(𝑖)_𝑙(𝑖)(𝑥,𝑢)); endmodule

Step 4. Describe the control input 𝑢_𝑖, 𝑖=1,2,…,𝑚, as follows: module input 𝑖 𝑢_𝑖∶[0..1] init 𝑢_𝑖(0); [PBN1] 𝑢_𝑖=0→(𝑢^′_𝑖=0) [PBN1] 𝑢_𝑖=0→(𝑢^′_𝑖=1) [PBN1] 𝑢_𝑖=1→(𝑢^′_𝑖=0) [PBN1] 𝑢_𝑖=1→(𝑢^′_𝑖=1) endmodule

Finally, consider solving Problem 2. For solving this problem, we use 𝙿𝚖𝚊𝚡 prepared in PRISM. For example, suppose that 𝑦=𝑦₁𝑦₂^⊤ and 𝑦_𝑓=01^⊤. Then, in PRISM, Problem 2 can be described by 𝙿𝚖𝚊𝚡=?𝙵<=𝙽(𝚢1=0)&(𝚢2=1). (8) Therefore, we see that Problem 2 can be solved using PRISM. The control input sequence 𝑢(0),𝑢(1),…,𝑢(𝑁−1) is obtained simultaneously, but in PRISM 4.0.3, the obtained control input sequence cannot be displayed except for the case of 𝑁=∞. In the case of 𝑁=∞, 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 𝑥(0)=𝑥₀, the initial Boolean function 𝑓^(𝑖)(0,𝑥(0),𝑢(0))=𝑓_𝑗^(𝑖)_0(𝑖)(𝑥(0),𝑢(0))(𝑗₀(𝑖)∈{1,2,…,𝑙(𝑖)}), and the control time 𝑁 are given (𝑢(0) is not given). Then, find a control input sequence 𝑢(0),𝑢(1),…,𝑢(𝑁−1) minimizing the cost function 𝐽=𝐸^𝑁−1_𝑘=0{𝑄𝑥(𝑘)+𝑅𝑢(𝑘)}+𝑄_𝑓𝑥(𝑁)∣𝑥(0)=𝑥₀, (9) where 𝑄,𝑄_𝑓∈ℛ^1×𝑛, 𝑅∈ℛ^1×𝑚 are weighting vectors whose element is a nonnegative real number and 𝐸[⋅∣⋅] 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 𝛿²=𝛿 holds. That is, in the cost function, the quadratic term such as 𝑥²_𝑖(𝑘) 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 𝑥_𝑖(𝑘)𝑥_𝑗(𝑘), 𝑖≠𝑗, is not necessary.

For a PBN, the authors have derived the following recursive representation of the expected value of the state: 𝐸𝑥_𝑖(𝑘+1)=^𝑙(𝑖)_𝑗=1𝑐_𝑗^(𝑖)𝑓_𝑗^(𝑖)(𝐸[𝑥(𝑘)],𝑢(𝑘)), (10) where the condition 𝑥(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 𝑥₀=111^⊤, 𝑓⁽¹⁾(0)=𝑓₁⁽¹⁾, 𝑓⁽²⁾(0)=𝑓₁⁽²⁾, and 𝑓⁽³⁾(0)=𝑓₁⁽³⁾, respectively.

Consider deriving the expected value of the state at time 𝑘=1. Suppose that 𝑢(0)=0. Since the Boolean function at time 𝑘=0 is given, the state at time 𝑘=1 is uniquely derived as 𝐸𝑥(1)∣𝑥(0)=𝑥₀,𝑢(0)=0=𝑥(1)=001^⊤. (11) Hereafter, the condition such as 𝑥(0)=𝑥₀, 𝑢(0)=0 in the expected value is omitted. Next, consider deriving the expected value of the state at time 𝑘=2. We remark that, for each 𝑥_𝑖, the discrete-time Markov chain such as Figure 2 can be obtained. For example, the probability that 𝑓₁⁽¹⁾ is selected at time 𝑘=1 is 0.6𝑞+(1−𝑞), and the probability that 𝑓₂⁽¹⁾ is selected at time 𝑘=1 is 0.4𝑞. Suppose that 𝑢(1)=1. Then, we can obtain 𝐸𝑥₁(2)=𝐸𝑓⁽¹⁾((1)(1=𝐸𝑓1,𝑥,𝑢))0.6𝑞+1−𝑞₁⁽¹⁾(𝑥(1)(1𝑓,𝑢))+0.4𝑞𝐸₂⁽¹⁾(𝑥(1)(1=⋅(),𝐸𝑥,𝑢))0.6𝑞+1−𝑞1−0⋅1+0.4𝑞⋅0=0.6𝑞+1−𝑞₂(2)=𝐸𝑓⁽²⁾((1)(1=𝐸𝑓1,𝑥,𝑢))0.7𝑞+1−𝑞₁⁽²⁾(𝑥(1)(1𝑓,𝑢))+0.3𝑞𝐸₂⁽²⁾(𝑥(1)(1=⋅(),𝐸𝑥,𝑢))0.7𝑞+1−𝑞1−0⋅1+0.3𝑞⋅0=0.7𝑞+1−𝑞₃(2)=𝐸𝑓⁽³⁾((1)(1=𝐸𝑓1,𝑥,𝑢))0.8𝑞+1−𝑞₁⁽³⁾(𝑥(1)(1𝑓,𝑢))+0.2𝑞𝐸₂⁽³⁾(𝑥(1)(1=⋅(),𝑢))0.8𝑞+1−𝑞0+1−0⋅1+0.2𝑞⋅1=1.0. (12) Finally, we remark that the probability that some Boolean function 𝑓_𝑗^(𝑖) is selected is time-varying. For example, the probability that 𝑓₁⁽¹⁾ is selected at time 𝑘=2 is (0.6𝑞+(1−𝑞))²+0.4𝑞⋅0.6𝑞, and the probability that 𝑓₂⁽¹⁾ is selected at time 𝑘=2 is (0.6𝑞+(1−𝑞))⋅0.4𝑞+0.4𝑞⋅(0.2𝑞+(1−𝑞)).

Next, consider a general case. From the observation of the above example, we can obtain the following recursive representation: 𝐸𝑥_𝑖(𝑘+1)=^𝑙(𝑖)_𝑗=1𝑑_𝑗^(𝑖)(𝑘)𝑓_𝑗^(𝑖)(𝐸[𝑥(𝑘(𝑘,𝑑)],𝑢))_(𝑖)(𝑘+1)=𝑃^(𝑖)𝑑^(𝑖)(𝑘), (13) where 𝑑^(𝑖)(𝑘)=𝑑₁^(𝑖)(𝑘)𝑑₂^(𝑖)(𝑘)⋯𝑑^(𝑖)_𝑙(𝑖)(𝑘)^⊤ and 𝑃^(𝑖)=𝑐₁^(𝑖)𝑞+1−𝑞𝑐₂^(𝑖)𝑞⋯𝑐^(𝑖)_𝑙(𝑖)𝑞𝑐₁^(𝑖)𝑞𝑐₂^(𝑖)𝑞+1−𝑞⋯𝑐^(𝑖)_𝑙(𝑖)𝑞𝑐⋮⋮⋱⋮_1(𝑖)𝑞𝑐₂^(𝑖)𝑞⋯𝑐^(𝑖)_𝑙(𝑖)𝑞+1−𝑞. (14) Therefore, Problem 4 can be reduced to the following polynomial optimization problem: find𝐸[𝑥(𝑘+1)]∈ℛ^𝑛,𝑢(𝑘)∈ℛ^𝑚,(9),()𝑥(0)𝑘=0,1,…,𝑁−1,minCostfunctionsubjecttoSystem13,𝑖=1,2,…,𝑛,=𝑥₀,𝑢_𝑖(𝑘)𝑢_𝑖(𝑘)−1=0. (15) The constraint 𝑢_𝑖(𝑘)(𝑢_𝑖(𝑘)−1)=0 guarantees that 𝑢(𝑘) 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 𝑥₁(𝑘+1)=¬𝑥₆(𝑘),𝑥₂(𝑘+1)=¬𝑥₂(𝑘)∧𝑥₄(𝑘)∧𝑥₆(𝑘)∨𝑥₂(𝑘)∧𝑥₄(𝑘)∨𝑥₆(𝑘),𝑥₃(𝑘+1)=¬𝑥₇(𝑘),𝑥₄(𝑘+1)=𝑥₄(𝑘),𝑥₅(𝑘+1)=𝑥₂(𝑘)∨¬𝑥₇(𝑘),𝑥₆(𝑘+1)=𝑥₃(𝑘)∨𝑥₄(𝑘),𝑥₇(𝑘+1)=¬𝑥₂(𝑘)∨𝑥₇(𝑘), (16) where the concentration level (high or low) of the gene WNT5A is denoted by 𝑥₁, the concentration level of the gene pirin by 𝑥₂, the concentration level of the gene S100P by 𝑥₃, the concentration level of the gene RET1 by 𝑥₄, the concentration level of the gene MART1 by 𝑥₅, the concentration level of the gene HADHB by 𝑥₆, and the concentration level of the gene STC2 by 𝑥₇. See [31] for further details.

Next, suppose that the control input 𝑢 is given by 𝑥₂ (the concentration level of the gene pirin), according to the discussion in [19]. By replacing 𝑥₂ and 𝑥₃,𝑥₄,…,𝑥₇ with 𝑢 and 𝑥₂,𝑥₃,…,𝑥₆, respectively, we can obtain the following model: 𝑥₁(𝑘+1)=𝑓_𝑑⁽¹⁾(𝑥(𝑘),𝑢(𝑘))=¬𝑥₅(𝑘),𝑥₂(𝑘+1)=𝑓_𝑑⁽²⁾(𝑥(𝑘),𝑢(𝑘))=¬𝑥₆(𝑘),𝑥₃(𝑘+1)=𝑓_𝑑⁽³⁾(𝑥(𝑘),𝑢(𝑘))=𝑥₃(𝑘),𝑥₄(𝑘+1)=𝑓_𝑑⁽⁴⁾(𝑥(𝑘),𝑢(𝑘))=¬𝑥₆(𝑘)(𝑘),𝑥∨𝑢₅(𝑘+1)=𝑓_𝑑⁽⁵⁾(𝑥(𝑘),𝑢(𝑘))=𝑥₂(𝑘)∨𝑥₃(𝑘),𝑥₆(𝑘+1)=𝑓_𝑑⁽⁶⁾(𝑥(𝑘),𝑢(𝑘))=𝑥₆(𝑘)∨¬𝑢(𝑘). (17)

Furthermore, we add the probabilistic behavior as follows: 𝑥_𝑖()=𝑓𝑘+1_𝑑(𝑖)(𝑥(𝑘),𝑢(𝑘)),withtheprobability𝑐₁,𝑥_𝑖(𝑘),withtheprobability𝑐₂, (18) where 𝑙(𝑖)=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 𝑐₁ and 𝑐₂ are given by 𝑐₁=0.5 and 𝑐₂=0.5, respectively. In the WNT5A network, it is important to inhibit the concentration level 𝑥₁ of the gene WNT5A [32]. From this fact, we set 𝑦=𝑥₁ and 𝑦_𝑓=0. The initial state is given by 𝑥₀=100100^⊤. The initial Boolean function is given by 𝑓_𝑑^(𝑖). In addition, we set 𝑁=5.

Next, we show the computation result. Then, we can obtain 𝑃_max=0.7215 for 𝑞=0.3, 𝑃_max=0.6587 for 𝑞=0.5, and 𝑃_max=0.6489 for 𝑞=0.7. It is desirable that 𝑃_max is close to 1. Hence, we see that the performance is degraded for a larger 𝑞.

3.3. Computational Result on Optimal Control

Consider solving Problem 4. For the PBN (18), we assume that 𝑐₁ and 𝑐₂ are given by 𝑐₁=0.8 and 𝑐₂=0.2, respectively. Since the concentration level 𝑥₁ must be inhibited, the weights 𝑄,𝑄_𝑓, and 𝑅 in Problem 4 are given by 𝑄𝑄=100000,𝑅=1,_𝑓=1000000, (19) respectively. The initial state is given by 𝑥₀=110100^⊤. The initial Boolean function is given by 𝑓_𝑑^(𝑖). In addition, we set 𝑁=5 and 𝑞=0.3.

Next, we show the computation result. By solving Problem 4, we can obtain 𝑢(0)=𝑢(1)=1, 𝑢(2)=𝑢(3)=𝑢(4)=0. The expected value of the state at each time is obtained as 𝐸[𝑥(1)]=110110^⊤,𝐸[𝑥(2)]=0.0610110^⊤,𝐸[𝑥(3)]=0.006110110^⊤,𝐸[𝑥(4)]=0.00080.2200.2210.9866^⊤,𝐸[𝑥(5)]=0.00010.044800.04480.33850.998^⊤. (20) Hence, we see that the concentration level 𝑥₁ of the gene WNT5A is inhibited with time.

In addition, the optimal value 𝐽^∗ of the cost function was 5.23. For 𝑞=0.5 and 𝑞=0.7, we can obtain 𝐽^∗=5.41 and 𝐽^∗=5.57, respectively. From these values, we see that the performance is degraded for a larger 𝑞.

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.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was partially supported by Grant-in-Aid for Young Scientists (B) 23760387.