Maximum-Likelihood State Estimators in Probabilistic Boolean Control Networks

This study addresses state estimation problems for probabilistic Boolean control networks (PBCNs). Compared with deterministic Boolean networks, PBCNs have the stochastic switching in logical update functions in the state equation. Consequently, statistical analysis is required to estimate unavailable states, which induces an optimization problem called maximum-likelihood estimation. This article mainly focuses on two scenarios: 1) state estimation from partially measured state and 2) state estimation from output data, meaning observer design. The resulting optimization problems are solved using efficient algorithms based on dynamic programming. Concurrently, Dijkstra-type algorithms, which solve equivalent shortest path problems, are also proposed using best-first search. Furthermore, both the proposed algorithms derive novel observer design methods for PBCNs. The proposed algorithms are evaluated with practical estimation problems aiming to the sensor reduction and applied to gene regulatory networks of apoptosis and Lac operon.

the verification of the observability is a difficult task [8], a basic idea behind the state estimation of BCN is simple; existing reports have examined whether the unknown state can be deterministically calculated using available measured data. Indeed, the future states of the system can be calculated recursively by a given current state, and previous states are also calculated similarly. By exploiting a matrix multiplication technique called semitensor product (STP) [9], the explicit and simple equation that a current state and a future state satisfy is obtained as a linear equation whose coefficient is calculated using a structure matrix [3]. Therefore, the goal of the state estimation in BCNs results in finding a solution of some matrix equations, as illustrated in the motivating examples of Section III. Therefore, one observed state at any timing is enough to calculate the whole state trajectory. Using this deterministic feature, the control problems, including synchronization [10]- [13], pinning control [14]- [16], optimal control [17]- [19], and Lyapunov method [20], [21], were well addressed.
Unlike the estimation problem in BCNs, that of PBCNs is somewhat complicated due to the stochastic behaviors caused by the random switching in a state equation. In fact, existing reports expanding concepts of BCNs to PBCNs have required complicated and detailed discussions [18], [22]- [24]. Similar to these existing reports on PBCNs, deterministic calculations as summarized in the previous section are not applicable because of the stochasticity in PBCNs, and the adaptation of the deterministic techniques with BCNs approximation cannot work, as discussed in the motivating examples of Section III. Motivated by the previous researches as summarized above, the estimation problems for both missing observation data and observer design are addressed in this article. The authors believe that studies regarding the PBCN estimation are limited.
The estimation problems of PBCNs have two challenging issues: 1) a reasonable criterion justifying an obtained estimator is required, whereas that of BCNs is not necessary and 2) an efficient algorithm to find an estimator satisfying the given criterion must be designed. For Point 1), this study employs the maximum-likelihood criterion, which results in an estimator that maximizes the likelihood of an estimated state trajectory. Using the STP technique [25] and its extension to PBCNs (see [26]), the formulation itself is also naturally given as matrix expressions with the STP. However, as for Point 2), because the estimation problem results in an optimization problem, this article proposes multiple algorithms to solve the problem efficiently. The proposed formulation can be converted into a multistage decision problem and provides an efficient optimization algorithm.
In addition, some variations of the maximum-likelihood estimation, including the estimation from the output signals, are addressed. In practical systems, especially biological systems, measuring all state variables requires many sensors; therefore, the number of measured variables is preferred to be reduced [27]. In Examples 4 and 5 of Section VI, the estimation problem using outputs of a system can address a situation with partially measured states, as illustrated in application to two well-known gene regulatory networks, an apoptosis network and Lac operon.
The main contributions of this article can be summarized as follows.
1) A new formulation of the estimation problem in PBCNs based on the maximum-likelihood estimation is introduced. Drawbacks of the conventional approach using deterministic BCN approximation are also highlighted by illustrating a simple example. 2) An equivalent reformulation of the estimation problem as an optimal control problem is deduced, and efficient optimization algorithms for the maximum-likelihood estimation are derived in the framework of the dynamic programming (DP) and the shortest path approach. The computation efficiency is examined using practical biological examples in detail.
As representative examples, this study addresses the gene regulatory networks of the apoptosis and the lac operon. 3) By expanding the state estimation problem summarized above, further results for the observer design just using the output of PBCNs are obtained. Algorithms for the resulting optimization problem are also proposed. By applying the proposed frameworks to the two biological systems-the apoptosis network and the lac operon network-the feasibility of the sensor reduction is validated.

II. NOTATIONS AND PRELIMINARIES
The notations used in this article are listed as follows. identified with vectors δ 1 2 and δ 2 2 , respectively. 5) The STP [9] of two matrices A ∈ M m×n and B ∈ M p×q is defined by In the equation above, ⊗ is the Kronecker product, and s = lcm(n, p) is the least common multiple of the integers n and p. In this article, the notation " " in an expression A B is omitted and simply written as AB.

III. PROBLEM FORMULATION
In PBCN, the ith state of a system is given by the following equation: Switching probabilities c (i) j , j = 1, . . . , n (i) , which represent an occurrence probability of each jth Boolean function f First, a theorem that provides the transition matrix of the PBCN is shown.
Theorem 1 [26]: Consider the PBCN in (1). Define a matrix K ∈ M n all ×n x as Obtain the structure matrix [25] of the logical function f (i) where n u +n x is the swapping matrix proposed in [25]. Then, the state variable X(k) satisfies the following equation: where E[X(k)] is the expectation of the state X(k).
In this article, the matrix T is called the structure matrix of the PBCN (1). The following theorem calculates transition probabilities.
Theorem 2 [28]: The transition probability of the PBCN (1) from the kth step to the k + 1th step is given by the following equation: where T is the structure matrix of the PBCN (1). Consider a situation that an open-loop input sequence u(0), . . . , u(N−1) is applied to the system (1), and the terminal state X(N) = x(N) is observed, whereas the other states x(k), k = 0, . . . , N − 1, which includes the initial state x(0), are missing data, that is, not observed. Fig. 1 shows an example with the length N = 3. In this article, a realization of a random variable X(k) is denoted by a small letter x(k), and  (14). [For the sake of clarity, the states δ 1 2 and δ 2 2 are colored in yellow and gray, respectively. Note that the states x(0), x (1), and x(2) in circles with dashed lines are unknown variables to be estimated, whereas the other variables u(0), u(1), u (2), and x(3) are given in advance and depicted by solid lines.] its generation probability Pr(X(k) = x(k)) is often denoted by Pr(x(k)) to simplify the notation. The generation probability of the states x(0), . . . , x(N) under the input u(k), . . . , u(N −1) is expressed by the product of two terms which depends on x(0), . . . , x(N − 1) and x(N), respectively, as follows: In (4), the second equality uses the causality of the first term, meaning that the states x(0), . . . , x(N − 1) are independent of u(N − 1), and the latter part is simplified based on the Markov property of the system, that is, x(N) depends only on x(N − 1) and u(N − 1). Use (4) similarly for k = 0, . . . , N − 1, and obtain the following form: Pr(x(k + 1)|x(k), u(k)).
In the equation above, the distribution of X(0), which is denoted as Pr(x(0)), can be expanded to Therefore, (5) is expressed as the likelihood function of x(0), . . . , x(N − 1), and the distribution of X(0), that is, . . , 2 n x as follows: Pr(x(k + 1)|x(k), u(k)). (7) In the equation above, by recalling that u(0), . . . , u(N − 1) and x(N) are given in advance with a fixed control period N, the arguments of the function L result in x(0), . . . , x(N − 1). Consequently, the maximum-likelihood estimator of L 0 is given as an optimal solution of the following maximization problem: The problem above has differences from conventional filtering problems with innovation, that is, the target system is not given as a differential or difference equation but a PBCN with Boolean variables with its randomness on the switching probabilities.
The optimal d i , i = 1, . . . , n x , are given by the following proposition.
Proposition 1: Consider the PBCNs (1) and the associated two optimization problems-Problem (8) and the following optimization problem: Then, the following statements hold. 1) Parts of an optimal solution corresponding to the arguments x(0), . . . , x(N − 1) in Problem (8) are given an optimal solution of Problem (9) denoted by x * (0), . . . , x * (N − 1). 2) The remaining part of the optimal solution corresponding to the arguments d 1 , . . . , d 2 nx in Problem (8) are given by for each i = 1, . . . , 2 n x . Proof: By recalling the expression of L 0 in (6) and letting L(x(0), . . . , x(N −1)) = N−1 k=0 Pr(x(k +1)|x(k), u(k)) in (7) to simplify the notation, the following reformulation is derived: (11) subject to 0 ≤ d i ≤ 1 for each i = 1, . . . , 2 n x with an equality constraint 2 nx i=1 d i = 1 and x(k) ∈ 2 nx for each k = 0, . . . , N − 1. Because x(0) ∈ 2 nx , for an index i, the following calculation holds: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. and its equality holds if d i = 1 for x(0) = δ i 2 nx and the other d j , j = i are all zeros. In (11) and its equality is satisfied by x * (0), . . . , x * (N − 1) of item 1) and d * 1 , . . . , d * 2 nx of item 2) in the statements of the theorem, and these compose an optimal solution of Problem (8).
Based on the proposition above, the optimization problem (9) is considered hereafter, although the main result of the study can be readily applied to a case that the initial distribution Pr(X(0) = δ i 2 nx ), i = 1, . . . , 2 n x , is given in advance (Remark 2). For the sake of comparison with previous reports, the next section addresses numerical examples revealing the drawback of conventional approaches based on the approximation as deterministic BCNs. Furthermore, a motivating example is examined. Consider the following simple PBCN with one state and input: An applied control sequence is assumed to be u(0) = F, u(1) = T, and u(2) = F, and the terminal state is observed as x(3) = δ 1 2 ( Fig. 1). Here, the estimation of the states x(0), x (1), and x(2) are considered, and two cases with different values of the selection probability c are examined; the estimation problems for c = 0.7, which is called case 1, and for c = 0.3, called case 2, are explored. In both cases, as a representative of existing approaches, a BCN approximation is also applied. Although the two approaches, meaning the BCN approximation and the proposed formulation, for case 1 result in the same estimation results, case 2 indicates that the BCN approximation cannot derive a solution. That is, the latter case 2 implies the necessity of the proposed framework, which exploits the maximum-likelihood estimation.
Consider case 1, meaning the case of c = 0.7, where the BCN approximation works correctly. The switching probabilities are given by Pr [f = f 1 ] = 0.7 and Pr [f = f 2 ] = 0.3. Since the system is mainly driven by the first logic function f 1 , a way to estimate the states is that the system is approximated by a BCN with a logical system Using the deterministic BCN (15), the following recursive calculations hold: Furthermore, a maximum-likelihood estimator, that is, a solution of Problem (9), is obtained by listing all the paths as depicted in Fig. 2. Fig. 2 shows that the two paths , obtained by the BCN approximation, are indeed optimal regarding Problem (9) with the optimal value 0.49. Fortunately, in this example, the BCN approximation derives the optimal path of Problem (9). However, the next example reveals the drawback of this approach.
Second, consider case 2, meaning c = 0.7, where the BCN approximation results in an incorrect result. The switching probabilities are given as Pr [f = f 1 ] = 0.3 and Pr [f = f 2 ] = 0.7. Similar to case 1, the system is approximated by the following deterministic BCN: If the state x(k + 1) and the input u(k) are given, x(k) is determined by the equation above. However, although equations x(3) = δ 1 2 and x(2) = δ 2 2 can be obtained, x(1) cannot be calculated as a solution of an equality x(2) = δ 2 [1, 1, 2, 1]u(1)x(1), meaning that there is no state x(1) ∈ 2 satisfying δ 2 2 = δ 2 [1,1]x(1). However, using the proposed formulation in (9) in the framework of maximum-likelihood estimation for probabilistic BCNs, a reasonable estimation as a solution of an optimization problem is available. To illustrate the problem simply, the structure matrices of the logical functions f 1 and f 2 in (14) are calculated as M 1 = δ 2 [1, 1, 1, 2] and M 2 = δ 2 [1, 1, 2, 1], respectively, and Theorem 1 calculates the structure matrix of the PBCN in (16) as and the transition probabilities are obtained by Theorem 2. Furthermore, a maximum-likelihood estimator, that is, a solution of Problem (9) is obtained by listing all the paths (Fig. 2). Fig. 2 shows that there are two optimal paths which achieve the optimal value 0.21.

IV. OPTIMIZATION ALGORITHMS FOR MAXIMUM-LIKELIHOOD ESTIMATION
The listing approach (Fig. 2) requires much computation cost. However, an efficient algorithm is available using the following theorem.
Theorem 3: Consider the maximum-likelihood estimation problem of the PBCN in (1) with unavailable states where the matrix T is the structure matrix of the PBCN in (1). Then, the optimal solution x * (k), k = 0, . . . , N − 1, is given by Proof: Take the logarithm of the object function L in (7) and obtain where the equation in Theorem 2, meaning Pr( , is used in the second equality. The following state transformation is applied: In the equation above, the variables ξ and v are a state variable and a control input, respectively, and the inversed time index s = N − 1 − k for the original index k = 0, . . . , N − 1 is introduced. In the first equation of (19), by replacing s with s + 1, the following calculation holds: Backward Calculation on s 4: for s = N − 1, . . . , 0 do 5: meaning the state equation of the state variable ξ(s) and the control input v(s) at the sth time period is obtained. Similarly, by replacing a time variable k in (18) and it is indeed equal to that of (17).
The application of the DP algorithm [29] to an optimal control problem has been addressed in existing papers [30], [31]. Based on the reports, DP for (17) is given by DP-MLSE (Algorithm 1).
Although the details of the DP algorithm itself cannot be thoroughly explained due to space limitation, the interpretations of variables in Algorithm 1 are summarized as follows. An interested reader can refer to existing reports [30], [31]. 1) μ s,i is an optimal policy, meaning an optimal control input v(s) at the sth time step, of a state variable ξ(s) = δ i 2 nx . Please note μ s,i is not a control input belonging to 2 nx but an index of a control input with a range μ s,i ∈ {1, . . . , 2 n x }.
2) A vector variable V s is a so-called value function, whose element V s,i is the minimum cost with an initial state ξ(s) = δ i 2 nx for accumulated stage cost from the sth time step to the Nth step. Consider the 6th and 7th lines of Algorithm 1 and examine the integrity of Algorithm 1 by referring to existing reports on the DP algorithm (see [30]). The first part of the right side can be reformulated as , that is, the term means the stage cost at the sth step. Furthermore Indeed, what the basic idea of the DP, well known as Bellman's principle of optimality, claims are summarized in the equations above; the cost starting from the state ξ(s) = δ i 2 nx at the sth step is the sum of the current stage cost in Section IV and the optimal cost starting from the s + 1th step with the state ξ(s + 1) = v(s).
Remark 1: In DP-MLSE (Algorithm 1), if the argument of log is given by a matrix, the logarithm of each element is calculated. When there is a zero element in the transition matrix Tu (N − 1 − s), a calculation rule log(0) = −∞ is applied formally, and a multiplication −∞ × 0 = 0 in the sixth and the seventh lines is also assumed. Based on these rules, an objective function value is correctly obtained.
The computation complexity is determined by iterations from the fourth line to the seventh line of DP-MLSE (Algorithm 1). The sixth and seventh lines require 2 n x evaluations of the right side for each iteration, and the numbers of iterations for the fourth and fifth lines are given by N and 2 n x , respectively. Therefore, multiplying them, the whole computation order is obtained as O(N · 2 2n x ).
Example 1: In case 1, the whole paths and corresponding likelihoods are calculated, as shown in Fig. 2. In this example, Algorithm 1 is applied, and a feedback law μ and a value function V are deduced (Table I). Using the obtained feedback μ, the 9th to 11th lines of Algorithm 1 provide two optimal paths (x * (0), x * (1), x * (2), x * (3)) = (δ 1 2 , δ 2 2 , δ 1 2 , δ 1 2 ) and (δ 2 2 , δ 1 2 , δ 1 2 , δ 1 2 ). Note that since the value of V is the logarithm of the likelihood function, a calculation exp(V 0,1 ) = exp(−1.5606) = 0.21 deduces the maximum-likelihood value (Fig. 2).  (6) and examining the state transform (19) in the same way as Theorem 3, the objective function is given as The expression above comprises the objective function in (17) Take the logarithm of both sides of the equation above, and separate a term depending on x(0), . . . , x(N −1) and the others as follows: where the function L is defined in (7). Since the first term of the right side depends only on one control variable x(0), which is equal to the terminal state of the transformed state ξ(N), the following modification of the value function V N similar to Remark 2 is performed: and Algorithm 1 is applicable. In the equation above, because the first term is independent of x(0), . . . , x(N−1), meaning it is a constant, only x(−1) and u(−1) affect the optimal solution. That is, the maximum-likelihood does not depend on the length of the data d.

V. SHORTEST PATH-BASED ESTIMATION ALGORITHM
The previous section proposed the optimal control problem (17) to examine DP-based algorithm (DP-MLSE, Algorithm 1). Since the stage cost function is nonnegative, an equivalent shortest path problem, which is illustrated in Proposition 3, is deduced. Before the proposition, the following definitions are introduced. if d v < +∞ then 13: if d v > d then 14: p

15:
HeapUpdate(Q, [v, d ]). 16: else 17: respectively. The numbers of elements are given by |V| = 2 n x · N + 1 and |E| = 2 n x (2 n x (N − 1) + 1). }, respectively. The proposition above implies that the Dijkstra algorithm for shortest path problems is applicable. Note that compared with conventional shortest path problems, SP-MLSE (Algorithm 2) can terminate if the set Q in SP-MLSE includes }, which is given as the set of the terminal vertices.
It is well known that the Dijkstra algorithm can exploit the priority queue implemented by the heap structure (see [32, Sec. 6.2] for further details). Simply speaking, the heap structure is the tree structure that the values of children exceed those of their parents; therefore, the minimum value appears as the root of the tree. Remark 4: The comparisons between DP-MLSE and SP-MLSE are conducted here. Although the computation complexity of SP-MLSE exceeds that of DP-MLSE due to the operations related to the priority queue and heap, the improvement of practical computation time can be expected since SP-MLSE does not necessarily calculate every cost of the edges. Note that DP-MLSE is a successive realization of bruteforce search, whereas SP-MLSE attempts effective search as best-first search. For example, refer to the numerical example (Example 2) and its result plot (Fig. 5). Table II presents the overview of the proposed algorithms. Although the numerical experiments of this article imply the advantage of SP-MLSE regarding the practical computation time, another advantage of DP-MLSE is the simplicity of the implementation. By exploiting the matrix expression given by the semitensor technique, the 6th and 7th lines of Algorithm 1 are given as simple matrix calculations, which can be implemented in a fast vectorized command for multicore processors, and the parallelization of the for loop for an index i in the 5th line is also applicable.
Like conventional optimal control problems and related studies (see [33]), the computation complexity is exponential regarding the dimension of the state n x , which is the curse of dimensionality. Therefore, the proposed framework is limited to the estimation problem whose dimension of the state is  not high. At least, as illustrated in the Lac operon network of Example 3, some classes of high-order systems can be handled by the proposed algorithms.
Remark 5: As illustrated in the pseudocodes of Algorithms 1 and 2, contrast to general optimization algorithms, there is no tuning parameter because the proposed algorithms search an optimal solution by listing feasible solutions efficiently.

Example 3 (Large-Scaled Problem: Lac Operon Network):
In this numerical example, the Boolean model of the lac operon derived in [36] is addressed. For details, please refer to [36] and references therein. The variables are given as follows: M(lac mRNA); B(β-galactosidase, LacZ); R(repressor protein, LacI); A(allolactose); L(lactose); P(transport protein, LacY, "lac permease"); and C(cAMP-CAP protein complex). A variable with a subscript "e" or "m" represents the extracellular and the least medium concentrations, respectively. The other variables represent the intracellular concentrations. Note that a Boolean variable for a biological quantity is True, which corresponds to 1 and δ 1 2 , if its concentration is high, and the Boolean variable is False, which is equivalent to 0 and δ 2 2 , otherwise. By summing up the notations above, the state variable is given by (x 1 , . . . , x 10 ) = (M, B, R, A, L, P, C, R m , A m , L m ), and the control input is (u 1 , u 2 , u 3 ) = (L e , L em , G e ). Boolean update functions for the states are given as follows: A graphical interpretation of the system is illustrated in Fig. 6.
Here, assume the measurement of the concentrations of the inputs u 1 = L e and u 2 = L em has stochastic uncertainties. In other words, flippings on u 1 and u 2 occur with probabilities In this numerical example, for the sake of the comparison of the practical computation time, the estimation problem (9) is computed by DP-MLSE (Algorithm 1) and SP-MLSE (Algorithm 2). Increasing an estimation period length N, open-loop input sequences u 1 (0), u 2 (0), u 3 (0), . . . , u 1 (N − 1), u 2 (N − 1), u 3 (N − 1) and the deterministic terminal state x(N) are randomly given, and each problem is solved. The flipping probabilities are assumed to be c 1 = c 2 = 0.1. The computations were performed using a laptop computer with an AMD Ryzen 3 PRO 3300U 2.1-GHz processor and 8.00-GB RAM in Python 3.7.9 environment. The first column of Fig. 7 depicts the number of computed edges. Obviously, that of DP-MLSE grows linearly, whereas that of SP-MLSE can reduce the computation of the edges. In addition, from the second column of Fig. 7, the computation time can be improved by SP-MLSE compared with DP-MLSE. Fig. 8. Histogram of computation time with randomized numerical experiment. The labels naive and heap stand for SP-MLSE (Algorithm 2) with and without the binary heap for the priority queue. Note that the fourth row depicts all the histograms of three different methods concurrently for the comparison, whereas the other first, second, and third rows depict that of each method separately.
Furthermore, for the sake of considering the stochastic uncertainty of the flipping probabilities c 1 and c 2 , a randomized numerical simulation where c 1 and c 2 are randomly sampled from independent uniform distribution, meaning c 1 and c 2 are uniform random numbers, are performed. Fig. 8, which is a histogram depicted by 500 trials with a problem instance N = 10, indicates that the computation time of the naive implementation of SP-MLSE (Algorithm 2) outperforms that of DP-MLSE (Algorithm 1).
Compared with DP-MLSE (Algorithm 1), the improvement of the computation time by SP-MLSE (Algorithm 2) is based on the reduction of the calculated edge number, depicted in the first row of Fig. 7. If the computation times are assumed to be normally distributed, conventional Welch's t-test static [37] between the computation time of an arbitrary pair of the algorithms is reduced to almost zero. Therefore, the improvement of the computation time is statistically significant. As illustrated in Proposition 4, the theoretical computation complexity of DP-MLSE surpasses that of SP-MLSE; however, SP-MLSE can reduce the number of calculated edges by the best-first search scheme.

VI. FURTHER RESULTS FOR STATE OBSERVER DESIGN
Consider the state estimation of the PBCN in (1) whose state X(k) cannot be observed directly, but outputs Y l (k) ∈ 2 , l = 1, . . . , n y , are obtained through output functions g (l) : 2 nx → 2 , l = 1, . . . , n y Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. i = 1, . . . , n x Y l (k) = g (l) X 1 (k), . . . , X n x (k) , l = 1, . . . , n y .
Similarly, the output logic functions g (l) , l = 1, . . . , n y , are represented by a structure matrix C, which derives an equation Y(k) = CX(k) with the state variable X(k) and a combined output variable Y(k) = Y 1 (k) . . . Y n y (k) ∈ 2 ny . Consider the observer design for the states x(0), . . . , x(N) and the initial distribution of X(0), that is, under the observations y(0), . . . , y(N).
Eventually, the following maximum-likelihood estimation problem derives the optimal estimation of x(0), . . . , x(N): The following definition is given prior to showing an optimization algorithm.
Definition 2: Assume the structure matrix of a logic function f : 2 nx → 2 ny is given by M. That is, if an equation y = Mx holds, then a set-valued function M + : 2 ny → 2 2 nx is defined as The following theorem presents an equivalent optimal control problem.
If the objective function takes −∞ value, it is not optimal, meaning 1 ≤ j ≤ 2 n x log Pr y(k)|δ is the set of feasible indices j, and the constraints x(k) ∈ C + [y(k)], k = 0, . . . , N, are required. The claim of the theorem is obtained by considering the system of the state z(k) and the input v(k) with the system equation z(k + 1) = v(k), and paying attention to the conditions x(k) ∈ C + [y(k)] and an equation Pr(x(k + 1)|x(k), u(k)) = x T (k + 1)Tu(k)x(k), which is followed by Theorem 2. Now, an optimization algorithm for Problem (24) is given as follows.
Remark 6: Consider the computation complexity of DP-MLSE-O (Algorithm 3). The number of evaluations of the right side in the ninth and the tenth line of Algorithm 3 is |I s |, whereas the sixth and seventh lines of DP-MLSE (Algorithm 1) require evaluations for all integers 1 ≤ j ≤ 2 n x . Since |I s | ≤ 2 n x , the upper bound of the computation complexity is given by O(N · 2 2n x ) like DP-MLSE (Algorithm 1).
Remark 7: SP-MLSE (Algorithm 2) can also be extended to solve Problem (24). Indeed, the seventh line of SP-MLSE (Algorithm 2), which calculates costs of edges starting from Here, this modified algorithm of SP-MLSE (Algorithm 2) is called SP-MLSE-O, and the detailed algorithm is omitted here because the modification is straightforward. SP-MLSE-O can improve the computation load related to heap operations by reducing the number of edges |E|, whereas SP-MLSE (Algorithm 2) searches an optimal path from all the edges.
Example 4 (Sensor Reduction in Apoptosis Network): The apoptosis network model in Example 2 is reconsidered. Assume the states X 2 and X 3 cannot be directly observed, and an output equation is given by Y(k) = X 1 (k), that is, the maximum-likelihood estimation of X 2 and X 3 from the observation of X 1 is considered. The target situation is illustrated in Fig. 9. An estimation horizon is set as N = 3, and the control and observation sequences are illustrated in Fig. 10. The optimal solution is given by The structure matrix of the output function g(x) = x 1 = IAP is C = δ 2 [1, 1, 1, 1, 2, 2, 2, 2], and the inverse is First, a computation result with the estimation period N = 10 is illustrated in Fig. 11. To simplify the figure, the Boolean values {T(True) = 1 δ 1 2 , F(False) = 0 δ 2 2 } are colored by blue and white, respectively. For given input and output sets illustrated in the first and second blocks, the estimated states are obtained as depicted in the third block. Fig. 11 implies plausible estimated states are available from the partially measured states (x 1 , x 2 , x 3 ) = (M, B, R), which are shown in the second block. Like Example 3, the stochastic uncertainty of the flipping probabilities c 1 and c 2 are examined here. A randomized numerical simulation where c 1 and c 2 are uniform random numbers are carried out. Fig. 12 is a histogram depicted by 500 trials with a problem instance N = 10. Fig. 12 shows that the computation time of the naive implementation of SP-MLSE-O outperforms that of the DP implementation.
Second, the comparisons of the computation time are illustrated in Fig. 13. The computation environment is the same as Example 3.

VII. CONCLUSION
In this study, the state estimation problems for PBCNs and their further application to observer design were addressed in the framework of maximum-likelihood estimation. To solve the optimization problems associated with the resulting likelihood functions, efficient algorithms were proposed by utilizing the DP and the shortest path approach. Furthermore, as typical examples of gene regulatory networks, the apoptosis network and the lac operon network were examined, and detailed comparisons of the theoretical computation complexity and the practical computation time were carried out. The numerical examples of the gene regulatory networks indicated that the proposed framework can provide plausible estimation in a reasonable computation time.

VIII. FUTURE RECOMMENDATIONS
This article explored the maximum-likelihood estimation problems and related efficient algorithms. Apart from some variants of the estimation problems, such as Remarks 2 and 3, other settings, inducing outliers and partial missing data, are also possible in practical situation. Further extensions for various estimation settings, including Bayes inference and resulting maximum a posteriori (MAP) estimation problems, can be considered. It is expected that the similar DP and shortest path-based approaches are applicable. Furthermore, from the viewpoint of the model structure, the estimation problems of other stochastic models for the BCNs, including BCNs with Markovian jump [38], are further examined.
For practical applications, although this article discussed biological examples, the application to another practical systems, including cyber-physical systems [39], can be further examined. If disturbances to a target system, such as cyberattacks, are modeled as a stochastic switching of the state and the system is simplified as BCNs, the considered system results in PBCNs and derives the same estimation problems discussed in this article. A more practical model with stochastic disturbance and resulting filter design problem are also open issues to be challenged.