On Robust State Estimation of Gene Networks

Abstract

Gene networks in biological systems are not only nonlinear but also stochastic due to noise corruption. How to accurately estimate the internal states of the noisy gene networks is an attractive issue to researchers. However, the internal states of biological systems are mostly inaccessible by direct measurement. This paper intends to develop a robust extended Kalman filter for state and parameter estimation of a class of gene network systems with uncertain process noises. Quantitative analysis of the estimation performance is conducted and some representative examples are provided for demonstration.

Introduction

Biochemical networks^1,2 involving metabolic networks, gene regulatory networks and signal transduction networks in biological systems play important roles in diagnosis of disease such as cancer and auto-immunity. Systems biology^2–4 is quite different from traditional biology; it has been developed recently to understand biological systems from system level. Researchers have been devoted to design and construct of biological models by engineering methods and molecular biology techniques. The microarray technology using high-throughput method to measure a large number of gene expression states has also been attempted in the recent decade. According to the measured data it is possible to reconstruct the structure of biochemical networks, quantitatively analyze and systematically design and simulate the system behavior in silico.

In the literature,^1,5,6 some mathematical models suitable for describing behaviors of the biological systems have been proposed. That kind of models can be classified into two major categories: logical model in the discrete-time domain and differential equation set in the continuous-time domain. Different from the deterministic case, biochemical networks of the real biological systems are generally non-ideal and are invariably noisy. For modeling accuracy, the influence resulting from noise contamination should not be ignored. In general, molecular noises involve intrinsic noises resulting from molecular birth and death and extrinsic noises due to environmental perturbations. ⁷ The stochastic model was developed to characterize the biological systems with intrinsic and extrinsic molecular noises. In, ⁸ the authors have presented a method for measuring performance robustness and presented two mechanisms to cope with the noise and process uncertainties.

With regard to control of biological systems, nonlinear feedback control strategies were applied to regulate the steady state of biological systems in.^9–11 Other issues received increasing attention for noisy biological systems are development of control strategies while ensuring robust stability and filtering ability. In,^12,13 a robust filtering circuit design has been developed by regulating parameters for gene networks to reduce intrinsic and extrinsic molecular noises. An adaptive control design method was also proposed in. ¹⁴

Before performing feedback control, the system states should all be available. However, the internal states of most biological systems can only be observed partially. Under the situation, a state estimator is appropriate to reconstruct the full state in the noisy environment. The Kalman filter (KF) ¹⁵ is the one for the purpose in the filed of engineering which has been well applied from system and control to signal processing for decades. However, applications of the extended KF (EKF) in state estimation of biochemical networks are rarely found. ¹⁶ Until recently, the EKF has been attempted to estimate parameters of the gene regulatory networks. ¹⁷ Moreover, a state observer was actually established using the EKF, based on the fluorescence probe model, a dynamic state model of the plant cell bioreactor, and online green fluorescent protein fluorescence measurement. ¹⁸

While there were a few papers dealing with the issue of state estimation for biological networks, most of the approaches were based on the traditional Kalman filtering theory which assumed that the noise covariances including the process noise and measurement noise have been precisely known as a priori.^16–18 The KF is shown to be the optimal state estimator against noise with Gaussian distributions. However, in the biological systems, the noise distribution may not be Gaussian, its autocorrelations may not be known exactly, or even uneasily to be precisely modeled. ¹⁹ In this paper, an EKF for robustly filtering the states and parameters of the noisy gene networks is introduced. Quantitative error analysis for that kind of systems is presented in details. On the basis of the results obtained one would be able to identify effect of the filtering gain and the sizes of noise uncertainties to estimation performance. Two numerical examples have been conducted to verify the proposed design.

To clarify the notation, throughout this paper, let the vector norm of the real vector x ∊ℝⁿ, denoted by ∥x∥, be defined as $\Vert x\Vert =\sqrt{E\hspace{0.17em}\left(x^{T}x\right)}$ with E(·) denoting the operation of expectation.

Lemma 1: The induced matrix norm ∥A∥ corresponding to the vector norm of $\Vert x\Vert =\sqrt{E\hspace{0.17em}\left(x^{T}x\right)}$ with x ∊ ℝⁿ is given by

$$\Vert A\Vert ={\sigma }_1$$

where σ₁ is the maximum singular value of A, ie, $\Vert A\Vert =\sqrt{{\lambda }_{\max }\hspace{0.17em}\left(A^{T}A\right)}$ .

Proof: By definition, the induced norm of the operator A is defined by

$$\Vert A\Vert =\underset{\Vert x\Vert =1}{\sup }\hspace{0.17em}\Vert Ax\Vert$$

Let the singular value decomposition of A is given by A = UΣV^T where Σ = diag(σ₁, σ₂, σ_n) with σ₁ ≥ σ₂ ≥ … ≥ 0, the unitary matrices U = [u₁ … u_n] and V = [ν₁ … ν_n]. Then

$$\begin{array}{c}\Vert Ax\Vert =\sqrt{E\left(x^T\hspace{0.17em}{A}^T\hspace{0.17em}Ax\right)}\\ =\sqrt{E\left(x^T\sum _{i=1}^n{\sigma }_i^2{\left(u_i{v}_i^T\right)}^T{u}_i{v}_i^{T}x\right)}\\ =\sqrt{E\left(x^T\sum _{i=1}^n{\sigma }_i^2{v}_i^T{v}_{i}x\right)}\\ \le \sqrt{E\left({\sigma }_1^2{x}^T\sum _{i=1}^n{v}_i^T{v}_{i}x\right)}\\ =\sqrt{{\sigma }_1^{2}E\left(x^{T}x\right)}\\ ={\sigma }_1\Vert x\Vert \end{array}$$

That is

$$\underset{\Vert x\Vert =1}{\sup }\hspace{0.17em}\Vert Ax\Vert \le {\sigma }_1$$ (1)

Alternately, if one chooses x′ = e₁ = [1 0 … 0]^T which obviously satisfies ∥x′∥ = 1, then

$$\underset{\Vert x\Vert =1,\forall x\in {ℝ}^n}{\sup }\hspace{0.17em}\Vert Ax\Vert \ge {\sigma }_1\Vert Ax\prime \Vert =\sqrt{E\left({\sigma }_1^{2}x{\prime }^T\hspace{0.17em}x\prime \right)}={\sigma }_1\sqrt{E\left(x{\prime }^T\hspace{0.17em}x\prime \right)}={\sigma }_1$$ (2)

From (1) and (2), it can be concluded that ∥A∥ = σ₁.

Problem Formulation

In real biological systems, biochemical networks such as signal transduction networks, gene regulatory networks, and metabolic networks are invariably noisy. As in engineering area of research, the system dynamic behavior can be mathematically described by stochastic models, which could be further used as a basis for the purpose of analysis and control.

System Modeling

S-system is a type of power-law formalism and is based on a particular type of ordinary differential equations in which the component processes are characterized by power-law functions:^1,5,6 $${\dot{x}}_i={\alpha }_i\prod _{j=1}^{\tilde{n}+\tilde{m}}{x}_j^{g_{ij}}-{\beta }_i\prod _{j=1}^{\tilde{n}+\tilde{m}}{x}_j^{h_{ij}},\hspace{1em}i=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,\hspace{0.17em}\tilde{n}$$ (3) where x _j (j = 1, …,ñ) are dependent variables such as intermediate metabolites and products, x _j $\left(j=\tilde{n}+1,\dots ,\tilde{n}+\tilde{m}\right)$ are independent variables such as substrates and enzymes, α_i ≥ 0 and β_i ≥ 0 are rate constants which denote, respectively, production and degradation effects for dependent and independent variables; g_ij and h_ij are kinetic orders. Gene j will active gene i when the values of kinetic orders are positive and gene j will inhibitive gene i when the values of kinetic orders are negative. Zero values of kinetic orders represent gene j won't affect gene i.

The system model can be expressed in the following generalized nonlinear biochemical dynamics with the stoichiometric equation described by

$$\dot{x}\left(t\right)=SV\left(x\left(t\right),p\right)$$ (4)

where x(t) ∊ℝ_n is a state vector which denotes the concentration of metabolite, mRNA or protein, p = p(α_i, β_i,g_ij, h_ij)∊ ℝ^mis the parameter vector which include rate constants and kinetic orders, S denotes the stoichiometric matrix, and V(·) is a nonlinear function of the reaction rate. It can be expressed in a more general form as

$$\dot{x}\left(t\right)=f\left(x\left(t\right),p\right)$$ (5)

where f(·) ∊ ℝⁿ is a generalized nonlinear function vector.

For biochemical reactions, suppose that there are M intrinsic noise sources and an extrinsic noise, the nonlinear stochastic dynamical system can be described by ¹³

$$\dot{x}\left(t\right)=f\left(x\left(t\right),p\right)+\sum _{i=1}^M{g}_i\left(x\left(t\right)\right)n_i\left(t\right)+w\left(t\right)$$ (6)

where x(t) ∊ ℝⁿ is the state which may denote a vector of protein concentrations of n genes, f(·) and g(·) are nonlinear function vectors, n_i(t) is an intrinsic noise source which is the white noise with zero mean and standard deviation σ_i, w(t) is an extrinsic noise source with zero mean and standard deviation σ_w Equation (6) implies that the system is suffered from intrinsic noise corruption due to M kinetic parameter fluctuations. In gene regulatory networks, the external and internal noises occur independently and they appear randomly. ¹⁹ More precisely, it has been proposed that the pattern of protein concentration growth is stochastic, exhibiting short bursts of variable numbers of proteins at varying time intervals.

Remark: The system with intrinsic noises in (6) can also be rewritten as the following Ito stochastic equation:

$$dx\left(t\right)=f\left(x\left(t\right),p\right)dt+\sum _{i=1}^M{g}_i\left(x\left(t\right),p\right)d{N}_i\left(t\right)+dW\left(t\right)$$ (7)

where N_i and W are standard Wiener processes or Brownian motions with dN_i(t) = n_i(t)dt and dW(t) = w(t)dt with the property $E\left[|N_i(t)-N_i(\tau ){|}^2\right]={\sigma }_i^2|t-\tau |$ and $E\left[|W(t)-W(\tau ){|}^2\right]={\sigma }_w^2|t-\tau |$ . The formulation is widely applicable to the general nonlinear gene network with n genes.

After the stochastic differential system in (6) or (7) is modeled to mimic the realistic behaviors of the object, some design specifications can be proposed for robust design of the system so that the desired behaviors or products can be achieved in spite of intrinsic parameter fluctuations and environmental disturbances.

Estimator Design

In practice, the internal states of biological systems may not be directly accessible. Biochemical process for gene regulatory networks is DNA to mRNA by transcription, mRNA to protein by translation and the generated protein regulates the gene. Not all protein concentrations are directly measurable, ¹⁸ however, one could access the status of individual protein by utilizing an appropriate state estimator. It is proposed here an EKF to estimate the internal states as well as the parameters of concern.

Fundamental Kalman Filter

The KF uses measurements for a dynamic system observed over time that contain noise, and produce estimated values that tend to be closer to the true values of the measurements.

A stochastic linear time-invariant system with measurement can be described by

$$\dot{x}\left(t\right)=Ax\left(t\right)+w_x\left(t\right)$$ (8)

and

$$y\left(t\right)=Cx\left(t\right)+v\left(t\right)$$ (9)

where x(t) is the state vector, y(t) is the output, A and C are, respectively, system and output matrices, and w_x(t) and ν(t) are, respectively, process and measurement noises which are zero mean Gaussian noises with covariance matrices Q and R, respectively.

With regard to the system (8)−(9), which x(t) is to be estimated, the corresponding KF is given by ¹⁵

$$\dot{\widehat{x}}\left(t\right)=A\widehat{x}\left(t\right)+K\left[y\left(t\right)-C\widehat{x}\left(t\right)\right]$$ (10)

where $\widehat{x}\left(t\right)$ is the estimated state and K is the optimal Kalman gain given by

$$K=P{C}^T{R}^{-1}$$ (11)

where P = P^T > 0 satisfies the following Riccati matrix equation:

$$\begin{array}{c}\dot{P}=AP+P{A}^T+Q-P{C}^T{R}^{-1}CP,\\ P_0=E\left({\tilde{x}}_0{\tilde{x}}_0^T\right)\end{array}$$ (12)

with the estimation error state $\tilde{x}\left(t\right)=x\left(t\right)-\widehat{x}\left(t\right)$ .

Extended Kalman Filter

For estimation of the stochastic nonlinear systems, the commonly applied filtering mechanism is EKF which evaluates the partial derivatives at the estimated state value and uses nonlinear functions on the estimate itself.

Suppose that the nonlinear system with measurement is given by

$$\dot{x}\left(t\right)=f\left(x\left(t\right)\right)+w_x\left(t\right)$$ (13)

and

$$y\left(t\right)=h\left(x\left(t\right)\right)+v\left(t\right)$$ (14)

where f(·)and h(·) are nonlinear function vectors with appropriate dimensions.

The EKF is given by ¹⁵

$$\dot{\widehat{x}}\left(t\right)=f\left(\widehat{x}\left(t\right)\right)+K\left(t\right)\left[y\left(t\right)-h\left(\widehat{x}\left(t\right)\right)\right]$$ (15)

where the Kalman gain is

$$K\left(t\right)=P{H}^T\hspace{0.17em}{R}^{-1}$$ (16)

with P = P^T > 0 satisfying

$$\begin{array}{c}\dot{P}=FP+P{F}^T+Q-P{H}^T{R}^{-1}HP,\\ P_0=E\left({\tilde{x}}_0{\tilde{x}}_0^T\right)\end{array}$$ (17)

with

$${F=\frac{\partial f\left(x\left(t\right)\right)}{\partial x\left(t\right)}|}_{x\left(t\right)=\widehat{x}\left(t\right)},{H=\frac{\partial h\left(x\left(t\right)\right)}{\partial x\left(t\right)}|}_{x\left(t\right)=\widehat{x}\left(t\right)}$$

Estimating states and parameters by EKF

In more general, we consider a generalized representation of the stochastic nonlinear dynamical model described by

$$\begin{array}{c}\dot{x}\left(t\right)=f\left(x\left(t\right),p\left(t\right)\right)+w_x\left(t\right),\hspace{1em}{x}_0=E\left[x\left(0\right)\right]\\ \dot{p}\left(t\right)=w_p\left(t\right),\hspace{1em}{p}_0=E\left[p\left(0\right)\right]\\ y\left(t\right)=h\left(x\left(t\right),p\left(t\right)\right)+v\left(t\right)\end{array}$$ (18)

where P_x0 = E[(x(0) − x₀)(x(0) − x₀)T], P_p0 = E[(p(0) − p₀) (p(0) − p₀)T], p(t) ℝ^m denotes the aggregated parameter vector, y(t) ∊ ℝ^r is the measurement output, f(·) and h(·) are nonlinear function vectors. Suppose that the uncorrelated extrinsic noise w_x(t), parameter noise w_p(t) and measurement noise ν(t) are white noises, uncorrelated and satisfy the following properties:

$$\begin{array}{l}E\left[w_x\left(t\right)\right]=E\left[w_p\left(t\right)\right]=E\left[v\left(t\right)\right]=0,\hfill \\ E\left[w_x\left(t\right)w_x^T\left(\tau \right)\right]=Q\delta \left(t-\tau \right),\hfill \\ E\left[w_p\left(t\right)w_p^T\left(\tau \right)\right]=Q_p\delta \left(t-\tau \right),\hfill \\ E\left[v\left(t\right)v^T\left(\tau \right)\right]=R\delta \left(t-\tau \right)\hfill \end{array}$$ (19)

where the noise uncertainties satisfy

$$\Vert Q-Q_0\Vert \le {\epsilon }_1,\hspace{1em}\Vert R-R_0\Vert \le {\epsilon }_2,\hspace{1em}\Vert Q_p-Q_{p0}\Vert \le {\epsilon }_3$$ (20)

where $Q=Q^T>0,\hspace{0.17em}{Q}_P=Q_P^T>0$ and R = R^T > 0 with Q₀, Q_p0 and R₀ being their corresponding nominal parts, and ε₁, ε₃ and ε₂ are finite constants characterizing the respective upper bound of the noise covariance.

For compactness, the matrix format of (18) can next be written in the state-space representation as follows

$$\left[\begin{array}{c}\dot{x}\left(t\right)\\ \dot{p}\left(t\right)\end{array}\right]=\left[\begin{array}{c}f\left(x\left(t\right),p\left(t\right)\right)\\ 0_{m\times 1}\end{array}\right]+\eta \left(t\right),\hspace{1em}y\left(t\right)=h\left(x\left(t\right),p\left(t\right)\right)+v\left(t\right)$$ (21)

where $\eta \left(t\right)={\left[\begin{array}{cc}{w}_x^T\left(t\right)& w_p^T\left(t\right)\end{array}\right]}^T$ . The idea of the EKF is that it operates by propagating the mean and error covariance of the state through time. The EKF for (18) in the matrix form can be represented as

$$\left[\begin{array}{c}\dot{\widehat{x}}\left(t\right)\\ \dot{\widehat{p}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}f\left(\widehat{x}\left(t\right),\widehat{p}\left(t\right)\right)\\ 0_{m\times 1}\end{array}\right]+K\left[y\left(t\right)-\widehat{y}\left(t\right)\right],\hspace{1em}\widehat{y}\left(t\right)=h\left(\widehat{x}\left(t\right),\widehat{p}\left(t\right)\right)$$ (22)

where $\widehat{x}\left(t\right)$ is the estimated state, $\widehat{p}$ is the estimated parameter, ŷ(t) is the estimator output, and K is the estimator gain.

Equations (21) and (22) can be further written as

$$\begin{array}{c}\dot{z}\left(t\right)=g\left(z\left(t\right)\right)+\eta \left(t\right),\hspace{1em}{z}_0={\left[x^T\left(0\right)\hspace{1em}{p}^T\left(0\right)\right]}^T\\ y\left(t\right)=h\left(z\left(t\right)\right)+v\left(t\right)\end{array}$$ (23)

and the state estimate equation

$$\begin{array}{c}\dot{\widehat{z}}\left(t\right)=g\left(\widehat{z}\left(t\right)\right)+K\left[y\left(t\right)-\widehat{y}\left(t\right)\right],\\ \widehat{y}\left(t\right)=h\left(\widehat{z}\left(t\right)\right)\end{array}$$ (24)

where

$$z\left(t\right)={\left[x^T\left(t\right)\hspace{1em}{p}^T\left(t\right)\right]}^T,\hspace{1em}\widehat{z}\left(t\right)={\left[{\widehat{x}}^T\left(t\right)\hspace{1em}{\widehat{p}}^T\left(t\right)\right]}^T,\hspace{1em}g\left(z\left(t\right)\right)=\left[\begin{array}{c}f\left(x\left(t\right),p\left(t\right)\right)\\ 0_{m\times 1}\end{array}\right],\hspace{1em}g\left(\widehat{z}\left(t\right)\right)=\left[\begin{array}{c}f\left(\widehat{x}\left(t\right),\widehat{p}\left(t\right)\right)\\ 0_{m\times 1}\end{array}\right]$$ (25)

The maximal covariances of the extrinsic, parameter and measurement noises can be specified by Q₀ + ε₁I_n, Q_p0 + ε₃I_m and R₀ + ε₂I respectively Define the estimation error vector $\tilde{z}\left(t\right)=z\left(t\right)-\widehat{z}\left(t\right)$ and the error covariance matrix $P\left(t\right)=E\left[\tilde{z}\left(t\right){\tilde{z}}^T\left(t\right)\right]$ . It can be proved that the error covariance matrix P and Kalman gain K satisfy, respectively, ²⁰

$$\dot{P}=GP+P{G}^T+L-P{H}^T{\left(R_0+{\epsilon }_2{I}_{n+m}\right)}^{-1}\hspace{0.17em}HP,\hspace{1em}{P}_0=E\left({\tilde{z}}_0{\tilde{z}}_0^T\right)$$ (26) $$K=P{H}^T{\left(R_0+{\epsilon }_2{I}_{n+m}\right)}^{-1}$$ (27)

where G and H are, respectively, the linearized system and output matrices given by

$$\begin{array}{c}G={\frac{\partial g\left(x\left(t\right),p\left(t\right)\right)}{\partial z\left(t\right)}|}_{x\left(t\right)=\widehat{x}\left(t\right),p\left(t\right)=\widehat{p}\left(t\right)}={\left[\begin{array}{cc}\frac{\partial f}{\partial x^T}& \frac{\partial f}{\partial p^T}\\ 0& 0\end{array}\right]|}_{x\left(t\right)=\widehat{x}\left(t\right),p\left(t\right)=\widehat{p}\left(t\right)}\\ H={\frac{\partial h\left(x\left(t\right),p\left(t\right)\right)}{\partial z\left(t\right)}|}_{x\left(t\right)=\widehat{x}\left(t\right),p\left(t\right)=\widehat{p}\left(t\right)}={\left[\begin{array}{cc}\frac{\partial h}{\partial x^T}& \frac{\partial h}{\partial p^T}\end{array}\right]|}_{x\left(t\right)=\widehat{x}\left(t\right),p\left(t\right)=\widehat{p}\left(t\right)}\end{array}$$

and

$$L=\left[\begin{array}{cc}{Q}_0+{\epsilon }^1{I}_m& 0_{m\times n}\\ 0_{m\times n}& Q_{p0}+{\epsilon }_3{I}_m\end{array}\right]$$

The linearization is performed around the estimated state $\widehat{x}\left(t\right)$ and parameter $\widehat{p}\left(t\right)$ respectively.

As it can be observed from (26) that the magnitude of the Kalman gain K is closely related to the amount of measurement noise reflected by the size of R₀ and the extent of the uncertain noise covariance specified by ε₂. The term L accounts for the increase of extrinsic noise and parameter noise and the term −PH^T(R₀ + ε₂I_n+m)^–1 HP reflects the decrease of uncertainty as a result of measurement.

When there is in the absence of extrinsic noise and a priori information in the biological system, by using the following matrix identity

$$\frac{d}{dt}{P}^{-1}=-P^{-1}\dot{P}{P}^{-1}$$

the continuous Riccati matrix equation (26) can be written in the linear equation in P⁻¹ as

$${\dot{P}}^{-1}=-P^{-1}G-G^T{P}^{-1}+H^T{\left(R_0+{\epsilon }_2{I}_{n+m}\right)}^{-1}\hspace{0.17em}H,\hspace{1em}{P}^{-1}\left(0\right)=P_0^{-1}$$ (28)

Clearly, either large measurement error (large R₀) or large uncertainty of the measurement covariance (large ε₂) cause the error covariance P to increase considerably whenever a measurement is utilized. This also results in a smaller Kalman gain for small state estimate errors and ease off the updating speed of state.

Performance Analysis

It follows from (23) and (24) that the estimation error dynamics is given by

$$\begin{array}{c}\dot{\tilde{z}}\left(t\right)=g\left(z\left(t\right)\right)+\eta \left(t\right)-g\left(\widehat{z}\left(t\right)\right)-Kh\left(z\left(t\right)\right)-Kv\left(t\right)+Kh\left(\widehat{z}\left(t\right)\right)\\ =\left(G-KH\right)\tilde{z}\left(t\right)+\hspace{0.17em}\Delta G\left(z,\widehat{z}\right)-K\Delta H\left(z,\widehat{z}\right)+\eta \left(t\right)-Kv\left(t\right)\end{array}$$ (29)

where ΔG(z,ẑ) = Δg(ẑ) − Δg(z) and ΔH(z,ẑ) = Δh(ẑ) − Δh(z) with Δg(·) = g(·) − G and Δh(·) = h(·) − H. It is reasonable to have

$$\begin{array}{c}\Vert \Delta G\left(z,\widehat{z}\right)\Vert \le {\rho }_1\Vert \tilde{z}\left(t\right)\Vert ,\forall \tilde{z}\left(t\right)\in {ℝ}^{n+m},\hspace{0.17em}t\ge 0\\ \Vert \Delta H\left(z,\widehat{z}\right)\Vert \le {\rho }_2\Vert \tilde{z}\left(t\right)\Vert ,\forall \tilde{z}\left(t\right)\in {ℝ}^{n+m},\hspace{0.17em}t\ge 0\end{array}$$

where ρ₁ and ρ₂ are finite positive constants and

$$\begin{array}{l}\Vert \eta \left(t\right)\Vert \le \Vert w_x\left(t\right)\Vert +\Vert w_p\left(t\right)\Vert \le \sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)},\hspace{1em}\forall t\ge 0,\hfill \\ \Vert v\left(t\right)\Vert \le \sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)},\hspace{1em}\forall t\ge 0,\hfill \end{array}$$ (30)

with tr(·) denoting the trace operation.

It can be observed from the estimation error dynamics (29) that its solution is

$$\tilde{z}\left(t\right)=\Phi \left(t\right){\tilde{z}}_0+{\int }_0^t\Phi \left(t-\tau \right)\left[\Delta G\left(z\left(\tau \right),\widehat{z}\left(\tau \right)\right)-K\Delta H\left(z\left(\tau \right),\widehat{z}\left(\tau \right)\right)\right]d\tau +{\int }_0^t\Phi \left(t-\tau \right)\left[\eta \left(\tau \right)-Kv\left(\tau \right)\right]d\tau$$ (31)

where Φ(t) = e^(G−KH)t. Or equivalently,

$$\tilde{z}\left(t\right)=\Phi \left(t\right){\tilde{z}}_0+{\int }_0^t\Phi \left(t-\tau \right)\left[I-K\right]\left[\begin{array}{c}\Delta G\left(z\left(\tau \right),\widehat{z}\left(\tau \right)\right)\\ \Delta H\left(z\left(\tau \right),\widehat{z}\left(\tau \right)\right)\end{array}\right]d\tau +{\int }_0^t\Phi \left(t-\tau \right)\left[I-K\right]\upsilon \left(\tau \right)d\tau$$

where $\upsilon \left(t\right)=\left[\begin{array}{c}\eta \left(t\right)\\ v\left(t\right)\end{array}\right]$ .

On the basis of the Kalman filtering theory, the matrix G − KH will be asymptotically stable if (G, H) is detectable around all estimated states ǩ(t). Therefore, there exist positive constants m, m₁ and β such that

$$\Vert \Phi \left(t\right)\Vert \le m{e}^{-\beta t},\forall t$$ (32)

and

$$\Vert \Phi \left(t\right)\left[I-K\right]\Vert \le m_1{e}^{-\beta t},\forall t$$ (33)

with the induced norms specified by Lemma 1, where β can be chosen to be $\underset{i}{\min }\left|\text{Re}\hspace{0.17em}{\lambda }_i\left(G-KH\right)\right|$ where

K = K(∊_1,∊_1,∊₁).

To proceed, by taking norms on both sides of (31) gives

$$\Vert \tilde{z}\left(t\right)\Vert \le \Vert \Phi \left(t\right)\Vert \Vert {\tilde{z}}_0\Vert +\left({\rho }_1+{\rho }_2\right){\int }_0^t\Vert \Phi \left(t-\tau \right)\left[I-K\right]\Vert \Vert \tilde{z}\left(\tau \right)\Vert d\tau +{\int }_0^t\Vert \Phi \left(t-\tau \right)\left[I-K\right]\Vert \left(\Vert \eta \left(\tau \right)\Vert +\Vert \nu \left(t\right)\Vert \right)d\tau$$ (34)

Or equivalently

$$\begin{array}{ll}\Vert \tilde{z}\left(t\right)\Vert e^{\beta t}\hfill & \le m\Vert {\tilde{z}}_0\Vert +m_1\left({\rho }_1+{\rho }_2\right){\int }_0^t\Vert \tilde{z}\left(\tau \right)\Vert e^{\beta \tau }d\tau +m_1{\int }_0^t\Vert \upsilon \left(\tau \right)\Vert e^{\beta \tau }d\tau \hfill \\ \hfill & \le m\Vert {\tilde{z}}_0\Vert +m_1\left({\rho }_1+{\rho }_2\right){\int }_0^t\Vert \tilde{z}\left(\tau \right)\Vert e^{\beta \tau }d\tau \hfill \\ \hfill & +m_1{\int }_0^t\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)\hspace{0.17em}{e}^{\beta \tau }d\tau \hfill \\ \hfill & =m\Vert {\tilde{z}}_0\Vert -\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)+\gamma \left(t\right)\hfill \\ \hfill & +m_1\left({\rho }_1+{\rho }_2\right){\int }_0^t\Vert \tilde{z}\left(\tau \right)\Vert e^{\beta \tau }d\tau \hfill \end{array}$$

where

$$\gamma \left(t\right)=\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)e^{\beta \tau }$$

Applying the Bellman-Gronwall inequality further yields

$$\Vert \tilde{z}\left(t\right)\Vert e^{\beta t}\le \left[m\Vert {\tilde{z}}_0\Vert -\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)\right]e^{m_1\left({\rho }_1+{\rho }_2\right)t}+\gamma \left(t\right)+m_1\left({\rho }_1+{\rho }_2\right){\int }_0^t\gamma \left(\tau \right)e^{m_1\left({\rho }_1+{\rho }_2\right)\left(t-\tau \right)}d\tau$$

Or

$$\Vert \tilde{z}\left(t\right)\Vert \le \left[m\Vert {\tilde{z}}_0\Vert -\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)\right]e^{-\left[\beta -m_1\left({\rho }_1+{\rho }_2\right)\right]t}+\gamma \left(t\right)e^{-\beta t}+m_1\left({\rho }_1+{\rho }_2\right)e^{-\left[\beta -m_1\left({\rho }_1+{\rho }_2\right)\right]t}{\int }_0^t\gamma \left(\tau \right)e^{m_1\left({\rho }_1+{\rho }_2\right)\left(t-\tau \right)}d\tau$$

Clearly, when

$$\beta >m_1\left({\rho }_1+{\rho }_2\right)$$ (35)

then

$$\Vert \tilde{z}\left(t\right)\Vert \le \frac{m_1}{\beta }(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}++m_1\left({\rho }_1+{\rho }_2\right)\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)e^{-\left[\beta -m_1\left({\rho }_1+{\rho }_2\right)\right]t}$$ $${\int }_0^t{e}^{\left[\beta -m_1\left({\rho }_1+{\rho }_2\right)\right]\tau }d\tau =\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)+m_1\left({\rho }_1+{\rho }_2\right)\frac{m_1}{\beta }\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)\frac{1}{\beta -m_1\left({\rho }_1+{\rho }_2\right)}$$

That implies that $\Vert \tilde{z}\left(t\right)\Vert$ wouldn't be diverged with its upper bound specified by

$$\underset{t}{\sup }\hspace{0.17em}\Vert \tilde{z}\left(t\right)\Vert \le \frac{m_1\left(\sqrt{n{\epsilon }_1+\text{tr}\left(Q_0\right)}+\sqrt{m{\epsilon }_3+\text{tr}\left(Q_{p0}\right)}+\sqrt{r{\epsilon }_2+\text{tr}\left(R_0\right)}\right)}{\beta -m_1\left({\rho }_1+{\rho }_2\right)}$$ (36)

Equations (32) and (35) specifies a crucial condition for determining convergence of the noisily corrupted system which is determined by the Kalman gain K and the values of ρ₁ and ρ₂. The condition is significantly affected by the amount of linearization and noise covariances. Inspection of the equations describing the behavior of the error covariance matrix reveals several observations which confirm engineers' intuition about the operation of the KF. As it can be observed from (36) that the estimation error is closely related to the upper bounds of the unreduced certain process and parameter noises and uncertain measurement noise. The larger the statistical parameters of the disturbances as reflected in the sizes of Q and Q _p , and the more pronounced effect of the disturbances as reflected in the size of R, the more rapidly the error covariance increases.

Larger Kalman gains will expedite the convergence of the estimation error. However, the estimation error increases considerably whenever there are larger noise uncertainties specified by larger ∊ _i, i = 1, 2, 3 and the linearization errors characterized by larger ρ₁ and ρ₂.

When there are process and measurement noises and uncertainties of noise covariances, from the observation of (26), (28) and (36), one should increase the magnitude of K for a larger stability margin β so as to assure convergence of the estimation error. However, as it was shown in (28), large measurement errors and the small error convariance P result in a small K. Thus, there is always a compromise between the optimal state estimation and stability robustness while designing the state estimator.

As for the implementation issue, possibility of the practical implementation of the estimator can be referred, for example, to, ¹⁸ which utilized the green fluorescent protein (GFP) as a reporter for real-time bioprocess sensing and GFP concentration and other important states in bioreactor culture of transgenic tobacco cells were successfully estimated. Application of the idea to the current estimator design deserves more attention and is worthy of further investigation.

Simulation Study

For demonstration, two examples for a class of noisy gene regulatory networks are illustrated.

Example 1 : Consider first a two-order system model for a real gene regulatory network given as follows: ⁶

$$\begin{array}{l}\left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\alpha }_1{x}_1^{g_{11}}\left(t\right)x_2^{g_{12}}\left(t\right)-{\beta }_1{x}_1^{h_{11}}\left(t\right)x_2^{h_{12}}\left(t\right)\\ {\alpha }_2{x}_1^{g_{21}}\left(t\right)x_2^{g_{22}}\left(t\right)-{\beta }_2{x}_1^{h_{21}}\left(t\right)x_2^{h_{22}}\left(t\right)\end{array}\right]+w_x\left(t\right),\hfill \\ {\left[\begin{array}{cc}{x}_{10}& x_{20}\end{array}\right]}^T={\left[\begin{array}{cc}1& 1.5\end{array}\right]}^T,\hfill \\ {\alpha }_i={\beta }_i=1,\hspace{1em}i-1,\hspace{0.17em}2,\hfill \\ g_{11}=0.268,\hspace{0.17em}{g}_{12}=-2.26,\hspace{0.17em}{g}_{21}=2.739,\hspace{0.17em}{g}_{22}=0.155,\hfill \\ h_{11}=0.465,\hspace{0.17em}{h}_{12}=0.359,\hspace{0.17em}{h}_{21}=0.197,\hspace{0.17em}{h}_{22}=0.281\hfill \end{array}$$

and

$$y=x_2+\nu \left(t\right)$$

where w_x(t) ~ (0, 0.1), w_p(t) ~ (0, 0.1) and ν(t) ~ (0, 0.2). States and parameters are both estimated. The parameters including 4 rate constants and 8 kinetic orders are treated as the states. Thus the state variables are extended from 2 to 14. In this network, for the first term of the first differential equation (the rate of change of x₁), ie, X^0.268₁ (t)x₂^−2.26 (t) with unit rate constant α₁, shows accumulation of gene product 1. Since the variable x₂ is raised to the power of the kinetic parameter −2.26 which reveals gene 2 will inhibit product of gene 1. On the other hand, for the second differential equation, the first term x₁^2.739(t) x₂^0.155(t) with unit rate constant α₂ reflects accumulation of gene product 2. Since the variable x₁ is raised to the power of 2.739 which reveals that gene product 1 will activate gene 2. The second terms −x₁^0.469 (t)x₂^0.359 (t) and −x₁^0.197 (t)x₂^0.281 (t) with unit rate constants β₁ and β₂ in the first and second differential equations reflecting degradation effect on gene products 1 and 2 respectively. Figure 1 illustrates the branch pathway of the two-dimensional S-system network.

The linearized system matrix ∂f/∂x^T based on the state estimation is

$${\frac{\partial f}{\partial x^T}|}_{x=\widehat{x}}=\left[\begin{array}{cc}0.268{\widehat{x}}_1^{-0.732}{\widehat{x}}_2^{-2.26}-0.469{\widehat{x}}_1^{-0.531}{\widehat{x}}_2^{0.359}& -2.26{\widehat{x}}_1^{-0.268}{\widehat{x}}_2^{-3.26}-0.359{\widehat{x}}_1^{0.46}{\widehat{x}}_2^{-0.641}\\ 0.739{\widehat{x}}_1^{1.739}{\widehat{x}}_2^{0.155}-0.197{\widehat{x}}_1^{-0.803}{\widehat{x}}_2^{0.281}& 0.155{\widehat{x}}_1^{2.739}{\widehat{x}}_2^{-0.845}-0.281{\widehat{x}}_1^{0.197}{\widehat{x}}_2^{-0.719}\end{array}\right],$$

and

$$H=\left[\begin{array}{ccc}0& 1& 0_{1\times 12}\end{array}\right]$$

Figure 1.

The gene regulatory network for Example 1.

For the initial error covariance P(0) = I₁₄, the results of dynamic simulation of the noise-free and estimated states and parameters for the noisy gene regulatory network using the proposed robust EKF given by (26) and (27) are shown as in Figure 2.

Figure 2.

Dynamic simulation of the noise-free and estimated gene states and parameters with P(0) = I₁₄; A) gene states; y-axis is concentration, x-axis is time,B–D) parameters.

Consider next the existence of uncertainties of the extrinsic noise and measurement noise with ∊₁ = 0.05, ∊₃ = 0.05 and ∊₂ = 0.1. The state and parameter responses are shown in Figure 3. The root mean square error (RMSE) was used to quantify the filtering performance with

$$RMSE=\sqrt{\frac{1}{T_2-T_1}{\int }_{T_1}^{T_2}{\left({\widehat{x}}_i\left(t\right)-x_{fi}\left(t\right)\right)}^2\hspace{0.17em}dt}$$

where T₁ ≤ t ≤ T₂, ${\widehat{x}}_i\left(t\right)$ (t) is the i-th estimation state and x_fi(t) is the corresponding noise-free state. The RMSE values for the noise-free and estimated states and parameters in 0 ≤ t ≤ 5 are listed in Table 1 which shows that the estimator is able to filter the extrinsic and measurement noises to retrieve the real state and parameter values.

Figure 3.

Dynamic simulation of the noise-free and estimated gene states and parameters with noise uncertainties (∊ = 0.05, ∊₃ = 0.05 and ∊₂ = 0.1);A) gene states; y-axis is concentration, x-axis is time,B–D) parameters.

Table 1.

Comparison of RMSE values for the system of Example 1 with and without noise uncertainties.

State	with noise uncertainties (∊1 = ∊3 = 0.05, ∊2 = 0.1)	without noise uncertainties
1	0.065	0.043
2	0.053	0.035
3	0.04	0.028
4	0.03	0.017
5	0.032	0.021
6	0.034	0.021
7	0.018	0.01
8	0.022	0.014
9	0.04	0.027
10	0.032	0.021
11	0.032	0.016
12	0.026	0.017
13	0.009	0.007
14	0.05	0.03

Example 2 : Consider a nonlinear gene regulatory of four genes described and shown in Figure 4: ¹³

$$\begin{array}{l}\left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\\ {\dot{x}}_3\left(t\right)\\ {\dot{x}}_4\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\alpha }_1{x}_2^{g_{12}}\left(t\right)-{\beta }_1{x}_1^{h_{11}}\left(t\right)\\ {\alpha }_2{x}_3^{g_{23}}\left(t\right)x_4^{g_{24}}\left(t\right)-{\beta }_2{x}_2^{h_{22}}\left(t\right)\\ {\alpha }_3{x}_1^{g_{31}}\left(t\right)x_2^{g_{32}}\left(t\right)-{\beta }_3{x}_3^{h_{33}}\left(t\right)\\ {\alpha }_4{x}_3^{g_{43}}\left(t\right)-{\beta }_4{x}_4^{h_{44}}\left(t\right)\end{array}\right]+w_x\left(t\right),\hfill \\ {\left[\begin{array}{cccc}{x}_{10}& x_{20}& x_{30}& x_{40}\end{array}\right]}^T=0.5{\left[\begin{array}{cccc}1& 1& 1& 1\end{array}\right]}^T\hfill \\ {\alpha }_1=1,\hspace{0.17em}{\alpha }_2=2,\hspace{0.17em}{\alpha }_3=1.5,\hspace{0.17em}{\alpha }_4=3.5,\hspace{0.17em}{\beta }_1=1,\hspace{0.17em}{\beta }_2=2,\hspace{0.17em}{\beta }_3=1.5,\hspace{0.17em}{\beta }_4=3.5,\hspace{0.17em}\hfill \\ g_{12}=1,\hspace{0.17em}{g}_{23}=-0.5,\hspace{0.17em}{g}_{24}=2,\hspace{0.17em}{g}_{31}=4,\hspace{0.17em}{g}_{32}=1,\hspace{0.17em}{g}_{43}=-0.5,\hfill \\ h_{11}=1,\hspace{0.17em}{h}_{22}=1,\hspace{0.17em}{h}_{33}=2,\hspace{0.17em}{h}_{44}=2\hfill \end{array}$$

where w_x ~ (0, 0.2) and w_p ~ (0, 0.2). The state variables are extended from 4 to 22 when 18 parameters including rate constants and kinetic orders are all treated as the state variables. The measurement model is given as

$$y\left(t\right)=x_2\left(t\right)+x_3\left(t\right)+\nu \left(t\right)$$

where y(t) is the measurement output and the measurement noise ν(t) ~ N(0, 0.5). For this gene regulatory network, gene product 1 activates gene 3, gene product 2 activates genes 1 and 3, gene product 3 represses genes 2 and 4, and gene product 4 activates gene 2.

Figure 4.

Another example of the gene regulatory network.

The linearized system matrix ∂f/∂x^T based on the state estimation for the EKF design can be obtained as

$${\frac{\partial f}{\partial x^T}|}_{x=\widehat{x}}=\left[\begin{array}{cccc}-1& 1& 0& 0\\ 0& -2& -{\widehat{x}}_3^{-1.5}\left(t\right){\widehat{x}}_4^2\left(t\right)& 4{\widehat{x}}_3^{-0.5}\left(t\right){\widehat{x}}_4\left(t\right)\\ 6{\widehat{x}}_1^3\left(t\right){\widehat{x}}_2\left(t\right)& 1.5{\widehat{x}}_{14}^4\left(t\right)& -3{\widehat{x}}_3\left(t\right)& 0\\ 0& 0& -17.5{\widehat{x}}_3^{-1.5}\left(t\right)& -7{\widehat{x}}_4\left(t\right)\end{array}\right],$$

and

$$\text{H}=\left[\begin{array}{cccc}0& 1& 1& 0_{1\hspace{0.17em}\times \hspace{0.17em}19}\end{array}\right]$$

For the initial covariance matrix P(0) = I₂₂, the results of dynamic simulation of the noise-free gene states and the estimated states of the noisy gene network are shown as in Figure 5. As it can be seen that the estimator tracked the noise-free case well while there were extrinsic noise and measurement noise.

Figure 5.

Dynamic simulation of the noise-free and estimated gene states and parameters with P(0) = I₂₂; A) gene states, y-concentration, x-axis is time,B–E) parameters.

Consider next the existence of uncertainties of the extrinsic noise and measurement noise with ∊₁ = 0.1, ∊₃ = 0.1 and ∊₂ = 0.25. The gene responses are shown in Figure 6. As in the previous example, the results exhibits larger estimation errors due to added noise uncertainties, however, the deviation is not significant while compared with magnitudes of the nominal state or parameter responses.

Figure 6.

Dynamic simulation of the noise-free and estimated gene states and parameters with noise uncertainties (∊₁ = 0.1, ∊₃ = 0.1 and ∊₂ = 0.25);A) gene states, y-axis is concentration, x-axis is time,B–E) parameters.

Conclusions

This paper proposes a continuous EKF to estimate internal states and parameters of a class of gene networks while there are extrinsic and intrinsic noises and parametric fluctuations. Quantitative performance analysis for state estimation of the EKF is presented. Numerical simulations have confirmed possibility of the proposed method in designing robust EKFs. This shows potential of the presented design method in bridging the engineering approach to solve for the estimation problem in biological systems.