Optimization-Based Approaches to Control of Probabilistic Boolean Networks

Abstract

Control of gene regulatory networks is one of the fundamental topics in systems biology. In the last decade, control theory of Boolean networks (BNs), which is well known as a model of gene regulatory networks, has been widely studied. In this review paper, our previously proposed methods on optimal control of probabilistic Boolean networks (PBNs) are introduced. First, the outline of PBNs is explained. Next, an optimal control method using polynomial optimization is explained. The finite-time optimal control problem is reduced to a polynomial optimization problem. Furthermore, another finite-time optimal control problem, which can be reduced to an integer programming problem, is also explained.

1. Introduction

It is important to understand biological behavior as a system. For this purpose, systems biology has been established as a new research field [1,2]. One of the practical applications in systems biology is gene therapy technologies (see, e.g., [3]). In gene therapy, the expression of certain genes must be controlled, and such a problem can be formulated as a control problem. In other words, control of gene regulatory networks corresponds to therapeutic interventions realized by radiation, chemotherapy, and so on. Thus, control theory of gene regulatory networks is important for developing medical services in the future. Furthermore, from the viewpoint of experimental study, several results on control of gene regulatory networks have been obtained in [4,5,6,7]. For example, in [5], it was shown that cellular behavior can be regulated by control. These results suggest that feedback control of biological systems can be implemented. Motivated by the above background, we are developing theory of control methods for gene regulatory networks.

In developing a control method of gene regulatory networks, first, it is important to consider what mathematical model is utilized. A gene regulatory network is modeled by discrete dynamical systems (e.g., Bayesian networks, Boolean networks (BNs) [8]), continuous dynamical systems (e.g., differential/difference equations), or hybrid dynamical systems (e.g., piecewise affine models) (see, e.g., [9] for further details). In control theory, BNs and piecewise affine models are frequently utilized (see, e.g., [10] for control of a biological system modeled by a piecewise affine model). The control problem of hybrid systems is generally reduced to a mixed integer programming problem, in which continuous variables and binary variables are decision variables, and the computation time to solve it is frequently too long. Hence, a class of gene regulatory networks is limited to those with a small network. A BN is a mathematical model whose interactions between genes are expressed by Boolean functions [8]. That is, gene expression is approximated by a binary value (0 or 1). Since a BN is too simple, there is a possibility that biological behavior is not precisely modeled (see, e.g., [11]). On the other hand, this model can be relatively applied to large-scale systems. This advantage is focused on, and a BN is utilized in the first step of developing control theory of gene regulatory networks. Furthermore, a probabilistic Boolean network (PBN) [12] and a context-sensitive PBN [13,14] have been proposed as the extended models of BNs. In a PBN, at each discrete time, a Boolean function is randomly chosen from the candidates of Boolean functions. In a context-sensitive PBN, the time of choosing a Boolean function is also selected randomly. Hence, a context-sensitive PBN is a generalized version of a PBN.

In the last decade, many results on control of BNs and PBNs have been obtained so far (see, e.g., [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]). Especially, analysis and control methods using the semi-tensor product (STP) method have been widely studied (see, e.g., [33,34,35,36,37,38,39,40,41,42,43]). In the review paper [44], many methods for PBNs, e.g., construction and inference methods, structural intervention methods, and control methods, have been explained in a comprehensive manner. Of course, it is difficult to understand details of individual methods by reading only [44]. It is also significant to restrict topics and explain details. In this review paper, our previously proposed methods [23,24] in optimal external control of PBNs, i.e., a polynomial optimization approach and an integer programming approach are explained. Comparing our methods with other existing methods, the main advantage of our methods is that a relatively large-scale PBN can be handled. In many existing control methods (e.g., [19,20,21,30,32,36]), transition probability matrices with $2^n\times 2^n$ must be manipulated (n is the number of genes). As a result, a class of PBNs is restricted (see also Section 5 for further details). In our methods, transition probability matrices are not used. Moreover, since several solvers are available, optimization-based approaches are useful. To the best of our knowledge, such approaches have not been studied sufficiently, and their development in the future is necessary.

First, we introduce an optimal control method using polynomial optimization. In this method, the expected value of the state in PBN is modeled by a polynomial system, which can be directly obtained by Boolean functions and information on probabilities. Then, the finite-time optimal control problem is reduced to a polynomial optimization problem. Using e.g., SparsePOP [45], which is available on MATLAB R2006b or later, the polynomial optimization problem can be easily described and can be solved. However, there is a weakness that the polynomial optimization problem for large-scale PBNs cannot be solved owing to memory (but, comparing this method with other methods, a relatively large-scale PBN can be handled). In addition, the computational complexity of the finite-time optimal control problem of a PBN is ${\Sigma }_2^p$-hard [17], which is harder than NP-hard (Non-deterministic Polynomial time hard). Next, we introduce an optimal control method using integer programming. Instead of the expected value of the cost function, we focus on upper and lower bounds of the cost function under a certain constraint on probabilities. Then, the optimal control problem is reduced to an integer linear programming problem.

This review paper is organized as follows. In Section 2, PBNs are briefly explained. In Section 3, a control method using polynomial optimization is explained. In Section 4, a control method using integer programming is explained. In Section 5, discussion on numerical experiments is done. In Section 7, we conclude this review paper.

Notation: Let $\mathcal{R}$ 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/vector X, let $X^{\top }$ denote the transpose of X. Let I and 0 denote the identity matrix with an appropriate size and the zero matrix/vector/scalar, respectively. For an event A, let $P\left(A\right)$ denote the probability that A occurs. For two events $A,B$, let $P\left(A\right|B)$ denote the conditional probability of A given B. For two events $A,B$, let $E\left[A\right|B]$ denote the conditional expected value of A given B.

2. Probabilistic Boolean Networks

In this section, we introduce a probabilistic Boolean network (PBN).

Consider the following PBN:

$$\left\{\begin{array}{c}{x}_1(k+1)=f^{\left(1\right)}(x\left(k\right),u\left(k\right)),\hfill \\ x_2(k+1)=f^{\left(2\right)}(x\left(k\right),u\left(k\right)),\hfill \\ ⋮\hfill \\ x_n(k+1)=f^{\left(n\right)}(x\left(k\right),u\left(k\right)),\hfill \end{array}\right.$$

(1)

where $x={[x_1{x}_2\cdots x_n]}^{\top }\in {\{0,1\}}^n$ is the state (e.g., the expression of genes), and $u={[u_1{u}_2\cdots u_m]}^{\top }\in {\{0,1\}}^m$ is the control input (e.g., the expression of genes), i.e., the value of u can be arbitrarily determined, $k\in \{0,1,2,\cdots \}$ is the discrete time. The function $f^{\left(i\right)}:{\{0,1\}}^n\times {\{0,1\}}^m\to {\{0,1\}}^1$ is a given Boolean function consisting of logical operators such as AND (∧), OR (∨), XOR (⊕), and NOT (¬). In deterministic Boolean networks, $x(k+1)$ is uniquely determined for given $x\left(k\right)$ and $u\left(k\right)$. In PBNs, the candidates of $f^{\left(i\right)}(x\left(k\right),u\left(k\right))$ are given, and for each $x_i$, selecting one Boolean function is probabilistically independent at each time. The candidates of $f^{\left(i\right)}(x\left(k\right),u\left(k\right))$ are denoted by $f_j^{\left(i\right)}(x\left(k\right),u\left(k\right))$, $j=1,2,\cdots ,l\left(i\right)$, and the probability that $f_j^{\left(i\right)}(x\left(k\right),u\left(k\right))$ is selected is denoted by

$$c_j^{\left(i\right)}=P\left(f^{\left(i\right)}(x\left(k\right),u\left(k\right))=f_j^{\left(i\right)}(x\left(k\right),u\left(k\right))\right)\in [0,1].$$

Then, the relation

$$\sum _{j=1}^{l\left(i\right)}{c}_j^{\left(i\right)}=1$$

(2)

must be satisfied.

First, we present a biological example of (deterministic) Boolean networks (BNs).

Example 1.

Consider the following BN of an apoptosis network (see also Figure 1):

$$\left\{\begin{array}{c}{x}_1(k+1)=\neg x_2\left(k\right)\wedge u\left(k\right),\hfill \\ x_2(k+1)=\neg x_1\left(k\right)\wedge x_3\left(k\right),\hfill \\ x_3(k+1)=x_2\left(k\right)\vee u\left(k\right),\hfill \end{array}\right.$$

(3)

where $x_1$ is the concentration level (high or low) of the inhibitor of apoptosis proteins (IAP), $x_2$ is the concentration level of the active caspase 3 (C3a), and $x_3$ is the concentration level of the active caspase 8 (C8a). The control input u corresponds to the concentration level of the tumor necrosis factor (TNF, a stimulus). This model is introduced in [46], and is a simplified version of an apoptosis network model in [47]. In this model, $x_2\left(k\right)=0$ implies cell survival, and $x_1\left(k\right)=0,x_2\left(k\right)=1$ imply cell death.

Next, we present a simple example of PBNs.

Example 2.

Consider the PBN with three states and single control input. Boolean functions and probabilities are given by

$$\begin{array}{ccc}\hfill f^{\left(1\right)}& =& \left\{\begin{array}{c}{f}_1^{\left(1\right)}=x_3\left(k\right)\wedge u\left(k\right),c_1^{\left(1\right)}=0.8,\hfill \\ f_2^{\left(1\right)}=\neg x_3\left(k\right),c_2^{\left(1\right)}=0.2,\hfill \end{array}\right.\hfill \\ \hfill f^{\left(2\right)}& =& f_1^{\left(2\right)}=x_1\left(k\right)\wedge \neg x_3\left(k\right),c_1^{\left(2\right)}=1.0,\hfill \\ \hfill f^{\left(3\right)}& =& \left\{\begin{array}{c}{f}_1^{\left(3\right)}=x_1\left(k\right)\wedge \neg x_2\left(k\right),c_1^{\left(3\right)}=0.7,\hfill \\ f_2^{\left(3\right)}=x_2\left(k\right)\wedge u\left(k\right),c_2^{\left(3\right)}=0.3,\hfill \end{array}\right.\hfill \end{array}$$

where $l\left(1\right)=2$, $l\left(2\right)=1$ and $l\left(3\right)=2$ hold, and we see that the relation (2) is satisfied. Next, consider the state trajectories. Then, for $x\left(0\right)={\left[000\right]}^{\top }$ and $u\left(0\right)=1$, we can obtain

$$\begin{array}{ccc}\hfill P\left(x\left(1\right)={\left[000\right]}^{\top }{|x\left(0\right)=\left[000\right]}^{\top }\right)& =& 0.8,\hfill \\ \hfill P\left(x\left(1\right)={\left[100\right]}^{\top }{|x\left(0\right)=\left[000\right]}^{\top }\right)& =& 0.2.\hfill \end{array}$$

In this example, the cardinality of the finite state set ${\{0,1\}}^3$ is given by $2^3=8$, and we can obtain the state transition diagram of Figure 2 by computing the transition from each state under $u\left(0\right)=0$. In Figure 2, the number assigned to each node denotes $x_1$, $x_2$, $x_3$ (elements of the state), and the number assigned to each arc denotes the transition probability from some state to other state. Note here that, for simplicity, the state transitions from only $x\left(k\right)={\left[000\right]}^{\top },{\left[001\right]}^{\top },{\left[010\right]}^{\top },{\left[110\right]}^{\top }$ are illustrated in Figure 2.

3. Polynomial Optimization Approach

In this section, a polynomial optimization approach to optimal control of PBNs is reviewed. The advantage of this approach is that implementation is easy by using e.g., MATLAB.

3.1. Problem Formulation

For PBN (1), consider the following finite-time optimal control problem.

Problem 1.

Suppose that, for PBN (1), the initial state $x\left(0\right)=x_0$ and the control time N are given. Then, find a control input sequence $u\left(0\right),u\left(1\right),\cdots ,u(N-1)$ minimizing the following cost function:

$$J=E\left[\sum _{k=0}^{N-1}\left\{Qx\left(k\right)+Ru\left(k\right)\right\}+Q_{f}x\left(N\right)|x\left(0\right)=x_0\right],$$

(4)

where $Q,Q_f\in {\mathcal{R}}^{1\times n}$, $R\in {\mathcal{R}}^{1\times m}$ are weighting vectors whose element is a non-negative real number.

In this problem, we use some linear function with respect to the state and the control input as a cost function. Noting that the relation ${\delta }^2=\delta$ holds for a binary variable δ, it is unnecessary to add the quadratic term such as $x_i^2\left(k\right)$. Moreover, one of the main purposes in control of gene regulatory networks is that the expression of a certain gene is inhibited or activated (see e.g., [20]). Hence, simple linear cost function (4) is sufficient.

3.2. Transformation of a PBN into a Polynomial System

Consider modeling the dynamics of the expected value of the state by a polynomial system. First, the relation between a Boolean function and a polynomial is explained as a preparation. Next, the result for a general PBN is explained. Finally, a simple example is shown.

As a preparation, the following lemma [48] is introduced.

Lemma 1.

Consider two binary variables ${\delta }_1$ and ${\delta }_2$. Then, the following relations hold:

(i) 

$\neg {\delta }_1$ is equivalent to $1-{\delta }_1$,

(ii) 

${\delta }_1\wedge {\delta }_2$ is equivalent to ${\delta }_1{\delta }_2$,

(iii)

${\delta }_1\vee {\delta }_2$ is equivalent to ${\delta }_1+{\delta }_2-{\delta }_1{\delta }_2$,

(iv)

${\delta }_1\oplus {\delta }_2$ is equivalent to ${\delta }_1+{\delta }_2-2{\delta }_1{\delta }_2$.

By Lemma 1, a given Boolean function can be transformed into a polynomial on the real number field. For example, ${\delta }_1\vee \neg {\delta }_2$ is equivalently transformed into ${\delta }_1+(1-{\delta }_2)-{\delta }_1(1-{\delta }_2)=1-{\delta }_2+{\delta }_1{\delta }_2$.

Next, a recursive representation of a general PBN is explained. Let ${\widehat{f}}^{\left(i\right)}$ denote the polynomial that can be obtained from the Boolean function $f^{\left(i\right)}$. Let ${\widehat{f}}_j^{\left(i\right)}$ denote the polynomial that can be obtained from the Boolean function $f_j^{\left(i\right)}$. Hereafter, for simplicity of notation, the condition in $E\left[x_i\left(k\right)\right|*]$ is omitted.

In [24], the following theorem was obtained.

Theorem 1.

Suppose that, for the PBN (1), the initial state $x\left(0\right)=x_0$ is given. Then, the expected value of the state, $E\left[x_i\left(k\right)\right]$, is expressed as the following polynomial system:

$${\Sigma }_i:E\left[x_i(k+1)\right]=\sum _{j=1}^{l\left(i\right)}{c}_j^{\left(i\right)}{\widehat{f}}_j^{\left(i\right)}(E\left[x\left(k\right)\right],u\left(k\right)).$$

(5)

The outline of the proof is explained as follows. Since a Boolean function is independently switched for each state, $E\left[x_i{x}_j\right]=E\left[x_i\right]E\left[x_j\right]$, $i\ne j$ holds. In addition, there exists no second-order term of $x_i$ in ${\widehat{f}}_j^{\left(i\right)}$. Then, noting that $E\left[u\right(k\left)\right]=u\left(k\right)$ holds, we can obtain the system (5).

Noting that ${\widehat{f}}_j^{\left(i\right)}(x\left(k\right),u\left(k\right))$ is a polynomial, the right-hand side of the system (5) is also a polynomial. Therefore, from Theorem 1, we see that time evolution of the expected value of the state is modeled by a polynomial system.

We present a simple example.

Example 3.

Consider the PBN in Example 2 again. Suppose that, for this PBN, $x\left(0\right)$ and $u\left(0\right)$ are given by $x\left(0\right)=x_0={\left[000\right]}^{\top }$ and $u\left(0\right)=0$, respectively. Then, from Figure 2, we can obtain $E\left[x_1\left(1\right)\right]=0.8·0+0.2·1=0.2$. We can also obtain it in other way. By Lemma 1, $x_3\left(k\right)\wedge u\left(k\right)$ and $\neg x_3\left(k\right)$ are transformed into $x_3\left(k\right)u\left(k\right)$ and $1-x_3\left(k\right)$, respectively. Hence, we can obtain

$$\begin{array}{ccc}\hfill E\left[x_1\left(1\right)\right]& =& 0.8\left(x_3\left(0\right)u\left(0\right)\right)+0.2(1-x_3\left(0\right))\hfill \\ \hfill & =& 0.2.\hfill \end{array}$$

In a similar way, $E\left[x_2\left(1\right)\right]$ and $E\left[x_3\left(1\right)\right]$ can be obtained as

$$\begin{array}{ccc}\hfill E\left[x_2\left(1\right)\right]& =& x_1\left(0\right)(1-x_3\left(0\right))\hfill \\ \hfill & =& 0,\hfill \\ \hfill E\left[x_3\left(1\right)\right]& =& 0.7\left\{x_1\left(0\right)(1-x_2\left(0\right))\right\}+0.3\left(x_2\left(0\right)u\left(0\right)\right)\hfill \\ \hfill & =& 0,\hfill \end{array}$$

respectively.

Next, suppose that $u\left(1\right)$ is given by $u\left(1\right)=0$, and consider deriving $E\left[x\right(2\left)\right]$. Then, noting that a Boolean function is independently switched for each state, we can obtain

$$\begin{array}{ccc}\hfill E\left[x_1\left(2\right)\right]& =& E\left[f^{\left(1\right)}(x\left(1\right),u\left(1\right))\right]\hfill \\ \hfill & =& 0.8E\left[f_1^{\left(1\right)}(x\left(1\right),u\left(1\right))\right]+0.2E\left[f_2^{\left(1\right)}(x\left(1\right),u\left(1\right))\right]\hfill \\ \hfill & =& 0.8E\left[x_3\left(1\right)u\left(1\right)\right]+0.2E[1-x_3\left(1\right)]\hfill \\ \hfill & =& 0.8E\left[x_3\left(1\right)\right]u\left(1\right)+0.2-0.2E\left[x_3\left(1\right)\right]\hfill \\ \hfill & =& 0.2,\hfill \\ \hfill E\left[x_2\left(2\right)\right]& =& E\left[f^{\left(2\right)}(x\left(1\right),u\left(1\right))\right]\hfill \\ \hfill & =& 1.0E\left[x_1\left(1\right)(1-x_3\left(1\right))\right]\hfill \\ \hfill & =& E\left[x_1\left(1\right)\right]-E\left[x_1\left(1\right)\right]E\left[x_3\left(1\right)\right]\hfill \\ \hfill & =& 0.2,\hfill \\ \hfill E\left[x_3\left(2\right)\right]& =& E\left[f^{\left(3\right)}(x\left(1\right),u\left(1\right))\right]\hfill \\ \hfill & =& 0.7E[x_1\left(1\right)-x_1\left(1\right)x_2\left(1\right)]+0.3E\left[x_2\left(1\right)u\left(1\right)\right]\hfill \\ \hfill & =& 0.14.\hfill \end{array}$$

By recursively repeating, we can obtain $E\left[x\right(k\left)\right]$, $k\ge 3$.

From this example, we see that the expected value of the state can be expressed as a polynomial system using Theorem 1.

3.3. Reduction to a Polynomial Optimization Problem

Consider reducing Problem 1 to a polynomial optimization problem. Then, the following theorem was obtained in [24].

Theorem 2.

Problem 1 is equivalent to the following polynomial optimization problem:

$$\begin{array}{ccc}\hfill \mathrm{find}& & E\left[x(k+1)\right]\in {\mathcal{R}}^n,u\left(k\right)\in {\mathcal{R}}^m,k=0,1,\cdots ,N-1\hfill \\ \hfill min& & \mathrm{Cost}\mathrm{function} \left(4\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& & \mathrm{System}{\Sigma }_{\mathrm{i}},i=1,2,\cdots ,n,\hfill \\ \hfill & & x\left(0\right)=x_0,\hfill \\ \hfill & & u_i\left(k\right)(u_i\left(k\right)-1)=0.\hfill \end{array}$$

This theorem is immediately obtained from the following two properties: (i) the relation $E\left[u\right(k\left)\right]=u\left(k\right)$ holds; (ii) the fact that $u_i\left(k\right)$ is a binary variable is equivalent to the constraint $u_i\left(k\right)(u_i\left(k\right)-1)=0$. We remark here that this constraint is non-convex. From the practical viewpoint, it is desirable to relax this constraint to $0\le u_i\left(k\right)\le 1$. Furthermore, from Theorem 1, we see that $E\left[x(k+1)\right]\in {[0,1]}^n$ always holds. Hence, we set $E\left[x(k+1)\right]\in {\mathcal{R}}^n$ in the above polynomial optimization problem. Finally, $E\left[x\right(k+1\left)\right]$ may be eliminated from decision variables in the above problem, but its elimination is not easy. Hence, $E\left[x\right(k+1\left)\right]$ is included in a decision variable of the above problem.

A polynomial optimization problem can be solved by using e.g., SparsePOP [45], which is available on MATLAB. In SparsePOP, polynomials can be directly described.

4. Integer Programming Approach

The effectiveness of the polynomial optimization approach is restricted because the computational complexity of Problem 1 is ${\Sigma }_2^p$-hard [17], which is harder than NP-hard. In this section, an integer programming approach is explained. Instead of the expected value of the cost function, we focus on upper and lower bounds of the cost function. Then, the optimal control problem is reduced to an integer linear programming (ILP) problem. This method is a simplified version of our previously proposed method in [23].

4.1. Problem Formulation

First, preparations are given. Suppose that for PBN (1), $\overline{j}(i,k)\in \{1,2,\cdots ,l\left(i\right)\}$ is given for fixed $i\in \{1,2,\cdots ,n\}$ and $k\in \{0,1,\cdots \}$. Then, we define

$${\pi }_{\overline{j}}\left(k\right):=\prod _{i=1}^n{c}_{\overline{j}(i,k)}^{\left(i\right)}.$$

The parameter ${\pi }_{\overline{j}}\left(k\right)$ implies the probability that the Boolean functions $f_{\overline{j}(1,k)}^{\left(1\right)},f_{\overline{j}(2,k)}^{\left(2\right)},\cdots ,f_{\overline{j}(n,k)}^{\left(n\right)}$ are selected at time k. Furthermore, we define

$${\pi }_{\overline{j}}(k_1,k_2):=\prod _{q=k_1}^{k_2}{\pi }_{\overline{j}}\left(q\right).$$

The parameter ${\pi }_{\overline{j}}(k_1,k_2)$ implies the probability that certain Boolean functions are selected in the time interval $[k_1,k_2]$.

Example 4.

Consider the PBN in Example 2 again. Suppose $\overline{j}(1,0)=1$, $\overline{j}(2,0)=1$, and $\overline{j}(3,0)=2$. In other words, $f_1^{\left(1\right)}$, $f_1^{\left(2\right)}$, $f_2^{\left(3\right)}$ are given as a Boolean function at time 0. Then, we can obtain

$${\pi }_{\overline{j}}\left(0\right)=c_1^{\left(1\right)}·c_1^{\left(2\right)}·c_2^{\left(3\right)}=0.8·1·0.3=0.24.$$

Furthermore, suppose that $\overline{j}(1,1)=2$, $\overline{j}(2,1)=1$, $\overline{j}(3,1)=1$ $\overline{j}(1,2)=1$, $\overline{j}(2,2)=1$, and $\overline{j}(3,2)=1$. Then, we can obtain

$$\begin{array}{ccc}\hfill {\pi }_{\overline{j}}\left(1\right)& =& 0.2·1·0.7=0.14,\hfill \\ \hfill {\pi }_{\overline{j}}\left(2\right)& =& 0.8·1·0.7=0.56,\hfill \\ \hfill {\pi }_{\overline{j}}(0,2)& =& {\pi }_{\overline{j}}\left(0\right)·{\pi }_{\overline{j}}\left(1\right)·{\pi }_{\overline{j}}\left(2\right)=0.018816.\hfill \end{array}$$

Under these preparations, two problems (i.e., the optimal control problem and the performance evaluation problem) are formulated. First, the optimal control problem is formulated as follows.

Problem 2.

Suppose that, for PBN (1), the initial state $x\left(0\right)=x_0$, ρ satisfying $0\le \rho \le 1$, the control time N, and the weighting vectors $Q,Q_f\in {\mathcal{R}}^{1\times n}$, $R\in {\mathcal{R}}^{1\times m}$ whose element is a non-negative real number are given. Then, solve the following two problems.

Problem 2-A:

Find a control input sequence $u\left(0\right)$, $u\left(1\right)$, ⋯, $u(N-1)$ minimizing the lower bound of the cost function

$$J=\sum _{k=0}^{N-1}\left\{Qx\left(k\right)+Ru\left(k\right)\right\}+Q_{f}x\left(N\right)$$

(6)

satisfying the following constraint

$${\pi }_{\overline{j}}(0,N-1)\ge \rho .$$

(7)

Problem 2-B:

Apply the control input sequence obtained in Problem 2-A to the PBN. Then, find an upper bound of the cost function (6).

Let ${\overline{J}}_s$ and ${\underline{J}}_s$ denote the obtained upper and lower bounds, respectively. In a similar way to Problem 1, a linear function with respect to x and u is considered as a cost function. Using the target value of the state (the offset vector) $x_d\in {\{0,1\}}^n$, $x\left(i\right)$ may be replaced to $\widehat{x}\left(i\right):=x\left(i\right)-x_d$. Then, the cost function (6) must be changed to $J={\sum }_{i=0}^{N-1}\left\{|Q\widehat{x}\left(i\right)|+Ru\left(i\right)\right\}+\left|Q_f\widehat{x}\left(N\right)\right|$. Furthermore, by the inequality constraint (7), combinations of Boolean functions selected with the low probability can be excluded.

Each element of the control input sequence obtained solving Problem 2 is in general different. In gene regulatory networks, it is difficult to implement such control input sequence. In this case, we must impose constraints for control input sequences. For example, one of the simple methods is that the control input is given by always 1 (active) or 0 (inactive). That is, Problem 2-A is modified as follows: find $u_0(=u\left(0\right)=u\left(1\right)=\cdots =u(N-1))$ minimizing the lower bound of the cost function (6) under the constraint (7).

We present a simple example.

Example 5.

Consider the PBN with a single control input. Boolean functions and probabilities are given by

$$\begin{array}{ccc}\hfill f^{\left(1\right)}& =& \left\{\begin{array}{c}{f}_1^{\left(1\right)}=x_3\left(k\right)\wedge \neg u\left(k\right),c_1^{\left(1\right)}=0.8,\hfill \\ f_2^{\left(1\right)}=\neg x_3\left(k\right)\wedge \neg u\left(k\right),c_2^{\left(1\right)}=0.2,\hfill \end{array}\right.\hfill \\ \hfill f^{\left(2\right)}& =& f_1^{\left(2\right)}=x_1\left(k\right)\wedge \neg x_3\left(k\right),c_1^{\left(2\right)}=1.0,\hfill \\ \hfill f^{\left(3\right)}& =& \left\{\begin{array}{c}{f}_1^{\left(3\right)}=x_1\left(k\right)\wedge \neg x_2\left(k\right),c_1^{\left(3\right)}=0.7,\hfill \\ f_2^{\left(3\right)}=x_2\left(k\right)\vee u\left(k\right),c_2^{\left(3\right)}=0.3,\hfill \end{array}\right.\hfill \end{array}$$

respectively, where $u\left(k\right)\in {\{0,1\}}^1$. The initial state and the cost function are given by

$$x\left(0\right)=x_0=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right],J=\sum _{k=0}^1\left\{x_3\left(k\right)+u\left(k\right)\right\}+x_3\left(2\right),$$

respectively. For simplicity of discussion, ρ is given by $\rho =0$, i.e., the constraint (7) is not imposed.

Consider calculating the cost function explicitly. First, $x_3\left(0\right)=1$ holds. Next, we can obtain the following relations at time 1:

$$\begin{array}{ccc}\hfill & & \left\{\begin{array}{c}P(x_1\left(1\right|x\left(0\right)=x_0)=\neg u\left(0\right)=0.8,\hfill \\ P(x_1\left(1\right)=0|x\left(0\right)=x_0)=0.2,\hfill \end{array}\right.\hfill \\ \hfill & & P(x_2\left(1\right)=0|x\left(0\right)=x_0)=1.0,\hfill \\ \hfill & & \left\{\begin{array}{c}P(x_3\left(1\right)=0|x\left(0\right)=x_0)=0.7,\hfill \\ P(x_3\left(1\right)=1|x\left(0\right)=x_0)=0.3.\hfill \end{array}\right.\hfill \end{array}$$

Furthermore, $x_3\left(2\right)$ can be derived as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}P(x_3\left(2\right)=\neg u\left(0\right)|x\left(0\right)=x_0)=0.7·0.8=0.56,\hfill \\ P(x_3\left(2\right)=0|x\left(0\right)=x_0)=0.7·0.2=0.14,\hfill \\ P(x_3\left(2\right)=u\left(1\right)|x\left(0\right)=x_0)=0.3·1.0=0.3.\hfill \end{array}\right.\end{array}$$

Noting that $\neg u\left(0\right)$ is equivalent to $1-u\left(0\right)$ (see Lemma 1), we can obtain

$$\begin{array}{ccc}\hfill & & P(J=2+u\left(1\right)|x\left(0\right)=x_0)=0.392,\hfill \\ \hfill & & P(J=1+u\left(0\right)+2u\left(1\right)|x\left(0\right)=x_0)=0.210,\hfill \\ \hfill & & P(J=1+u\left(0\right)+u\left(1\right)|x\left(0\right)=x_0)=0.098,\hfill \\ \hfill & & P(J=3+u\left(1\right)|x\left(0\right)=x_0)=0.168,\hfill \\ \hfill & & P(J=2+u\left(0\right)+2u\left(1\right)|x\left(0\right)=x_0)=0.090,\hfill \\ \hfill & & P(J=2+u\left(0\right)+u\left(1\right)|x\left(0\right)=x_0)=0.042.\hfill \end{array}$$

First, consider solving Problem 2. We can obtain the optimal solution of Problem 2-A as $u\left(0\right)=u\left(1\right)=0$. The minimum value of the lower bound of the cost function is ${\underline{J}}_s=1$. The minimum value of the upper bounds can be obtained as ${\overline{J}}_s=3$ (this is the solution of Problem 2-B). Next, consider the case of $u\left(0\right)=u\left(1\right)=1$ as an example of other constant input. Then, the upper and lower bound of the cost function can be obtained by 2 and 5, respectively. From this result, we see that a solution for Problem 2-A should be applied to the system.

Problem 2 for a simple PBN can be easily solved. However, solving this problem for a large-scale PBN is hard. Here, we introduce an integer programming approach for solving Problem 2.

4.2. Solution Method for Problem 2-A

In this subsection, a solution method for Problem 2-A is derived. The solution method for Problem 2-B will be explained in the last part of this subsection.

As a preparation, we explain one lemma. The product of multiple binary variables such as ${\delta }_1{\delta }_2$ can be linearized by using the following lemma [49].

Lemma 2.

Suppose that binary variables ${\delta }_j\in \{0,1\}$, $j\in \mathcal{J}$ are given, where $\mathcal{J}$ is some index set. Then, $z={\prod }_{j\in \mathcal{J}}{\delta }_j$ is equivalent to the following linear inequalities

$$\sum _{j\in \mathcal{J}}{\delta }_j-z\le \left|\mathcal{J}\right|-1,-\sum _{j\in \mathcal{J}}{\delta }_j+\left|\mathcal{J}\right|z\le 0,$$

where $\left|\mathcal{J}\right|$ is the cardinality of $\mathcal{J}$.

From Lemmas 1 and 2, we see that any Boolean function can be equivalently transformed into a pair of linear functions and linear inequalities. See [48,49] for further details.

Consider transforming Problem 2-A into an ILP problem. Using Lemma 1, $f_j^{\left(i\right)}(x\left(k\right),u\left(k\right))$ in the PBN (1) can be equivalently transformed into a certain polynomial. Let ${\widehat{f}}_j^{\left(i\right)}(x\left(k\right),u\left(k\right))$ denote the obtained polynomial. Then, consider the following system:

$$x_i(k+1)=\sum _{j=1}^{l\left(i\right)}\left\{{\delta }_{i,j}\left(k\right){\widehat{f}}_j^{\left(i\right)}(x\left(k\right),u\left(k\right))\right\},$$

(8)

where $i=1,2,\cdots ,n$ and ${\delta }_{i,j}\left(k\right)\in {\{0,1\}}^1$. For ${\delta }_{i,j}\left(k\right)$, we impose the following constraint:

$$\sum _{j=1}^{l\left(i\right)}{\delta }_{i,j}\left(k\right)=1.$$

(9)

Next, using $c_j^{\left(i\right)}$ and ${\delta }_{i,j}\left(k\right)$, we define

$$\begin{array}{ccc}\hfill L_i& :=& \left[\begin{array}{cccc}\ln c_1^{\left(i\right)}& \ln c_2^{\left(i\right)}& \cdots & \ln c_{l\left(i\right)}^{\left(i\right)}\end{array}\right],\hfill \\ \hfill L& :=& \left[\begin{array}{cccc}{L}_1& L_2& \cdots & L_n\end{array}\right],\hfill \\ \hfill {\delta }_i\left(k\right)& :=& {\left[\begin{array}{cccc}{\delta }_{i,1}& {\delta }_{i,2}& \cdots & {\delta }_{i,l\left(i\right)}\end{array}\right]}^{\top }\in {\{0,1\}}^{l\left(i\right)},\hfill \\ \hfill \delta \left(k\right)& :=& {\left[\begin{array}{cccc}{\delta }_1\left(k\right)& {\delta }_2\left(k\right)& \cdots & {\delta }_n\left(k\right)\end{array}\right]}^{\top }\in {\{0,1\}}^l,\hfill \end{array}$$

where $l:={\sum }_{i=1}^{n}l\left(i\right)$. Then, the following lemma was obtained in [23].

Lemma 3.

Problem 2-A is equivalent to the following problem.

Problem 2-C:

$$\begin{array}{ccc}\hfill \mathrm{find}& & u\left(k\right),\delta \left(k\right),k=0,1,\cdots ,N-1\hfill \\ \hfill min& & \mathrm{Cost}\mathrm{function} \left(6\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& & \mathrm{System}\left(8\right),x\left(0\right)=x_0,\hfill \\ \hfill & & \mathrm{Inequality}\mathrm{constraint}{\displaystyle \sum _{k=0}^{N-1}L\delta \left(k\right)}\ge \ln \rho ,\hfill \\ \hfill & & \mathrm{Equality}\mathrm{constraint} \left(9\right).\hfill \end{array}$$

The outline of the proof is explained as follows. By deciding $\delta \left(k\right)$, Boolean functions at each time are uniquely determined. Then, we can obtain $\ln {\pi }_{\overline{j}}\left(k\right)=L\delta \left(k\right)$. That is, using the logarithm of the probability, ${\pi }_{\overline{j}}\left(k\right)$ can be expressed by a linear form with respect to $\delta \left(k\right)$. Hence, we can obtain

$$\ln {\pi }_{\overline{j}}(0,N-1)=\sum _{k=0}^{N-1}\ln {\pi }_{\overline{j}}\left(k\right)=\sum _{k=0}^{N-1}L\delta \left(k\right).$$

Thus, Problem 2-A is equivalent to Problem 2-C.

We consider further rewriting Problem 2-C. Since the system (8) is a polynomial system, the system (8) can be equivalently transformed into a linear system by using Lemma 2. Hence, the system (8) can be transformed into the following system:

$$\left\{\begin{array}{c}x(k+1)=Ax\left(k\right)+B_{1}u\left(k\right)+B_2\delta \left(k\right)+B_{3}z\left(k\right),\hfill \\ Cx\left(k\right)+D_{1}u\left(k\right)+D_2\delta \left(k\right)+D_{3}z\left(k\right)\le E,\hfill \end{array}\right.$$

(10)

where $z\left(k\right)\in {\{0,1\}}^p$ is an auxiliary variable, and the dimension p can be determined from the number of the product of binary variables. For simplicity of notation, we define $B:=\left[B_1{B}_2{B}_3\right]$, $v\left(k\right):={\left[u^{\top }\left(k\right){\delta }^{\top }\left(k\right)z^{\top }\left(k\right)\right]}^{\top }$. In addition, under $x\left(0\right)=x_0\in {\{0,1\}}^n$ and $v\left(k\right)\in {\{0,1\}}^{m+l+p}$, $x\left(k\right)$ becomes a binary variable automatically. Hence, we set $x\left(k\right)\in {\mathcal{R}}^n$, $k\ge 1$.

We present a simple example.

Example 6.

Consider the PBN in Example 5 again. Then, the system (8) can be derived as

$$\begin{array}{ccc}\hfill x_1(k+1)& =& {\delta }_{1,1}\left(k\right)\left\{x_3\left(k\right)(1-u\left(k\right))\right\}+{\delta }_{1,2}\left(k\right)\left\{(1-x_3\left(k\right))(1-u\left(k\right))\right\},\hfill \\ \hfill x_2(k+1)& =& {\delta }_{2,1}\left(k\right)\left\{x_1\left(k\right)(1-x_3\left(k\right))\right\},\hfill \\ \hfill x_3(k+1)& =& {\delta }_{3,1}\left(k\right)\left\{x_1\left(k\right)(1-x_2\left(k\right))\right\}+{\delta }_{3,2}\left(k\right)\left\{x_2\left(k\right)+u\left(k\right)-x_2\left(k\right)u\left(k\right)\right\}.\hfill \end{array}$$

The equality constraint (9) can be derived as

$$\sum _{j=1}^2{\delta }_{1,j}\left(k\right)=1,{\delta }_{2,1}\left(k\right)=1,\sum _{j=1}^2{\delta }_{3,j}\left(k\right)=1.$$

The state equation in the system (10) can be derived as

$$\begin{array}{ccc}\hfill x_1(k+1)& =& z_{1,1}\left(k\right)-z_{1,2}\left(k\right)+{\delta }_{1,2}\left(k\right)-z_{1,3}\left(k\right)-z_{1,4}\left(k\right)+z_{1,5}\left(k\right),\hfill \\ \hfill x_2(k+1)& =& x_1\left(k\right)-z_{2,1}\left(k\right),\hfill \\ \hfill x_3(k+1)& =& z_{3,1}\left(k\right)-z_{3,2}\left(k\right)+z_{3,3}\left(k\right)+z_{3,4}\left(k\right)-z_{3,5}\left(k\right),\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill z_{1,1}& =& {\delta }_{1,1}{x}_3,z_{1,2}={\delta }_{1,1}{x}_{3}u,z_{1,3}={\delta }_{1,2}u,z_{1,4}={\delta }_{1,2}{x}_3,z_{1,5}={\delta }_{1,2}{x}_{3}u,\hfill \\ \hfill z_{2,1}& =& x_1{x}_3,\hfill \\ \hfill z_{3,1}& =& {\delta }_{3,1}{x}_1,z_{3,2}={\delta }_{3,1}{x}_1{x}_2,z_{3,3}={\delta }_{3,2}{x}_2,z_{3,4}={\delta }_{3,2}u,z_{3,5}={\delta }_{3,2}{x}_{2}u.\hfill \end{array}$$

Using Lemma 2, these relations can be equivalently transformed into a set of linear inequalities. Then, $p=11$ holds. Thus, the linear inequality in the system (10) can be derived.

Using the system (10), consider rewriting Problem 2-C as an ILP problem. First, equality constraint (9) is embedded in the inequality constraint in the system (10). From the state equation in the system (10), we can obtain

$$x\left(k\right)=A^k{x}_0+\sum _{i=1}^k{A}^{i-1}Bv(k-i).$$

Using this expression, we can obtain

$$\overline{x}=\overline{A}{x}_0+\overline{B}\overline{v},$$

(11)

where $\overline{x}:={[x^{\top }\left(0\right)x^{\top }\left(1\right)\cdots x^{\top }\left(N\right)]}^{\top }$, $\overline{v}:={[v^{\top }\left(0\right)v^{\top }\left(1\right)\cdots v^{\top }(N-1)]}^{\top }$ and

$$\begin{array}{c}\hfill \overline{A}=\left[\begin{array}{c}I\\ A\\ A^2\\ ⋮\\ A^N\end{array}\right],\overline{B}=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ B& 0& \cdots & 0\\ AB& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ A^{N-1}B& \cdots & AB& B\end{array}\right].\end{array}$$

From the linear inequality in the system (10), we can obtain

$$\overline{C}\overline{x}+\overline{D}\overline{v}\le \overline{E},$$

(12)

where

$$\overline{C}=\left[\begin{array}{cccc}C& & & 0\\ & \ddots & & ⋮\\ & & C& 0\end{array}\right],\overline{D}=\left[\begin{array}{ccc}D& & \\ & \ddots & \\ & & D\end{array}\right],\overline{E}=\left[\begin{array}{c}E\\ ⋮\\ E\end{array}\right].$$

Next, consider the linear inequality ${\sum }_{k=0}^{N-1}L\delta \left(k\right)\ge \ln \rho$ in Problem C. From $\delta \left(k\right)=Wv\left(k\right)$, $W:=\left[0_{l\times m}{I}_l{0}_{l\times m}\right]$, this inequality can be rewritten as

$$-\overline{L}\overline{W}\overline{v}\le -\ln \rho ,$$

where $\overline{L}=[L\cdots L]$ and

$$\overline{W}=\left[\begin{array}{ccc}W& & \\ & \ddots & \\ & & W\end{array}\right].$$

Cost function (6) can be rewritten as

$$J=\overline{Q}\overline{x}+\overline{R}\overline{v},$$

(13)

where $\overline{Q}:=[Q\cdots Q{Q}_f]$ and $\overline{R}:=[R\cdots R]$. By substituting the state Equation (11) into the inequality (12) and the cost function (13), we can obtain the ILP problem. In [23], this result was obtained as the following theorem.

Theorem 3.

Problem 2-C is equivalent to the following ILP problem.

Problem 2-D:

$$\begin{array}{ccc}\hfill \mathrm{find}& & \overline{v}\in {\{0,1\}}^{(m+l+p)N}\hfill \\ \hfill min& & (\overline{R}+\overline{Q}\overline{B})\overline{v}+\overline{Q}\overline{A}{x}_0\hfill \\ \hfill \mathrm{subject}\mathrm{to}& & \left[\begin{array}{c}\overline{C}\overline{B}+\overline{C}\\ -\overline{L}\overline{W}\end{array}\right]\overline{v}\le \left[\begin{array}{c}\overline{E}-\overline{C}\overline{A}{x}_0\\ -\ln \rho \end{array}\right].\hfill \end{array}$$

Problem 2-D can be solved by using a suitable free/commercial solver.

Finally, we explain a solution method for Problem 2-B. In Problem 2-B, we assume that $u\left(k\right)$ (i.e., certain elements of $\overline{v}$) is given. Under this assumption, Problem 2-B can be transformed into an ILP problem by replacing “min” with “max” in Problem 2-D.

5. Discussion

In this section, from the viewpoint of numerical experiments, we compare two introduced methods [23,24] with other existing methods.

In many existing methods, transition probability matrices with $2^n\times 2^n$ in discrete-time Markov chains must be manipulated (see, e.g., [19,20,21,30,32,36]). In Example 2, the transition probability matrix with $u\left(0\right)=0$ can be obtained as

$$\left[\begin{array}{cccccccc}0.8& 0& 0& 0& 0.2& 0& 0& 0\\ 0.2& 0& 0& 0& 0.8& 0& 0& 0\\ 0.56& 0.24& 0& 0& 0.14& 0.06& 0& 0\\ 0.14& 0.06& 0& 0& 0.56& 0.24& 0& 0\\ 0& 0& 0.24& 0.56& 0& 0& 0.06& 0.14\\ 0.06& 0.14& 0& 0& 0.24& 0.56& 0& 0\\ 0& 0& 0.56& 0.24& 0& 0& 0.14& 0.06\\ 0.14& 0.06& 0& 0& 0.56& 0.24& 0& 0\end{array}\right].$$

See also the state transition diagram in Figure 1. In recent years, analysis and control methods for BNs and PBNs using the semi-tensor product approach have been widely studied (see, e.g., [34,36]). In this approach, matrices with $2^n\times 2^n$ must be manipulated. It is hard to manipulate matrices with such size. In MATLAB, matrices with $2^{16}\times 2^{16}$ cannot be generated. Although matrices with $2^{15}\times 2^{15}$ can be generated, computation of the product of matrices requires a long time. This is a weakness of many existing methods.

In two methods introduced in this review paper, matrices with such size are not manipulated. As a result, we can compute the optimal control problem. In [24], we presented a numerical experiment for the polynomial optimization approach. Then, we reported that for a PBN with 15 states and three control inputs, the computation times for solving Problem 1 in the cases of $N=2,3,4,5$ were 0.6 s ($N=2$), 19.6 s ($N=3$), 218.8 s ($N=4$), and 334.4 s ($N=5$), respectively. In this numerical experiments, we used SparsePOP [45] for solving a polynomial optimization problem. The computer environment is CPU: Intel Core i7 1.2 GHz, Memory: 4 GB, Windows 7 Professional 64 bit. In the integer programming approach, we presented in [23] a numerical experiment for a context-sensitive PBN with 15 states and three control inputs. Then, we reported that the computation time for solving Problem 2-A with $N=10$ was 18.96 s. In this numerical experiments, we used ILOG CPLEX 11.0 for solving an ILP problem. The computer environment is CPU: Intel Core 2 Duo CPU 3.0 GHz, Memory: 4 GB, Windows 7 Professional 64 bit. We remark that it is hard to solve the optimal control problem of a PBN with such size by using the other existing methods such as [20,21,34,36]. Thus, two methods introduced here enable us to compute the optimal control input for relatively large-scale PBNs.

6. Open Problems in Control Theory of Probabilistic Boolean Networks

We have explained our previously proposed methods, and we have compared these methods with other existing methods. In this section, we suggest open problems in control theory of PBNs. Here, we consider not only optimal control methods but also other control methods.

Although our previously proposed methods can handle large-scale PBNs, these methods are not scalable. Developing scalable control methods is important. Then, decomposition of networks is important (see, e.g., [50,51]). For BNs, we consider only one network structure (directed graph) such as Figure 1. On the other hand, for PBNs, we must consider a set of directed graphs. Hence, the problem of decomposing given networks in a PBN is more difficult than that in a BN. This is one of the challenging problem in control of PBNs. Moreover, the optimal control problems for PBNs obtained from decomposed graphs can be solved fast, but their solutions may not be an optimal solution for the optimal control problem of a given PBN. As a result, the optimality of the cost function is degraded. Now, some analysis of the approximation accuracy is required.

It is also important to consider what can be obtained from graph structure alone. In [22], controllability of BNs is characterized using a time sequence of graph structure. In [16], singleton attractors (fixed points) of BNs are characterized from graph structure only. In [52,53], a feedback vertex set of a graph is focused. A feedback vertex set of a graph is a set of vertices whose removal results an acyclic graph. For a given graph, several feedback vertex sets can be obtained, but in particular, a minimum feedback vertex set in which the number of vertices is minimum is important. A minimum feedback vertex set is also not given uniquely. It was shown in [52,53] that in the system given by a differential equation, a certain property of all states can be characterized by analyzing the states relating to a certain feedback vertex set. However, the related results for BNs and PBNs have not been obtained so far.

In our previously proposed methods explained in this paper, a controller must be implemented by the personal computer that can solve a polynomial optimization problem and an integer programming problem (of course, actuators for interacting genes are necessary, but this topic is out of scope of this review paper). We may consider the case where a controller must be described by a certain Boolean function. In such a case, the algebraic representation using the STP and the matrix-based representation [54,55] are useful. Details are one of the open problems.

Since BNs and PBNs are synchronous models, an extension to asynchronous models is important. Then, it is appropriate to utilize a Petri net representation. A Petri net is well known as a model of asynchronous behavior (see, e.g., [56] for further details), and is also used as a model of biological networks (see, e.g., [9,57]). In [58], the authors have proposed a Petri net representation for modeling asynchronous BN. In the case of PBNs, a stochastic Petri net may be used as a model of asynchronous PBNs, but further discussion is required.

Finally, the state and the control input in BNs and PBNs are binary. There is a possibility that biological behavior is not precisely modeled (see, e.g., [11]). Then, it may be appropriate that some of states/inputs are given by continuous-valued variables. A dynamical system with binary variables and continuous variables is called a hybrid dynamical system (see, e.g., [59]). Further research on modeling and control methods is needed.

Thus, in control theory of PBNs, there are many remaining open problems. Further development of many control methods is required.

7. Conclusions

The optimal control problem of PBNs is one of the significant problems in systems biology. This problem can be rewritten as the optimal control problem of discrete-time Markov chains, which has been extensively studied. However, the size of Markov chains exponentially grows with the state of PBNs. Hence, it is important to consider a solution method in which Markov chains are not computed. In this review paper, our previously proposed methods (i.e., the polynomial optimization approach and the integer programming approach) in control of PBNs were explained as computation methods that do not require Markov chains. In addition, we compared these methods with the other existing methods. Finally, we explained open problems in control theory of PBNs.