Particle Filters for Partially-Observed Boolean Dynamical Systems

TL;DR

This paper develops particle filter-based algorithms for estimating states and parameters in partially-observed Boolean dynamical systems, enabling scalable analysis of complex Boolean processes like gene regulatory networks from noisy data.

Contribution

It introduces approximate MMSE filtering and smoothing algorithms using auxiliary particle filters for large POBDS, along with a specialized EM-based smoother for parameter estimation.

Findings

Algorithms successfully applied to gene regulatory network data

Improved computational efficiency over exact methods

Effective state and parameter estimation demonstrated

Abstract

Partially-observed Boolean dynamical systems (POBDS) are a general class of nonlinear models with application in estimation and control of Boolean processes based on noisy and incomplete measurements. The optimal minimum mean square error (MMSE) algorithms for POBDS state estimation, namely, the Boolean Kalman filter (BKF) and Boolean Kalman smoother (BKS), are intractable in the case of large systems, due to computational and memory requirements. To address this, we propose approximate MMSE filtering and smoothing algorithms based on the auxiliary particle filter (APF) method from sequential Monte-Carlo theory. These algorithms are used jointly with maximum-likelihood (ML) methods for simultaneous state and parameter estimation in POBDS models. In the presence of continuous parameters, ML estimation is performed using the expectation-maximization (EM) algorithm; we develop for this purpose a special smoother which reduces the computational complexity of the EM algorithm. The resulting particle-based adaptive filter is applied to a POBDS model of Boolean gene regulatory networks observed through noisy RNA-Seq time series data, and performance is assessed through a series of numerical experiments using the well-known cell cycle gene regulatory model.

Tables

Table 1. Table 1: Derivatives of Q ​ ( θ , θ ( n ) ) 𝑄 𝜃 superscript 𝜃 𝑛 Q(\theta,\theta^{(n)}) with respect to different parameters needed in Numerical Experiment.

Parameter	derivation
$p$	$\frac{1}{N}\sum _{s=1}^T\sum _{i=1}^N\left(\frac{{\Vert {\widetilde{𝐱}}_{s,i}\oplus 𝐟\left({\widetilde{𝐱}}_{s-1,i}\right)\Vert }_1-dp}{p\left(1-p\right)}\right)$
$s$	$\frac{1}{Ns}\sum _{s=1}^T\sum _{i=1}^N\sum _{j=1}^d\left(\frac{{\varphi }_j\left({𝐘}_s\left(j\right)-s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)\right)}{{\varphi }_j+s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)}\right)$
$\mu$	$\frac{1}{N}\sum _{s=1}^T\sum _{i=1}^N\sum _{j=1}^d\left(\frac{{\varphi }_j\left({𝐘}_s\left(j\right)-s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)\right)}{{\varphi }_j+s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)}\right)$
${\delta }_j$	$\frac{1}{N}\sum _{s=1}^T\sum _{i=1}^N\left(\frac{{\widetilde{𝐱}}_{s,i}\left(j\right){\varphi }_j\left({𝐘}_s\left(j\right)-s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)\right)}{{\varphi }_j+s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)}\right)$
${\varphi }_j$	$\frac{1}{N}\sum _{s=1}^T\sum _{i=1}^N\left(\frac{{\Gamma }^{\prime }\left({𝐘}_s\left(j\right)+{\varphi }_j\right)}{\Gamma \left({𝐘}_s\left(j\right)\right)+{\varphi }_j)}-\frac{{\Gamma }^{\prime }\left({\varphi }_j\right)}{\Gamma \left({\varphi }_j\right)}+\frac{s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)-{𝐘}_s\left(j\right){\varphi }_j}{{\varphi }_j\left(s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)+{\varphi }_j\right)}+\log \frac{{\varphi }_j}{s\exp \left(\mu +{\delta }_j{\widetilde{𝐱}}_{s,i}\left(j\right)\right)+{\varphi }_j}\right)$

Table 2. Table 2: Parameter values for numerical experiments using the Cell-Cycle gene regulatory network.

Parameter	Value
Length of time series $T$	50, 100
Number of genes $d$	10
Initial distribution $P\left({𝐗}_0={𝐱}^i\right),i=1:2^{10}$	$1/2^{10}$
Number of particles $N$	200, 1000, 5000
Bias $b_i,i=1,\dots ,10$	-1/2
Transition noise intensity $p$	0.01, 0.05
Sequencing depth $s$	1.02 (1K-50K reads)
Baseline expression $\mu$	0.1
Differential expression ${\delta }_i,i=1,\dots ,10$	2
Inverse dispersion ${\varphi }_i,i=1,\dots ,10$	1, 5
APF-CPMLA-BKS stopping criterion $\epsilon$	${10}^{-4}$

Table 3. Table 3: Experiment 1: Average rates of correct state estimation over 1000 independent runs for a time series with length 100.

$p$	$\varphi$	$N$	BKF	APF-BKF	BKS	APF-BKS
0.01	5	200	93.9	85.4	96.6	88.1
		1000		92.1		95.0
		5000		93.2		95.7
	1	200	83.8	74.6	90.7	80.4
		1000		80.6		88.3
		5000		82.1		89.8
0.05	5	200	82.9	75.0	93.4	82.3
		1000		80.3		91.3
		5000		81.9		92.6
	1	200	58.5	50.1	70.8	62.3
		1000		55.1		68.2
		5000		56.9		69.9

Table 4. Table 4: Experiment 2 results. Average accuracy rates for estimation of the gene interaction parameter a 42 subscript 𝑎 42 a_{42} .

$n$	$p$	$N$	High Disp.	Low Disp.
			$\varphi =1$	$\varphi =5$
30	0.01	200	0.45	0.62
		1000	0.58	0.71
		5000	0.59	0.73
	0.05	200	0.52	0.71
		1000	0.64	0.74
		5000	0.67	0.75
60	0.01	200	0.82	0.86
		1000	0.87	0.93
		5000	0.89	0.93
	0.05	200	0.86	0.91
		1000	0.89	0.95
		5000	0.92	0.96

Equations

$${\mathbf{X}}_k$$

$${\mathbf{Y}}_k$$

$$\mathrm{MSE}({\hat{\mathbf{X}}}_{k|r}\mid {\mathbf{Y}}_{1\text{mathchar}58r})=E\left[||{\hat{\mathbf{X}}}_{k|r}-{\mathbf{X}}_k|{|}^2\mid {\mathbf{Y}}_{1\text{mathchar}58r}\right].$$

$${\hat{\mathbf{X}}}_{k|r}^{\mathrm{MS}}=\arg \underset{{\hat{\mathbf{X}}}_{k|r}\in X_{k|r}}{\min }\mathrm{MSE}({\hat{\mathbf{X}}}_{k|r}\mid {\mathbf{Y}}_{1\text{mathchar}58r}),$$

$${\hat{\mathbf{X}}}_{k|r}^{\mathrm{MS}}=\overline{E\left[{\mathbf{X}}_k\mid {\mathbf{Y}}_{1\text{mathchar}58r}\right]},$$

$$\mathrm{MSE}({\hat{\mathbf{X}}}_{k|r}^{\mathrm{MS}}\mid {\mathbf{Y}}_{1\text{mathchar}58r})$$

$$\displaystyle\>=\,\big{|}\big{|}\min\big{\{}E\left[{\bf X}_{k}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}k}\right],E\left[{\bf X}_{k}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}k}\right]^{c}\!\big{\}}\,\big{|}\big{|}_{1},$$

$${\mathbf{\Pi }}_{k|r}(i)=P\left({\mathbf{X}}_k={\mathbf{x}}^i\mid {\mathbf{Y}}_{1\text{mathchar}58r}\right),i=1,\dots ,2^d,$$

$${\hat{\mathbf{X}}}_{k|r}^{\mathrm{MS}}=\overline{A{\mathbf{\Pi }}_{k|r}},$$

$$\mathrm{MSE}({\hat{\mathbf{X}}}_{k|r}^{\mathrm{MS}}\mid {\mathbf{Y}}_{1\text{mathchar}58r})=||\min \{A{\mathbf{\Pi }}_{k|r},(A{\mathbf{\Pi }}_{k|r}{)}^c\}|{|}_1.$$

$$(M_k{)}_{ij}=P({\mathbf{X}}_k={\mathbf{x}}^i\mid {\mathbf{X}}_{k-1}={\mathbf{x}}^j)$$

$$=P\left({\mathbf{n}}_k=\mathbf{f}({\mathbf{x}}^j)\oplus {\mathbf{x}}^i\right),i,j=1,\dots ,2^d.$$

$${\left(T_k({\mathbf{Y}}_k)\right)}_{ii}=p\left({\mathbf{Y}}_k\mid {\mathbf{X}}_k={\mathbf{x}}^i\right),i=1,\dots ,2^d,$$

$${\mathbf{\Delta }}_{k|s}(i)=p\left({\mathbf{Y}}_{s+1},\dots ,{\mathbf{Y}}_T\mid {\mathbf{X}}_k={\mathbf{x}}^i\right),i=1,\dots ,2^d,$$

$${\mathbf{\Pi }}_{k|T}=\frac{{\mathbf{\Pi }}_{k|k-1}\bullet {\mathbf{\Delta }}_{k|k-1}}{||{\mathbf{\Pi }}_{k|k-1}\bullet {\mathbf{\Delta }}_{k|k-1}|{|}_1},$$

$$P({\mathbf{X}}_k,{\zeta }_k\mid {\mathbf{Y}}_k)\propto p({\mathbf{Y}}_k\mid {\mathbf{X}}_k)P({\mathbf{X}}_k\mid {\mathbf{x}}_{k-1,{\zeta }_k})W_{k-1,{\zeta }_k},$$

$$P({\mathbf{X}}_k,{\zeta }_k\mid {\mathbf{Y}}_k)\propto p({\mathbf{Y}}_k\mid {\mu }_{k,{\zeta }_k})P({\mathbf{X}}_k\mid {\mathbf{x}}_{k-1,{\zeta }_k})W_{k-1,{\zeta }_k},$$

$${\mu }_{k,i}$$

$$=\text{Mode}[\mathbf{f}({\mathbf{x}}_{k-1,i})\oplus {\mathbf{n}}_k]=\mathbf{f}({\mathbf{x}}_{k-1,i}),$$

$$V_{k,i}=p({\mathbf{Y}}_k\mid {\mu }_{k,i})W_{k-1,i},$$

$${\mathbf{x}}_{k,i}={\mu }_{k,{\zeta }_{k,i}}\oplus {\mathbf{n}}_k\sim P({\mathbf{X}}_k\mid {\mathbf{x}}_{k-1,{\zeta }_{k,i}}),$$

$${\tilde{W}}_{k,i}=\frac{p({\mathbf{Y}}_k\mid {\mathbf{x}}_{k,i})}{p({\mathbf{Y}}_k\mid {\mu }_{k,{\zeta }_{k,i}})}.$$

$$||{\hat{\mbox{\boldmath$\beta$}}}_{k}||_{1}\,=\,\left(\frac{1}{N}\,\sum_{i=1}^{N}V_{k,i}\right)\,\left(\frac{1}{N}\,\sum_{i=1}^{N}\tilde{W}_{k,i}\right)\,.$$

$$E[{\mathbf{X}}_k\mid {\mathbf{Y}}_{1\text{mathchar}58k}]\approx {\mathbf{z}}_k=\sum _{i=1}^N{W}_{k,i}{\mathbf{x}}_{k,i}.$$

$${\hat{\mathbf{X}}}_{k|k}^{\mathrm{MS}}=\overline{{\mathbf{z}}_k},$$

$$\mathrm{MSE}({\hat{\mathbf{X}}}_{k|k}^{\mathrm{MS}},{\mathbf{Y}}_{1\text{mathchar}58k})=||\min \{{\mathbf{z}}_k,{\mathbf{z}}_k^c\}|{|}_1.$$

$$\{{\mathbf{x}}_{k,i}^u,W_{k,i}^u{\}}_{i=1}^{F_k}\stackrel{\text{Unique}}{\leftarrow }\{{\mathbf{x}}_{k,i},W_{k,i}{\}}_{i=1}^N,k=0,\dots ,T.$$

$$P({\mathbf{X}}_s,$$

$$=$$

$$P({\mathbf{X}}_s\mid$$

$$\times \sum _{{\mathbf{X}}_{s+1}}\frac{P({\mathbf{X}}_{s+1}\mid {\mathbf{X}}_s)P({\mathbf{X}}_{s+1}\mid {\mathbf{Y}}_{1\text{mathchar}58T})}{P({\mathbf{X}}_{s+1}\mid {\mathbf{Y}}_{1\text{mathchar}58s})}.$$

$$P({\mathbf{X}}_T\mid {\mathbf{Y}}_{1\text{mathchar}58T})\approx \sum _{i=1}^{F_T}{W}_{T|T,i}{\delta }_{{\mathbf{x}}_{T,i}^u}.$$

Full text

,

1 Introduction

Partially-observed Boolean dynamical systems consist of a Boolean state process, also known as a Boolean network, observed through an arbitrary noisy mapping to a measurement space [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Instances of POBDSs abound in fields such as genomics [12], robotics [13], digital communication systems [14], and more. The optimal recursive minimum mean-square error (MMSE) state estimators for this model are called the Boolean Kalman Filter (BKF) [1] and the Boolean Kalman Smoother (BKS) [15]. These filters have many desirable properties; in particular, it can be shown that the MMSE estimate of the state vector provides both the MMSE and the maximum-a-posteriori (MAP) estimates of each state vector component. Notice that the software tool “BoolFilter” [16] is available under R library for estimation and identification of partially-observed Boolean dynamical systems.

However, for large systems with large number of state variables, the computation of both the BKF and BKS becomes impractical due to large computational and memory requirements. In [3], an approximate sequential Monte-Carlo (SMC) algorithm was proposed to compute the BKF using sequential importance resampling (SIR). By contrast, we develop here SMC algorithms for both the BKF and fixed-interval BKS based on the more efficient auxiliary particle filter (APF) algorithm [17].

The BKF and BKS require for their application that all system parameters be known. In the case where noise intensities, the network topology, or observational parameters are not known or only partially known, an adaptive scheme to simultaneously estimate the state and parameters of the system is required. An exact adaptive filtering framework to accomplish that task was proposed recently in [18], which is based on the BKF and BKS in conjunction with maximum-likelihood estimation of the parameters. In this paper, we develop an accurate and efficient particle filtering implementation of the adaptive filtering framework in [18], which is suitable for large systems.

In the case where the parameter space is discrete (finite), the adaptive filter corresponds to a bank of particle filters in parallel, which is reminiscent of the multiple model adaptive estimation (MMAE) procedure for linear systems [19]. If the parameter space is continuous, then a particle-based version of the Expectation Maximization (EM) algorithm [20] is developed. The computational complexity of EM method arises from three main parts:

The computational complexity of applying smoothing at the E-step. 2. 2.

Memory necessary to store the required matrices and vectors (e.g. the posterior probability vectors) from the E-Step to the M-Step. 3. 3.

The complexity of each iteration in the M-step in which several function evaluations are required.

Our proposed particle-based implementation addresses each of the above issues.

Our application of interest in this paper is to model Boolean gene regulatory networks [12, 22] observed through a single time series of RNA-seq data [23]. Using the POBDS model, we employ the proposed approximate adaptive ML algorithm to estimate the gene expression state simultaneously to the inference of the network topology and noise and expression parameters. Performance is assessed through a series of numerical experiments using the well-known cell-cycle gene regulatory model [24]. The influence of transition noise, expression parameters, and RNA-seq measurement noise (data dispersion) on performance is studied, and the consistency of the adaptive ML filter (i.e., convergence to true parameter values) is empirically established.

The article is organized as follows. In Section 2, the POBDS signal model and the Boolean Kalman Filter and Boolean Kalman Smoother are reviewed, while in Section 3, a detailed description of the APF-based filtering and smoothing algorithms proposed in this paper is provided. In Section 4, the particle-based ML adaptive filter is developed for discrete and continuous parameter spaces. A POBDS model for gene regulatory networks observed though RNA-seq measurements is reviewed in Section 5. Results for the numerical experiments with the cell-cycle network are presented in Section 6. Finally, Section 7 contains concluding remarks.

2 Optimal State Estimators for POBDS

In this section, we review the POBDS model and exact algorithms for computation of its optimal state estimators. For more details see [1, 4, 5, 15]. For a proof of optimality of the BKF, see [18].

We assume that the system is described by a state process $\{{\bf X}_{k};k=0,1,\ldots\}$, where ${\bf X}_{k}\in\{0,1\}^{d}$ is a Boolean vector of size $d$. The sequence of states is observed indirectly through the observation process $\{{\bf Y}_{k};k=1,2,\ldots\}$, where ${\bf Y}_{k}$ is a vector of (typically non-Boolean) measurements. The states are assumed to be updated and observed at each discrete time through the following nonlinear signal model:

[TABLE]

for $k=1,2,\ldots$ Here, ${\bf n}_{k}\in\{0,1\}^{d}$ is Boolean transition noise, “$\oplus$” indicates componentwise modulo-2 addition, ${\bf{f}}_{k}\mathrel{\mathop{\mathchar 58\relax}}\{0,1\}^{d}\rightarrow\{0,1\}^{d}$ is a Boolean function, called the network function, whereas ${\bf{h}}_{k}$ is a general function mapping the current state and observation noise ${\bf v}_{k}$ into the measurement space, for $k=1,2,\ldots$. The noise processes $\{{\bf n}_{k},{\bf v}_{k};k=1,2,\ldots\}$ are assumed to be “white” in the sense that the noises at distinct time points are independent random variables. It is also assumed that the noise processes are independent from each other and from the initial state ${\bf X}_{0}$; their distribution is otherwise arbitrary.

We will assume further that the Boolean process noise ${\bf n}_{k}$ is zero-mode, i.e., ${\bf n}_{k}={\bf 0}$ is the most probable value of the noise vector at each time $k$. This implies that the most likely value of ${\bf X}_{k}$ at each time $k$ is ${\bf{f}}_{k}({\bf X}_{k-1})$ — this could be seen as the counterpart of the zero-mean noise assumption in continuous state-space models. As is the case with nonzero mean noise, nonzero mode noise introduces a systematic error component, which can always be removed by moving it into the function ${\bf{f}}_{k}$. Hence, the state model in (1) is a general model for a first-order Markov Boolean stochastic process. For a specific example, which will be adopted in Section 6, one has $P({\bf n}_{k}(i)=1)=p$, for $i=1,\ldots,d$, and $k=1,2,\ldots$ independently for $i\neq j$. In this case, ${\bf n}_{k}$ is zero-mode if and only if $p\leq 1/2$. The systematic bias introduced in the case $p>1/2$ can be removed by considering the equivalent state process with ${\bf{f}}_{k}^{\prime}={\bf{1}}-{\bf{f}}_{k}$, where ${\bf{1}}$ is the vector with all components equal to 1, and $p^{\prime}=1-p$, i.e., by moving the systematic bias into the model. Therefore, one effectively needs only to consider the case $p\leq 1/2$.

A Boolean estimator ${\bf\hat{X}}_{k|r}$ predicts the Boolean state ${\bf X}_{k}$ based on the sequence of observations ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}r}=({\bf Y}_{1},\ldots,{\bf Y}_{r})$. The estimator ${\bf\hat{X}}_{k|r}$ is called a filter, smoother, or predictor according to whether $k=r$, $k<r$, or $k>r$, respectively. The set of all Boolean estimators for a given $k$ and $r$ shall be denoted by ${\mathcal{}X}_{k|r}$. The (conditional) mean-square error (MSE) of ${\bf\hat{X}}_{k|r}$ given ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}r}$ is:

[TABLE]

We would like to obtain the Boolean MMSE estimator, i.e., a Boolean estimator ${\bf\hat{X}}_{k|r}^{\rm MS}$ such that

[TABLE]

at each value of ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}r}$ (so that it also minimizes the frequentist expected MMSE over all possible realizations of ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}r}$). For a Boolean vector ${\bf{v}}\in\{0,1\}^{d}$, define the binarized vector $\overline{{\bf{v}}}$, such that $\overline{{\bf{v}}}(i)=1$ if ${\bf{v}}(i)>1/2$ and $\overline{{\bf{v}}}(i)=0$ otherwise, for $i=1,\ldots,d$, the complement vector ${\bf{v}}^{c}$, such that ${\bf{v}}^{c}(i)=1-{\bf{v}}(i)$, for $i=1,\ldots,d$, and the $L_{1}$-norm $||{\bf{v}}||_{1}=\sum_{i=1}^{d}|{\bf{v}}(i)|$. It can be proved [18] that the solution to (3) is given by

[TABLE]

with optimal MMSE

[TABLE]

where the minimum is computed componentwise.

The optimal Boolean MMSE estimator will be called a Boolean Kalman Filter (BKF), Boolean Kalman Smoother (BKS), or Boolean Kalman Predictor (BKP), according to whether $k=r$, $k<r$, or $k>r$, respectively. The terminology is justified by the fact that we seek the MMSE estimator for a nonstationary process, as in the case of the classical Kalman Filter, as opposed to, say, the Maximum-A-Posteriori (MAP) estimator, which is more common in discrete settings. Interestingly, we can show that each component of the MMSE estimator ${\bf\hat{X}}_{k|r}^{\rm MS}(i)$ is both the MMSE and the MAP estimator of the corresponding state variable ${\bf X}_{k}(i)$, for $i=1,\ldots,d$. Perhaps surprisingly, the MAP estimator ${\bf\hat{X}}_{k|r}^{\rm MAP}$ does not enjoy in general the property that ${\bf\hat{X}}_{k|r}^{\rm MAP}(i)$ is the MAP estimator of ${\bf X}_{k}(i)$, for $i=1,\ldots,d$. In cases where optimal estimation performance for each component of ${\bf\hat{X}}_{k|r}$ is required (e.g., estimating the state of each gene in a gene regulatory network), then this is an important distinction.

Let $({\bf x}^{1},\ldots,{\bf x}^{2^{d}})$ be an arbitrary enumeration of the possible state vectors, define the state conditional probability distribution vector ${\bf\Pi}_{k|r}$ of length $2^{d}$ via

[TABLE]

and let $A\,=\,\left[{\bf x}^{1}\cdots{\bf x}^{2^{d}}\right]$ be a matrix of size $d\times 2^{d}$. Then it is clear that $E\left[{\bf X}_{k}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}r}\right]=A{\bf\Pi}_{k|r}$, so it follows from (4) and (5) that

[TABLE]

with optimal MSE

[TABLE]

The distribution vector ${\bf\Pi}_{k|r}$ can be computed by a matrix-based procedure similar to the “forward-backward” algorithm [25]. Briefly, let ${M_{k}}$ of size $2^{d}\times 2^{d}$ be the transition matrix of the Markov chain defined by the state model:

[TABLE]

Additionally, given a value of the observation vector ${\bf Y}_{k}$ at time $k$, let $\,{T_{k}}({\bf Y}_{k})$ be a diagonal matrix of size $2^{d}\times 2^{d}$ defined by:

[TABLE]

where $p$ is either a probability density or a mass function, according to the nature of the measurement ${\bf Y}_{k}$.

At this point, we distinguish two cases. The first is a recursive implementation of the Boolean Kalman Filter (BKF), which does not need a backward iteration, and can be iterated forward as new observations arrive, for as long as desired. In this case, we use (7) and (8) with $r=k$ to get the optimal filter estimator and its minimum MSE [1]. The entire procedure is given in Algorithm 1.

The second case is a fixed-interval Boolean Kalman Smoother, where a fixed batch of observations ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T}=({\bf Y}_{1},\ldots,{\bf Y}_{T})$ of length $T$ is available, and it is desired to obtain estimates of the state at all points in the interval $k=1,\ldots,T$. In this case, a backward iteration will be needed (unless $k=T$). Define the probability distribution vector ${\bf\Delta}_{k|s}$ of length $2^{d}$ via

[TABLE]

for $s=0,\ldots,T$, where ${\bf\Delta}_{T|T}$ is defined to be ${\bf 1}_{d\times 1}$, the vector with all components equal to 1. It can be shown that

[TABLE]

where “$\,\bullet\,$” denotes componentwise vector multiplication. We then use (7) and (8) with $r=T$ to get the optimal smoothed estimator and its minimum MSE [15]. The entire procedure is given in Algorithm 2.

3 Particle Filters for State Estimation

When the number of states is large, the exact computation of the BKF and the BKS becomes intractable, due to the large size of the matrices involved, which each contain $2^{2d}$ elements, and approximate methods must be used, such as sequential Monte-Carlo methods, also known as particle filter algorithms which have been successfully applied in various fields [26, 27, 28, 29, 30]. In the next subsections, we describe particle filter implementations of the BKF and BKS.

3.1 Auxiliary Particle Filter Implementation of the BKF (APF-BKF)

The basic algorithm to perform particle filtering is called sequential importance resampling (SIR). Importance sampling is used when direct sampling of the target distribution is difficult. The idea is to approximate the target distribution $p({\bf x})$ using sample points (“particles”) $\{{\bf x}_{i}\}_{i=1}^{N}$ drawn from a proposal distribution $q({\bf x})$, which is easier to sample than the target distribution. The discrepancy created by sampling from $q({\bf x})$ instead of $p({\bf x})$ is compensated by weighting each particle. After a few iterations of the algorithm, a condition in usually reached where only few of the particles have significant weights, whereas most particles have negligible weight. To address this degeneracy problem, SIR performs resampling of the particles, whereby a fresh set of particles is drawn (with replacement) from the approximate current posterior distribution.

The original particle filtering implementation of the BKF in [3] was based on the SIR algorithm; we therefore call it the SIR-BKF algorithm. We present here a more sophisticated implementation based on the Auxiliary Particle Filter (APF) of [17]. The APF algorithm can be seen as a variation of SIR, and is thus also known as auxiliary SIR (ASIR). Basically, APF is a look-ahead method that at time step $k-1$ tries to predict the location of particles with high probability at time $k$, with the purpose of making the subsequent resampling step more efficient. Without the look-ahead, the basic SIR algorithm blindly propagates all particles, even those in low probability regions. As put in [31], “it is natural to ask whether it is possible to employ knowledge about the next observation before resampling to ensure that particles which are likely to be compatible with that observation have a good chance of surviving.”

The APF algorithm augments the state vector to $({\bf X}_{k},\zeta_{k})$, where $\zeta_{k}$ is an auxiliary variable. Particles are drawn from the filtering distribution $P({\bf X}_{k},\zeta_{k}\mid{\bf Y}_{k})$ (to be specified below), and the auxiliary variable is simply dropped to obtain particles from $P({\bf X}_{k}\mid{\bf Y}_{k})$. Given particles $\{{\bf x}_{k-1,i}\}_{i=1}^{N}$ at time $k-1$, with associated weights $\{W_{k-1,i}\}_{i=1}^{N}$, the APF algorithm defines

[TABLE]

for $\zeta_{k}=1,\ldots,N$. The auxiliary variable functions thus as an index for the particles at the previous time point. As will be seen below, sampling from (13) will have the effect of “selecting” the particles that are compatible with the observation at time $k$.

One can sample from (13) by using SIR on the following approximation:

[TABLE]

for $\zeta_{k}=1,\ldots,N$, where $\mu_{k,i}$ is a characteristic of ${\bf X}_{k}$ given ${\bf x}_{k-1,i}$, which can be the mean, the mode or even a sample from $P({\bf X}_{k}\mid{\bf x}_{k,i})$ [17]. In our implementation, we use the mode:

[TABLE]

for $i=1,\ldots,N$, where we used (1). The approximation on the right is accurate as long as the “noise intensity” is low, i.e., the probability of nonzero ${\bf n}_{k}$ is small.

Sampling from (13) is done in two steps. In the first step, $\{\mu_{k,i}\}_{i=1}^{N}$ is obtained from the particles $\{{\bf x}_{k-1,i}\}_{i=1}^{N}$ using (15) and the first-stage weights $\{V_{k,i}\}_{i=1}^{N}$ are computed as:

[TABLE]

for $i=1,\ldots,N$. In the second step, the auxiliary variables $\{\zeta_{k,i}\}_{i=1}^{N}$ (i.e., the indices of the selected particles) are obtained as a sample from the discrete distribution defined by $\{V_{k,i}\}_{i=1}^{N}$ (after proper normalization). For example, if $N=4$ and $V_{k,1}=V_{k,2}$, $V_{k,3}=V_{k,4}$, and $V_{k,1}=2V_{k,3}$, then the indices $\zeta_{k,0},\ldots,\zeta_{k,4}$ will be independent and each will be twice as likely to be 1 or 2 than 3 or 4. We denote this by $\{\zeta_{k,i}\}_{i=1}^{N}\sim{\rm Cat}(\{V_{k,i}\}_{i=1}^{N})$, where “Cat” stands for the categorical (discrete) distribution.

Finally, the new particles $\{{\bf x}_{k,i}\}_{i=1}^{N}$ and associated second-stage weights $\{\tilde{W}_{k,i}\}_{i=1}^{N}$ can be obtained as follows:

[TABLE]

[TABLE]

It can be shown that the unbiased estimator of the unnormalized posterior probability at each time step can be obtained by [32]

[TABLE]

This quantity will be needed in Section 4 when the particle filter for maximum-likelihood adaptive estimation is discussed.

Given the normalized second-stage weights $W_{k,i}=\tilde{W}_{k,i}/\sum_{j=1}^{N}\tilde{W}_{k,j},\,i=1,\ldots,N$, one can write

[TABLE]

From (4) and (5), it follows that the MMSE state estimate and conditional MSE at time step $k$ are approximated as:

[TABLE]

with optimal MMSE

[TABLE]

The entire procedure of APF-BKF is summarized in Algorithm 3.

3.2 Auxiliary Particle Filter Implementation of the BKS (APF-BKS)

There are a few different approximate Monte-Carlo smoothing methods in the literature of nonlinear and non-Gaussian systems [33, 29, 34]. It should be noted that some of these particle smoother methods suffer from degeneracy problems or can only be applied in a few special conditions (such as MC with good forgetting properties). We follow an approach similar to the well-known fixed-interval smoother of [35] to approximate the Boolean Kalman Smoother.

As described in Section 2, a fixed-interval smoother is a forward-backward method, such that the filtering distributions ${\bf\Pi}_{k|k}$ for $k=0,1,\ldots,T$ are computed in the forward step, and the smoothed distributions ${\bf\Pi}_{k|T}$ are found in a backward step. The forward process is obtained here by running the APF-BKF algorithm of Section 3.1, while the backward process is performed by correcting the filtering weights in the backward iteration. We explain next how the backward step is applied efficiently.

First, assume $\{{\bf x}_{k,i},W_{k,i}\}$, $k=0,\ldots,T$ are the forward particles and weights obtained by the APF-BKF algorithm for the sequence of measurements ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T}$. Due to the finite number of states in the POBDS, one can compute unique particles and their associated weights at different time steps as:

[TABLE]

where ${F_{k}}$ is the number of unique forward particles, and ${\bf x}^{u}_{k,i}$ is the $i$th unique particle with aggregated weight $W^{u}_{k,i}$, both at time step $k$.

The backward process is based on the following equation:

[TABLE]

where $s<T$ and $P({\bf X}_{s+1}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T})$ is the smoothed distribution at time step $s+1$. The summation over ${\bf X}_{s+1}$ in both sides of equation (24) results in

[TABLE]

As we mentioned before, the filter and smoother estimate at final time $T$ are the same. Therefore, the smoothed weights $W_{T|T,i}$ are defined in the same way as the forward unique weights $W_{T,i}^{u}$, so that

[TABLE]

Now, using equation (25), the the smoothed weights at time $s<T$ can be obtained as:

[TABLE]

The smoothed weights are obtained by solving equation (27) in a backward fashion using the terminal condition $W_{T|T,j}=W_{T,j}^{u}$, $j=1,\ldots,F_{T}$. The computational complexity of equation (27) is of order $O(F_{s}\times F_{s+1})$ which can be much smaller than $O(N\times N)$ in practice.

Using the smoothed weights, we can write

[TABLE]

From (4) and (5), it follows that the MMSE state estimate and conditional MSE at time step $k$ are approximated as:

[TABLE]

with optimal MMSE

[TABLE]

The entire procedure is summarized in Algorithm 4.

This particle smoother is an efficient method for state estimation, as will be shown in Section 6, but it is not appropriate for parameter estimation, as we will argue in the next section. A different particle smoother will be used in the next section to perform continuous parameter estimation.

4 Particle Filters For Maximum-Likelihood Adaptive Estimation

Suppose that the nonlinear signal model in (1) is incompletely specified. For example, the deterministic functions ${\bf{f}}_{k}$ and ${\bf{h}}_{k}$ may be only partially known, or the statistics of the noise processes ${\bf n}_{k}$ and ${\bf v}_{k}$ may need to be estimated. By assuming that the missing information can be coded into a finite-dimensional vector parameter $\theta\in\Theta$, where $\Theta$ is the parameter space, we propose next particle filtering approaches for simultaneous state and parameter estimation for POBDS. For simplicity and conciseness, we consider two cases: a Boolean Kalman Filter algorithm with finite (i.e., discrete) $\Theta$ and a Boolean Kalman Smoother algorithm with $\Theta\subseteq R^{m}$, but the algorithms can be modified and even combined to perform other estimation tasks. Exact algorithms for such filters can be found in [18].

4.1 APF Implementation of the Discrete-Parameter ML Adaptive BKF (APF-DPMLA-BKF)

In this case, $\Theta=\{\theta_{1},\theta_{2},\ldots,\theta_{M}\}$. Given the observations ${\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}k}=\{{\bf Y}_{1},\ldots,{\bf Y}_{k}\}$ up to time $k$, the log-likelihood function can be written as

[TABLE]

for $i=1,\ldots,M$, where \|{\mbox{\boldmath\beta}}_{k}^{\theta_{i}}\|_{1}=p_{\theta_{i}}({\bf Y}_{k}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}k-1}) can be approximated by running the APF-BKF algorithms discussed in the previous section tuned to parameter vector $\theta_{i}$.

The approximate log-likelihood is updated via

[TABLE]

with $\hat{L}_{0}(\theta_{i})=0$, for $i=1,\ldots,M$, and the ML estimator for both parameter and state at time $k$ can be directly obtained by running $M$ particle filters in parallel, each tuned to a candidate parameter $\theta_{i}$, for $i=1,\ldots,M$:

[TABLE]

[TABLE]

for $k=1,2,\ldots$

The computation in (32)–(34) is parallelized, on-line, and entirely recursive: as a new observation at time $k+1$ arrives, the ML estimator can be updated easily without restarting the computation from the beginning. The procedure is summarized in Figure 1 and Algorithm 5.

4.2 APF Implementation of the Continuous-Parameter ML Adaptive BKS

(APF-CPMLA-BKS)

Here, $\Theta\subseteq R^{m}$, and the approach developed in the last subsection for discrete parameter spaces is not directly applicable. There are two options: 1) discretize the parameters using a suitable quantization grid and apply an approach similar to the one in the last section; 2) attempt to obtain a good approximation of the MLE in the continuous parameter space directly. In this section, we describe how to implement the second option using the expectation-maximization (EM) algorithm for a particle filter implementation of a fixed-interval Boolean Kalman Smoother.

In our case, maximum likelihood estimation attempts to find the value of $\theta$ that maximizes the “incomplete” log-likelihood function $L_{k}(\theta)=\log\,p_{\theta}({\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T})$. The EM algorithm considers instead the “complete” log-likelihood function $\log p_{\theta}({\bf X}_{0\mathrel{\mathop{\mathchar 58\relax}}T},{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T})$, which includes the unknown state sequence, the assumption being that maximising the complete log-likelihood is easier than maximising the incomplete one (the reader is referred to [18] for more information on the EM algorithm for POBDS).

The EM algorithm obtains a sequence of parameter estimates $\{\theta^{(n)};n=0,1,\ldots\}$. Given the current estimate $\theta^{(n)}$, the algorithm obtain the next estimate $\theta^{(n+1)}$ in the sequence by computing (E-step) the function (see [18]):

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and then maximizing (M-step) this function:

[TABLE]

In [18] this computation is carried out exactly. For large systems, this is impractical, for the following reasons:

The E-Step (computing the $Q$ function) requires performing a Boolean Kalman Smoother, which is too expensive computationally. 2. 2.

The transition matrix and filtered and smoothed posterior probability vectors at all time steps must be stored, demanding large amounts of memory. 3. 3.

In certain cases, such as when ${\bf{f}}_{k}$ and ${\bf{h}}_{k}$ are linear in the parameter vector $\theta$, it is possible to maximize ${Q}(\theta,\theta^{(n)})$ using closed-form expressions (e.g. [34]). However, in general, one needs to resort to gradient-based optimization methods in the M-step. This requires evaluating ${Q}(\theta,\theta^{(n)})$ and computing its derivatives, which is analytically intractable.

To address all these issues, we develop our EM algorithm based on the Forward Filter Backward Simulation [27] method. This method tries to capture the most probable state trajectories and use them to find the smoothing particles. The method contains two main steps: 1) Forward Step: the APF-BKF algorithm is employed to obtain the particles and their weights from time 0 to $T$ ($\{{\bf x}_{1\mathrel{\mathop{\mathchar 58\relax}}T,i},W_{1\mathrel{\mathop{\mathchar 58\relax}}T,i}\}_{i=1}^{N}$). 2) Backward Step: the backward simulation procedure, which is explained in detail in the sequel, computes $N$ trajectories $\{\tilde{\bf{x}}_{0\mathrel{\mathop{\mathchar 58\relax}}T,i}\}_{i=1}^{N}\,\sim\,P({\bf X}_{0\mathrel{\mathop{\mathchar 58\relax}}T}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}T})$, where

[TABLE]

Based on equation (40), smoothed particles can be obtained using the FFBSi method, by means of the following backward procedure:

[TABLE]

for $i=1,\ldots,N$ and $s=T-1,\ldots,0$, where $\{\tilde{\bf{x}}_{s,i}\}_{i=1}^{N}$ are the smoothed particles at time step $s$.

The backward process starts by resampling $N$ particles $\{\tilde{\bf{x}}_{T,i}\}_{i=1}^{N}$ from the unique forward particles $\{{\bf x}^{u}_{T,i}\}_{i=1}^{F_{T}}$ at time step $T$ using the forward weights $\{W_{T,i}^{u}\}_{i=1}^{F_{T}}$ as:

[TABLE]

Now, to obtain $N$ smoothed particles at time step $s<T$, let $\{\tilde{\bf{x}}_{s+1,i}\}_{i=1}^{N}$ be the smoothed particles at time $s+1$, and let

[TABLE]

where $S_{s+1}$ specifies the number of unique smoothed particles at time $T$, and $\xi_{s+1}^{j}$ contains the indexes of the $j$-th unique smoothed particles before shrinkage and reordering. Notice that $N=\sum_{j=1}^{S_{s+1}}|\xi_{T}^{j}|$. For the $j$-th unique smoothed particle at time step $s+1$, one can use the fact that $P({\bf X}_{s}\mid{\bf X}_{s+1},{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}s})\propto P({\bf X}_{s+1}\mid{\bf X}_{s})\,P({\bf X}_{s}\mid{\bf Y}_{1\mathrel{\mathop{\mathchar 58\relax}}s})$ to compute the following weights

[TABLE]

for $i=1,\ldots,F_{s}$, and draw $|\xi_{s+1}^{j}|$ particles as:

[TABLE]

Repeating the above process for $j=1,\ldots,S_{T}$ and $s=T-1,\ldots,0$ results in $N$ trajectories from the joint smoothed distribution $\{\tilde{\bf{x}}_{0\mathrel{\mathop{\mathchar 58\relax}}T,i}\}_{i=1}^{N}$. Notice that the computational complexity of this method is of order $O(S_{s}\times F_{s})$ at time step $s$, which can be much smaller than the computational complexity of the optimal smoother which is $O(N\times N)$.

Given the $N$ trajectories $\{\tilde{\bf{x}}_{0\mathrel{\mathop{\mathchar 58\relax}}T,i}\}_{i=1}^{N}$ obtained by running the Forward Filter Backward Simulation tuned to parameter $\theta^{(n)}$, equations (36)-(38) can be approximated as:

[TABLE]

Thus, one can approximate the $Q$ function in equation (35) as:

[TABLE]

In [36], similar equations are derived for the Hammerstein-Wiener model structure. The computational complexity of evaluating the $\hat{Q}$ function is only of order $O(Nk)$, which results in large savings in the computation of the gradient of $\hat{Q}$ in the M-Step of the EM algorithm. In Section 5, expressions for the gradients are given in the special case where the observations consist of RNA sequencing data. Finally, as regards to memory, the only values that must be stored from the E-Step to be used in the M-Step are the $N$ smoothed trajectories (storing filter weights or particles is not necessary). In Section 6, the effect of the value of $N$ on performance will be discussed.

The steps of the EM adaptive filter are as follows. Initially, $N$ smoothed trajectories are obtained using the developed FFBSi method tuned to a initial parameter guess $\theta^{(0)}$ to compute $\hat{Q}(\theta,\theta^{(0)})$ (E-Step). The aforementioned gradient-descent procedure is applied to find the best parameter $\theta^{(1)}$ that maximizes $\hat{Q}(\theta,\theta^{(0)})$ with $\theta^{(0)}$ fixed (M-Step). The obtained parameter vector is set as the parameter for the particle smoother for the next run, and the process continues until there is no significant change in parameter estimates between two consecutive steps, yielding the final parameter estimate $\theta^{\rm ML}$. Then the smoothed state estimates can be obtained by performing an APF-BKS tuned to parameter $\theta^{\rm ML}$. The procedure is summarized in Figure 2 and Algorithm 6.

5 Gene Regulatory Network and

RNA-Seq Measurement Models

The algorithms developed in the previous section apply to the general partially-observed Boolean dynamical systems signal model in (1). In this section, we describe a specific instance of that model, which allows the application of the methodology to Boolean gene regulatory networks observed through next-generation sequencing measurements. The gene regulatory network model corresponds to the state model in the general POBDS model, whereas the RNA-seq measurement model corresponds to the observation model.

5.1 Gene Regulatory Network Model

This model is motivated by gene pathway diagrams commonly encountered in biomedical research. The network function in (1) is assumed to be time-invariant and expressed as ${\bf{f}}=(f_{1},\ldots,f_{d})$, where each component $f_{i}\mathrel{\mathop{\mathchar 58\relax}}\{0,1\}^{d}\rightarrow\{0,1\}$ is a Boolean function given by

[TABLE]

where $a_{ij}$ and $b_{i}$ are system parameters. The former can take three values: $a_{ij}=+1$ if there is positive regulation (activation) from gene $j$ to gene $i$; $a_{ij}=-1$ if there is negative regulation (inhibition) from gene $j$ to gene $i$; and $a_{ij}=0$ if gene $j$ is not an input to gene $i$. The latter specifies regulation biases and can take two values: $b_{i}=+1/2$ or $b_{i}=-1/2$. The network function is depicted in Figure 3, where the threshold units are step functions that output 1 if the input is nonnegative, and 0, otherwise.

The process noise ${\bf n}_{k}$ in (1) is assumed to have independent components with $P({\bf n}_{k}(i)=1)=p$, for $i=1,\ldots,d$, and $k=1,2,\ldots$ The noise parameter $0\leq p\leq 0.5$ gives the amount of “perturbation” to the Boolean state process; the closer it is to $p=0.5$, the more chaotic the system will be, while a value of $p$ close to zero means that the state trajectories are nearly deterministic, being governed tightly by the network function.

5.2 RNA-Seq Measurement Model

Next-generation sequencing (NGS) technologies are able to sequence millions of short DNA fragments in parallel; the length and number of reads vary with the specific technology [37]. The application of NGS technology to transcriptional profiling is called RNA-seq, which records how frequently each transcript is represented in a sequence sample [38]. RNA-seq is a probe-free approach that can capture any relevant transcript present in a sample, without the need of prior knowledge about the target sequence.

Let ${\bf Y}_{k}=({\bf Y}_{k}(1),\ldots,{\bf Y}_{k}(d))$ be a vector containing the RNA-seq data at time $k$, for $k=1,2,\ldots$ such that ${\bf Y}_{k}(j)$ is the read count corresponding to transcript $j$ in a single-lane platform, for $j=1,\ldots,d$. We assume conditional independence of the transcript counts given the state,

[TABLE]

and adopt the negative binomial model for each count,

[TABLE]

where $\Gamma$ denotes the Gamma function, and $\phi_{j},\lambda_{j}>0$ are the real-valued inverse dispersion parameter and mean read count of transcript $j$, respectively, for $j=1,\ldots,d$. The inverse dispersion parameter models observation noise; the smaller it is, the more variable the measurements are.

Now, recall that, according to the Boolean state model, there are two possible states for the abundance of transcript $j$: high, if ${\bf x}(j)=1$, and low, if ${\bf x}(j)=0$. Accordingly, we model the parameter $\lambda_{j}$ in log-space as:

[TABLE]

where the parameter $s$ is the sequencing depth (which is instrument-dependent), $\mu>0$ is the baseline level of expression in the inactivated transcriptional state, and $\delta_{j}>0$ expresses the effect on the observed RNA-seq read count as gene $j$ goes from the inactivated to the activated state, for $j=1,\ldots,d$.

Based on equations (51)–(53), given the particles $\tilde{\bf{x}}_{k,i}$, one can compute $P({\bf Y}_{k}\mid{\bf X}_{k}=\tilde{\bf{x}}_{k,i})$ as :

[TABLE]

The RNA-seq measurement model parameters are thus the sequencing depth $s$, the baseline expression level $\mu$, the transcript-dependent differential expression levels $\delta_{j}$, for $j=1,\ldots,d$, and the transcript-dependent inverse dispersion parameters $\phi_{j}$, for $j=1,\ldots,d$. These are all continuous parameters.

6 Numerical Experiments

In this section, we carry out detailed numerical experiments to assess the performance of the developed particle-based methods. We base our experiments on the well-known Mammalian Cell-Cycle network [39]. The pathway diagram for this network is presented in Figure 4. The state vector is ${\bf x}=$ (CycD, Rb, p27, E2F, CycE, CycA, Cdc20, Cdh1, UbcH10, CycB). The gene interaction parameters $a_{ij}$ can be read off Figure 4 easily. As an example, Rb is activated by p27, and is inactivated by CycD, CycE, CycA, CycB. These interactions can be expressed in terms of interaction parameters as: $a_{21}=-1$, $a_{22}=0$, $a_{23}=+1$, $a_{24}=0$, $a_{25}=-1$, $a_{26}=-1$, $a_{27}=0$, $a_{28}=0$, $a_{29}=0$ and $a_{2\,10}=-1$.

In all numerical experiments to follow, we assume the same fixed set of “true” values for the system parameters, summarized in Table 2.

6.1 Experiment 1: State Estimation

In this section, the state estimation performance of the APF-BKF and APF-BKS is compared to that of the exact BKF and BKS, respectively. Given that the cell-cycle network comprises 10 genes, the size of the transition and update matrices required by both the BKF and BKS is $2^{10}\times 2^{10}$. As a result, the computational cost of the BKF and BKS is high. Table 3 shows the average rate of correct state estimation over 1000 independent runs for a time series with length 100.

As expected, the performance of both the AFP-BKF and APF-BKS is higher for large number of particles. However, the improvement is significantly larger by moving from 200 to 1000 particles in comparison to moving from 1000 to 5000 particles. One can also see the reduction in performance of all filters and smoothers as process noise or dispersion in measurements increases, which both make the estimation process more challenging. Also as expected, the BKS and APF-BKS outperform the BKF and APF-BKF, respectively, due to the availability of more data for estimation.

The average time to run the various algorithms for time series of length 100 and different number of particles is displayed in Figure 5. Here, $p=0.05$ and $\phi=5$. The average computational time is measured on a PC with an Intel Core i7-4790 CPU@3.60 GHz clock and 16 GB of RAM. The results show the very large computational savings afforded by the APF-BKF and APF-BKS in comparison to the exact algorithms.

6.2 Experiment 2: Incomplete Network Topology

In this experiment, we assume that the interaction between genes Rb and E2F, or equivalently the gene interaction parameter $a_{42}$, is unknown, and all other parameters are known. Since this is a discrete parameter estimation problem, the APF-DPMLA-BKF is run, which in this case consists of three APF-BKFs running in parallel — one for each possible kind of interaction (activation, inhibition, no interaction).

Table 4 displays the average accuracy rate in the estimation of the interaction type between Rb and E2F over 100 different runs. We can observe that performance increases with longer time series and larger number of particles, as expected. The performance is better for larger transition noise. The reason is that large transition noise gets the system out of its attractors more often and, as a result, helps the estimation process.

The evolution of estimated state and discrete-parameter for a single sample run of the experiment, for transition noise $p=0.05$, $\phi=5$, and $N=5000$ is displayed in Figure 6. We can observe that the discrete parameter $a_{42}$ is estimated correctly after less than $20$ time step. In addition, we can see that the state estimator of each gene eventually converges to the true state value.

6.3 Experiment 3: Unknown Noise and Expression Parameters

In the final experiment, the Boolean network topology (gene interaction parameters $a_{ij}$ and biases $b_{i}$) is assumed to be completely known, whereas the transition noise parameter $p$, the baseline expression $\mu$, and differential expression $\delta_{i}$, $i=1,\ldots,10$, are unknown. The inverse dispersion parameters are assumed to be $\phi_{i}=5$, $i=1,\ldots,10$.

In order to assess continuous-parameter estimation accuracy, we define the relative distance between estimated and true parameter values as

[TABLE]

where $\theta^{*}$ is the true parameter value, and $R(\theta)$ is the range of parameter $\theta$ assumed in the M-step of the APF-CPMLA-BKS algorithm. Here, the range is $R(p)=[0,0.5]$ for the transition noise $p$, $R(\mu)=[0,2]$ for the baseline expression $\mu$, and $R(\delta_{i})=[0.1,10]$ for the differential expression $\delta_{i}$, $i=1,\ldots,10$.

A new version of the “augmented Lagrange method” [40] is used for optimization in the M-Step of the particle-based EM algorithm. The gradient vector at each step is computed based on Table 1. The procedure terminates when the maximum of the absolute values of the changes in the parameter estimates in two consecutive iterations gets smaller than $10^{-4}$.

The average relative distance between estimated and true parameter values over 100 independent runs for different inverse dispersion parameters and time series lengths are plotted in Figure 7. As expected, the performance of APF-CPMLA-BKS improves steadily as time goes on. Performance improves by increasing the number of particles; however, the curves get close to each other as the length of the time series increases. All curves show a decreasing trend, which indicates that the parameter estimates become arbitrarily close to the true values for sufficiently long time.

7 Conclusion

In this paper, we introduced approximate particle-based algorithms for state and simultaneous state and parameter estimation for large partially-observed Boolean dynamical systems. For approximate state estimation, filtering and smoothing methods based on auxiliary particle filtering (APF) were developed to approximate the optimal BKF and BKS. These algorithms are called APF-BKF and APF-BKS, and are original contributions of this work.

Moreover, we considered the case where some of the parameters may not be known. In the discrete parameter case, an adaptive filter scheme is developed based on APF-BKF algorithms running in parallel. For continuous parameter problems, a particle-based EM algorithm for POBDS was presented.

The methodology was applied to a model of Boolean gene regulatory networks observed through RNA sequencing data. The numerical experiments with a cell-cycle Boolean network demonstrated the ability of the proposed methodologies to efficiently estimate the state and also the parameters of the large Boolean regulatory network observed through noisy measurements.

Acknowledgment

The authors acknowledge the support of the National Science Foundation, through NSF award CCF-1320884.