title: Stochastic filters based on hybrid approximations of multiscale stochastic reaction networks
authors: Zhou Fang Ankit Gupta Mustafa Khammash
text: 
  
  
  
  Efficiency gain and bidirectional operation of quantum engines with decoupled internal levels
      Daniel Jonathan
      December 30, 2023
  =============================================================================================
  
  empty
  empty
  
  
  
  
  
  We consider the problem of estimating the dynamic latent states of an intracellular multiscale stochastic reaction network from time-course measurements of fluorescent reporters. We first prove that accurate solutions to the filtering problem can be constructed by solving the filtering problem for a reduced model that represents the dynamics as a hybrid process. The model reduction is based on exploiting the time-scale separations in the original network, and it can greatly reduce the computational effort required to simulate the dynamics. This enables us to develop efficient particle filters to solve the filtering problem for the original model by applying particle filters to the reduced model. We illustrate the accuracy and the computational efficiency of our approach using a numerical example.
  
  
  
  
  
  § INTRODUCTION
  
  
  In the past few decades, scientists' ability to look into the dynamic behaviors of a living cell has been greatly improved by the fast development of fluorescent technologies and advances in microscopic techniques.
  Despite this big success, light intensity signals observed in a microscope can only report the dynamics of a small number of key components in a cell, such as fluorescent proteins and mRNAs, and, therefore, leave other dynamic states of interest, e.g., gene (on/off) state, indirectedly observable.
  As a result, it is urgent to build efficient stochastic filters to estimate latent states of intracellular biochemical reaction systems from these partial observations.
  
  Although it is possible to solve the filtering (and inference) problem accurately for intracellular systems with particle filters (see <cit.> for particle filters), the required computational effort to simulate the system in the sampling step usually prevents it from being used on-line due to the high complexity of practical models.
  Motivated by this obstacle, some methods are proposed to obtain computationally efficient filters by approximating the underlying Markov jump process by a more tractable model, such as a diffusion process <cit.> or a block model <cit.>.
  All these methods are shown to be efficient for broad classes of bio-chemical reaction systems, however, their drawbacks are also clear — the diffusion process strategy loses its validity in low copy number scenarios, and the block model methodology may not work well for highly interconnected nonlinear biological circuits.
  Unfortunately, some reaction systems involve species with low copy numbers and densely interconnected biological circuits, which precludes the utilization of previous methods to these systems and require researchers to build new solutions. 
  
  In this paper, we propose another strategy to obtain efficient particle filters based on hybrid approximations of multiscale stochastic reaction networks.
  Specifically, the model reduction technique introduced in <cit.> is applied to approximating the underlying stochastic reaction system by a piecewise deterministic model, which can greatly reduce the computational effort required to simulate the system.
  Also, we  prove that the solution to the filtering problem for the original system can be constructed by solving the filtering problem for the reduced model if some mild conditions are satisfied. 
  Both of these facts enable us to develop efficient particle filters to solve the filtering problem for the original model by applying particle filters to the reduced model. 
  
  
  It is worth noting that the idea of applying time-scale separation techniques to reduce underlying models to make particle filters computationally feasible has already been proposed in <cit.> for diffusion type stochastic models, and its efficiency is proven in the literature <cit.>.
  Compared with these works, our paper considers a different type of underlying model, the Markov jump process, and provides a much-simplified proof for the convergence of the approximate filters.
  All these references and our paper show the efficiency of applying the time-scale separation technique to the filtering problem for multiscale systems. 
  
  The rest of this paper is organized as follows.
  In Section <ref>, we introduce some basic notations and review the chemical reaction network theory and the filtering theory.
  In Section <ref>, we present our main results. 
  We first prove that the accurate solution to the filtering problem can be accurately approximated by the solution to the filtering problem for the reduced model under some mild condition.
  Then, based on this result, we establish an efficient particle filter to solve the filtering problem for the original model by applying the particle filter to the reduced model.
  For the sake of readability, we put the proof of our main result in the appendix.
  In Section <ref>, a numerical example is presented to illustrate our approach.
  Finally, Section <ref> concludes this paper.
  
  
  
  
  
  
  § PRELIMINARY
   
  
  
  
   §.§ Notations
  
  We first denote a natural filtered probability space by ( Ω,  ℱ, {ℱ_t}_t≥ 0, ℙ), where Ω is the sample space, ℱ is the σ-algebra, {ℱ_t}_t≥ 0 is the filtration, and ℙ is the natural probability.
  Also, we term ℕ_>0 as the set of positive integers, ℝ^n with n being a positive integer as the space of n-dimensional real vectors, · as the Euclidean norm, |·| as the absolute value notation, s∧ t (for any s,t ∈ℝ ) as min(s,t), and s∨ t as max(s,t).
  For any positive integer n and any t>0, we term 𝔻_ℝ^ n[0,t] as the Skorokhod space that consists of all ℝ^n valued cadlag functions on [0,t] and 𝔻_ℝ^ n[0,∞) as the Skorokhod space that consists of all ℝ^n valued cadlag functions on [0,∞).
  
  
  
  
   §.§ Stochastic chemical reaction networks
   
  
  We consider an intracellular system undergoing reactions
  
      v_1,j S_1 + … + v_n,j S_n  ->[k_j] v'_1,j S_1 + … + v'_n,j S_n,    j=1,…, r,
  
  where S_i (i=1,…,n) are distinguished species in the systems, v_i,j and v'_i,j are non-negative integers, called stoichiometric coefficients, k_j are the reaction constants, and r is the number of reactions.
  Also, we name the linear combination of species (e.g, v_1,j S_1 + … + v_n,j S_n).
  The topology of the above reaction equations can be fully represented by a triplet {𝒮, 𝒞, ℛ}, called a chemical reaction network <cit.>, where
  
  	
    * 𝒮 is the species set {S_1, S_2, …, S_n} whose elements represent distinct substances,
  	
    * 𝒞 is the complex set ⋃_j=1^r{v_· j , v'_· j}, where v_· j≜(v_1,j, …, v_n,j)^⊤ and v'_· j≜(v'_1,j, …, v'_n,j)^⊤, indicating substrate complexes and product complexes,
  	
    * ℛ is the reaction set {v_·,1→ v'_·1, …, v_· r→ v'_· r} showing the connections between the complexes.
  
  Let X(t)=(X_1(t),X_2(t), …, X_n(t))^⊤ be the numbers of molecules of these species at time t, then the system's dynamics following mass-action kinetics can be expressed as <cit.>
  
      X(t)= X(0)+∑_j=1^r (v'_· j-v_· j) R_j( ∫_0^t λ_j (X(s)) s )
  
  where R_j(t) (j=1,…,r) are mutually independent unit rate Poisson random processes, λ_j (x)≜ k_j∏_i=1^n x_i!/(x_i-v_i,j)! 1_{x_i≥ v_i,j} with 1_{·} the indicator function, and the initial condition X(0) is subject to a particular known distribution.
  
  
  
  
   §.§ Model reduction via the time scale separation
  
  In a practical bio-chemical reaction system, different species can vary a lot in abundance, and rate constants can also vary over several orders of magnitude. 
  To normalize these quantities, we term N as a scaling factor, α_i (for i=1,…, n) as the magnitude of the i-th species such that X^N_i(t)≜ N^-α_iX_i(t) = O(1), β_j (for j=1,…,r) as the magnitude of the reaction constant k_j such that k'_j≜ k_j N^-β_j=O(1),
  γ as the time scale parameter, and X^N,γ(t) ≜ X^N(tN^γ).
  Then, in the new coordinate, the dynamic equation can be expressed as <cit.>
  
      X^N,γ(t) =    X^N,γ(0)   +∑_j=1^rΛ^N (v'_· j-v_· j) 
                   × R_j( ∫_0^t N^γ+ρ̃_jλ^N_j(X^N,γ(s))  s )
  
  where Λ^N= diag(N^-α_1,…,N^-α_s), α=(α_1,…,α_n), λ^N_j (x)≜ k'_j ∏_i=1^n∏_ℓ=0^v_i,j-1  (x_i-ℓ N^-α_i) 1_{x_i≥ N^-α_iv_i,j} are of order O(1), and ρ̃_j=β_j + v^⊤_· jα.
  For this system, we denote its transition kernel by K^N,γ_1_t(x,A)≜ℙ ( X^N,γ_1(t)∈ A | X^N,γ_1(0)=x ).
  
  Following the notations in <cit.>, we further term Γ_i^+≜{ j | v'_i,j>v_i,j} as the set of reactions that result in an increase of the i-th species, Γ_i^-≜{ j | v'_i,j<v_i,j} as the set of reactions that result in a decrease of the i-th species.
  Then, the constant α_i-max_j∈Γ_i^+∪Γ_i^-(β_j+v_· j^⊤α) is the time scale of the i-th substance, i.e., the time scale where X^N,γ_i (t) evolves at the rate of O(1).
  Moreover, we term γ_1≜min_i
  ( α_i-max_j∈Γ_i^+∪Γ_i^-(β_j+v_· j^⊤α)) as the parameter of the fastest time scale of the system and
  D^α̃≜diag(… 1_{α_i=α̃}…) as a diagonal matrix indicating whether a species is at the scale of N^α̃.
  In the rest of this paper, we focus on the dynamic behaviors of a reaction system at the fastest time scale γ_1.
  
  For a large scaling factor N, by neglecting all slow reactions (γ_1+ρ̃_j<0) and approximating fast reactions (γ_1+ρ̃_j>0) by a continuous process, we can arrive at a simplified piecewise deterministic process as follows.
  
  
      X^γ_1(t)=   lim_N→∞ X^N(0) 
         
      + ∑_j : γ_1+ρ̃_j>0∫_0^t λ'_j( X^γ_1(s))  D^γ_1+ρ̃_j(v'_· j-v_· j) s  
          +  ∑_j : γ_1+ρ̃_j=0 
      R_j(∫_0^t λ'_j( X^γ_1(s))  s ) D^0(v'_· j-v_· j).
  
  We denote the transition kernel of the above process by K^γ_1_t(x,A)≜ℙ ( X^γ_1(t)∈ A | X^γ_1(0)=x ).
  If we further assume the initial conditions to satisfy that 
  
      lim_N→∞ X^N(0)  exists,
  
  and that the above systems are non-explosive, i.e.,
  
      lim_c→∞τ^N,γ_1_c=∞,     and   lim_c→∞τ^γ_1_c=∞,     a.s.,
  
  where τ^N,γ_1_c≜inf{t | X^N,γ_1(t)≥ c } and τ^γ_1_c≜inf{t | X^γ_1(t)≥ c},
  then the original model X^N,γ_1(·) converges in distribution to the reduced model X^γ(·) on any finite time interval. 
  
  
  	On any finite time interval [0,t], if (<ref>) and (<ref>) hold, then X^N,γ_1(· )⇒ X^γ_1(·) in the sense of the Skorokhod topology.
  
  
  	It follows from <cit.> and the portmanteau lemma.
  
  
  Also, the finite-dimensional distribution of X^N,γ_1(·) converges to the finite-dimensional distribution of X^γ_1(·).
  
  
  	Let i be a positive integer and {t_j}_0<j≤ i be a sequence of scalars. If (<ref>) and (<ref>) hold, then X^N,γ_1_t_1:t_i⇒ X^γ_1_t_1:t_i where X^N,γ_1_t_1:t_i≜(X^N,γ_1(t_1), …, X^N,γ_1(t_i)) and X^γ_1_t_1:t_i≜(X^γ_1(t_1), …, X^γ_1(t_i)).
  
  
  	It follows immediately from Proposition <ref> and <cit.>.
  
  
  Note that the convergence result does not always hold on the infinite time horizon, because the limit model may miss bi-stability or some other phenomenon that the full model has (see <cit.>).
  
  In the above results, the non-explosivity condition (<ref>) is not a very strong condition; there exist checkable sufficient conditions to ensure non-explosivity (see <cit.>) that can deal with it efficiently.
  In this paper, we do not provide specific conditions for non-explosivity in order to make the result general and avoid the deviation from this paper's main concern.
  
  
  	The computational complexity to simulate X^γ_1(·) can be greatly lower than the complexity to simulate X^N, γ_1(·), because the former avoids the exact simulation of fast reactions (γ+ρ̃_j >0), which consume a lot of computational resources to update the system at a rate proportional to N^γ_1+ρ̃_j.
  
  
  
  
   §.§ Filtering problems, the change of measure method, and the particle filter.
   
  
  In this paper, we assume that m channels of light intensity signals of fluorescent reporters in a single cell can be observed discretely in time from a microscopy platform (e.g., the one in <cit.>).
  The discrete-time property is mainly caused by the cell segmentation and tracking, which is required prior to the measuring step and can be time-consuming. 
  Suppose that {t_i}_i∈ℕ_>0 are time points when the observation comes, and the value of each observation, Y^N, γ_1(t_i), satisfies 
  
      Y^N,γ_1_ℓ(t_i)= h_ℓ(X^N,γ_1(t_i))  + W_ℓ(t_i)     ℓ=1,…,m,
  
  where h_ℓ are bounded Lipschitz continuous functions indicating the relation between the observation and the reaction process, and W_ℓ(t_i) (ℓ =1,…,m and i∈ℕ_>0) are mutually independent standard Gaussian random variables which are also independent of R_j(·) (j=1,…, r).
  
  Our task is to provide optimal estimates to latent states of the system, e.g., ϕ(X^N,γ_1(t_i)) where ϕ is a known measurable function, in the sense of mean square error based on observations up to the time t_i; in other words, we are going to calculate the conditional expectation π^N,γ_1_t_i(ϕ) ≜𝔼 _ℙ[ ϕ(X^N,γ_1(t_i))  | 𝒴^N,γ_1_t_i], where 𝒴^N,γ_1_t_i is the filtration generated by the process {Y^N,γ_1(t_j)}_0 < j≤ t_i.
  
  Similarly, we can define artificial readouts for the reduced model by 
  
      Y^γ_1_ℓ(t_i)= h_ℓ(X^γ_1(t_i))  + W_ℓ(t_i)     ℓ=1,…,m,
  
  and attempt to solve the corresponding filtering problem, i.e., calculating the conditional expectation π^γ_1_t_i(ϕ) ≜𝔼 _ℙ[ ϕ(X^γ_1(t_i))  | 𝒴^γ_1_t_i], where 𝒴^γ_1_t_i is the filtration generated by the process {Y^γ_1(t_j)}_0 < j≤ t_i.
  Note that in practical experiments, one can only observe {Y^N,γ_1_ℓ(t_i)}_i∈ ℕ_>0 but never get {Y^γ_1_ℓ(t_i)}_i∈ ℕ_>0.
  
  The change of measure method is a powerful tool to deal with filtering problems, whose core idea lies in constructing a reference probability measure that orthogonalizes the underlying process and the observation.
  For our problem, we construct an auxiliary random process Z^N,γ_1(t_i)= ∏_j=1^i g(X^N,γ_1(t_j),Y^N,γ_1(t_j)), where g(x,y)=exp(∑_ℓ=1^m h_ℓ(x) y_ℓ-h^2_ℓ(x)/2), and whose reciprocal is a martingale with respect to the filtration ℱ_t_i under ℙ.
  Based on this martingale, we define a reference probability measure, ℙ^N,γ_1, by the Radon-Nikodym derivation . ℙ^N,γ_1/ℙ|_ℱ_t_i= (Z^N,γ_1(t_i))^-1, under which {Y^N,γ_1(t_i)}_i∈ℕ_>0 are mutually independent standard m-dimensional Gaussian random variables, the law of the process X^N,γ_1(·) is the same as its law under ℙ, and the process X^N,γ_1(·) and observation Y^N,γ_1(·) are independent of each other. (These results can be easily shown by checking the joint characteristic function of these processes.)
  Furthermore, by the Kallianpur-Striebel formula <cit.>, for any measurable function ϕ such that ϕ(X^N,γ_1(t_i)) is integrable under ℙ, there holds
  
      π^N,γ_1_t_i(ϕ) = 𝔼_ℙ^N,γ_1[. Z^N,γ_1(t_i) ϕ(X^N,γ_1(t_i)) |  𝒴^N,γ_1_t_i]/𝔼 _ℙ^N,γ_1[ . Z^N,γ_1(t_i) |  𝒴^N,γ_1_t_i]
  
  ℙ-a.s. and ℙ^N,γ_1-a.s..
  Recall that π^N,γ_1_t_i(ϕ) is the notation for 𝔼 _ℙ[ ϕ(X^N,γ_1(t_i))  | 𝒴^N,γ_1_t_i].
  Here, we name ρ^N,γ_1_t_i (ϕ)≜𝔼_ℙ^N,γ_1[. Z^N,γ_1(t_i) ϕ(X^N,γ_1(t_i)) |  𝒴^N,γ_1_t_i] as the unnormalized conditional expectation and ρ^N,γ_1_t_i (1) as the normalization factor.
  Note that both π^N,γ_1_t_i(ϕ) and ρ^N,γ_1_t_i(ϕ) are 𝒴^N,γ_1_t_i-measurable.
  Therefore, for any bounded ϕ(·), there exist measurable functions, f̂^N,γ_1_ϕ,t_i(·) and ĝ^N,γ_1_ϕ,t_i(·) from the domain of the observations to ℝ such that f̂^N,γ_1_ϕ,t_i(Y^N,γ_1_t_1:t_i)=π^N,γ_1_t_i(ϕ) and ĝ^N,γ_1_ϕ,t_i(Y^N,γ_1_t_1:t_i)=ρ^N,γ_1_t_i(ϕ), where Y^N,γ_1_t_1:t_i≜(Y^N,γ_1(t_1),…,Y^N,γ_1(t_i)) is the observation process up to the time t_i.
  
  Based on (<ref>), an algorithm called the particle filter (also known as sequential Monte Carlo method) can be constructed as Algorithm <ref> to numerically solve the filtering problem. 
  In Algorithm <ref>, the particles x_j(·) mimic the dynamic state X^N,γ_1(·), and the weights w_j(·) mimic the normalized Z^N,γ_1(·); therefore, the computed filter can approximate the true filter by (<ref>).
  A resampling step is executed in each iteration to remove non-significant particles so that sample impoverishment can be avoided <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Similarly, for the filtering problem of the reduced model, we can define another auxiliary function Z^γ_1(t_i)= ∏_j=1^i g(X^γ_1(t_j),Y^γ_1(t_j)) and a reference probability ℙ^γ_1 by . ℙ^γ_1/ℙ|_ℱ_t_i= (Z^γ_1(t_i))^-1, under which {Y^γ_1(t_i)}_i∈ℕ_>0 are mutually independent standard m-dimensional Gaussian random variables, the law of the process X^γ_1(·) is the same as its law under ℙ, and the process X^γ_1(·) and observation Y^γ_1(·) are independent of each other.
  By the Kallianpur-Striebel formula, for any measurable function ϕ such that ϕ(X^γ_1(t_i)) is integrable under ℙ, there holds
  
      π^γ_1_t_i(ϕ) = 𝔼_ℙ^γ_1[. Z^γ_1(t_i) ϕ(X^γ_1(t_i)) |  𝒴^γ_1_t_i]/𝔼 _ℙ^γ_1[ . Z^γ_1(t_i) |  𝒴^γ_1_t_i]
  
  ℙ-a.s. and ℙ^γ_1-a.s., in which we further denote the numerator by ρ^γ_1_t_i (ϕ) and the denominator by ρ^γ_1_t_i (1).
  Note that both π^γ_1_t_i(ϕ) and ρ^γ_1_t_i(ϕ) are 𝒴^γ_1_t_i-measurable.
  Therefore, for any bounded ϕ(·), there exist measurable functions, f̂^γ_1_ϕ,t_i(·) and ĝ^γ_1_ϕ,t_i(·) from the domain of the observations to ℝ, such that f̂^γ_1_ϕ,t_i(Y^γ_1_t_1:t_i)=π^γ_1_t_i(ϕ) and ĝ^γ_1_ϕ,t_i(Y^γ_1_t_1:t_i)=ρ^γ_1_t_i(ϕ), where Y^γ_1_t_1:t_i≜(Y^γ_1(t_1),…,Y^γ_1(t_i)) is the observation process up to the time t_i.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § MAIN RESULTS
   
  
  Recall that our task is to accurately and computationally efficiently solve the filtering problem for the multi-scale system (<ref>), i.e., to calculate {π^N,γ_1_t_i(ϕ)}_i∈ℕ_>0.
  A straightforward idea to solve this problem is to use the particle filter where the transition kernel of the full model (i.e., K^N,γ_1(·,·)) and observations Y^N,γ_1(·) are inserted. 
  We denote this particle filter by π̅^N,γ_1_M,t_i(ϕ).
  Although it can accurately solve the filtering problem if the resampling method is properly chosen (see <cit.>), this algorithm involves the simulation of the full model (<ref>) and, therefore, can be computationally expensive (see Remark <ref>).
  Consequently, we need to figure out a smarter way to approach this problem.
  
  We first show that the exact solution to the filtering problem of the original model can be constructed by solving the filtering problem for the reduced model. 
  
  
  	Assume that (<ref>) and (<ref>) are satisfied.
  	Then, for any i∈ℕ_>0 and any bounded continuous function ϕ(·), 
  	
  	there holds  π^N,γ_1_t_i(ϕ)-π̃^N,γ_1_t_i(ϕ) ≜f̂_ϕ,t_i^N,γ_1(Y^N,γ_1_t_1:t_i) - f̂_ϕ,t_i^γ_1(Y^N,γ_1_t_1:t_i) → 0, as N goes to infinity, in ℙ-probability.
  
  
  	The proof is in the appendix.
  
  
  Recall that f̂_ϕ,t_i^γ_1(·) is a mapping that maps the observation to the solution of the filtering problem for the reduced model.  
  The above theorem suggests that the probability of the error of these filters below a given threshold is very close to 1 if the scaling factor, N, is large, and, therefore, that π̃^N,γ_1_t_i (ϕ)≜f̂_ϕ,t_i^γ_1(Y^N,γ_1_t_1:t_i) is a good approximation to the exact filter π^N,γ_1_t_i (ϕ).
  
  
  
  Let us constructed a particle filter π̅^γ_1_M,t_i(ϕ) where the transition kernel of the reduced model (i.e., K^γ_1(·,·)) and observations Y^N,γ_1(·) are inserted.
  Clearly, this particle filter is an approximation of the filter f̂_ϕ,t_i^γ_1(Y^N,γ_1_t_1:t_i).
  Based on Theorem <ref>, we can further show that this particle filter π̅^γ_1_M,t_i(ϕ) can also accurately approximate the true filter π^N,γ_1_t_i(ϕ) with high probability.
  
  
  
  
  
  
  
  
  	Assume that (<ref>) and (<ref>) are satisfied, and the resampling step in the particle filter (Algorithm <ref>) adopts the multinomial branching method.
  	Then for any i∈ℕ_>0, any positive δ, and any bounded continuous function ϕ(·), there holds
  	
      lim_N→∞lim_M→∞ℙ(|π̅^γ_1_M,t_i(ϕ)-π^N,γ_1_t_i(ϕ)| >δ) =0.
  
  
  
  	It is a direct consequence of Corollary 2.4.4 in the literature <cit.> and Theorem <ref>.
  
  
  
  Recall that the required computational effort to simulate the reduced model X^γ(·) is much lower than the effort to simulate X^N,γ(·) (see Remark <ref>).
  As a consequence, particle filter π̅^γ_1_M,t_i(ϕ), which applies the kernel function K^γ_1(·,·) and, therefore, is required to simulate the reduced model X^γ_1(·) in the sampling step, has a much lower computational cost compared with the particle filter π̅^N,γ_1_M,t_i(ϕ), which applies kernel function K^N, γ_1(·,·) and is required to simulated the full model.
  This fact, together with the accuracy result Corollary <ref>, indicates that the filter π̅^γ_1_M,t_i(ϕ) based on the reduced model is a better candidate than π̅^N,γ_1_M,t_i(ϕ) for solving a practical filtering problem for multiscale intracellular reaction systems.
  
  
  
  
  We need to emphasize that the boundedness condition of ϕ is not generally satisfied in biological applications, where concentrations of most target species have no theoretical upper bounds apart from some exceptions, e.g., gene copies. To tackle this problem, we can truncate a target quantity by a large number beyond which the conditional probability is comparatively low, and the tail event contributes little to the conditional expectation. Then, an estimation of the truncated quantity generated by the particle filtering applying the reduced model can provide an accurate approximation to the conditional expectation of the target quantity.
  
  
  
  
  
  
  § NUMERICAL EXAMPLE
  
  
  In this section, we illustrate our approach using a numerical example.
  We consider a gene expression model involving four species: a gene having on and off states (denoted respectively as S_2 and S_1), the mRNA it transcribes (denoted as S_3), and a fluorescent protein it expresses (denoted as S_4), and six reactions: S_1 ->[k_1] S_2 (gene activation), S_2 ->[k_2] S_1  (gene deactivation), S_2 ->[k_3] S_3 + S_2 (gene transcription), S_3	 ->[k_4]  S_3+S_4 (translation), S_3 ->[k_5]∅ (mRNA degradation), and S_4  ->[k_6]∅ (protein degradation).
  
  
      S_1    ->[k_1] S_2       Gene activation 
      
      S_2    ->[k_2] S_1       Gene deactivation
      
      S_2    ->[k_3] S_3 + S_2       Gene transcription, 
      
      S_3	    ->[k_4]  S_3+S_4        Translation
      
      S_3    ->[k_5]∅      mRNA degradation,
      
      S_4     ->[k_6]∅      Protein degradation.
  
  
  
  For this system, we take N=100, α_4=1, and α_i=0 for 1≤ i ≤ 3, i.e., the cellular system consists of hundreds of fluorescent protein molecules but very few copies of other molecules. 
  The values of reaction constants and their scaling exponents are shown in Table <ref>, where the unit of both reaction constants and scale rates is “minute^-1".
  Also, we assume X^N_1(0) to have a binary distribution with mean 1/3, X^N_2(0) to satisfy X^N_2(0)=1- X^N_1(0), X^N_3(0) to have a Poisson distribution with mean 2, and X^N_4(0) to also have a Poisson distribution with mean 2.
  Moreover, we assume that all initial conditions except X_2(0) are independent of each other.
  
  With this setting, we can easily calculate that γ_1=0, i.e., the fastest time scale is 0, and the full dynamic model satisfies
  
  
      X^N,γ_1(t) 
         = X^N,γ_1(0) + 
      	ζ_1
      	R_1 ( k'_1  ∫_0^t X^N,γ_1_1 (s)  s ) 
            + 
          ζ_2
      	R_2(k'_2∫_0^t X^N,γ_1_2(s)  s)
      	+		ζ_3
      	R_3(k'_3∫_0^t X^N,γ_1_2(s)  s)
           +
      	ζ_4N^-1R_4(k'_4  N ∫_0^t X^N,γ_1_3(s)  s)
          +
      	ζ_5 R_5(k'_5∫_0^t X^N,γ_1_3(s)  s) 
            +
      	ζ_6 N^-1R_6(k'_6 N  ∫_0^t X^N,γ_1_4(s)  s)
  
  where ζ_1=(-1,1,0,0)^⊤, ζ_2=(1,-1,0,0)^⊤, ζ_3=(0,0,1,0)^⊤, ζ_4=(0,0,0,1)^⊤, ζ_5=(0,0,-1,0)^⊤, and ζ_6=(0,0,0,-1)^⊤;
  the reduced dynamics satisfies
  
  
      X^γ_1(t) 
         = lim_N→∞X^N,γ_1(0)+ 
      		R_1 ( k'_1  ∫_0^t X^N,γ_1_1 (s)  s ) 
            + 
      	ζ_2
      	R_2(k'_2∫_0^t X^N,γ_1_2(s)  s)
      	+		ζ_3
      	R_3(k'_3∫_0^t X^N,γ_1_2(s)  s)
            + ζ_4∫_0^t k'_4 X^N,γ_1_3(s)  s  +
      	ζ_5 R_5(k'_5∫_0^t X^N,γ_1_3(s)  s) 
            
      	 +
      	ζ_6 ∫_0^t k'_6  X^N,γ_1_4(s)  s.
  
  We simulate the full dynamics by a modified next reaction method <cit.> and the reduced model by Algorithm 2 in <cit.>.
  
  
  
  
  
  
      X^N,γ_1_1(t)=    X^N,γ_1_1(0)- R_1 ( k'_1  ∫_0^t X^N,γ_1_1 (s)  s )  
         + R_2(k'_2∫_0^t X^N,γ_1_2(s)  s)
      
      X^N,γ_1_2(t)=    X^N,γ_1_2(0)+ R_1 ( k'_1  ∫_0^t X^N,γ_1_1 (s)  s ) 
         
      - R_2(k'_2∫_0^t X^N,γ_1_2(s)  s)  
      
      X^N,γ_1_3(t)=   X^N,γ_1_3(0)+R_3(k'_3∫_0^t X^N,γ_1_2(s)  s)  
         -R_5(k'_5∫_0^t X^N,γ_1_3(s)  s) 
       
      X^N,γ_1_4(t)=    X^N,γ_1_4(0) + N^-1R_4(k'_4  N ∫_0^t X^N,γ_1_3(s)  s)
           - N^-1R_6(k'_6 N  ∫_0^t X^N,γ_1_4(s)  s)
  
  and the reduced dynamics satisfies
  
      X^γ_1_1(t)=   lim_N→∞X^N,γ_1_1(0)- R_1 ( k'_1  ∫_0^t X^γ_1_1 (s)  s ) 
         + R_2(k'_2∫_0^t X^γ_1_2(s)  s)
      
      X^γ_1_2(t)=   lim_N→∞ X^N,γ_1_2(0)+ R_1 ( k'_1  ∫_0^t X^γ_1_1 (s)  s )  
         - R_2(k'_2∫_0^t X^γ_1_2(s)  s)  
      
      X^γ_1_3(t)=   lim_N→∞ X^N,γ_1_3(0) 
      +R_3(k'_3∫_0^t X^γ_1_2(s)  s) 
         -R_5(k'_5∫_0^t X^γ_1_3(s)  s) 
      
      X^γ_1_4(t)=   lim_N→∞ X^N,γ_1_4(0) 
       +   ∫_0^t k'_4 X^N,γ_1_3(s)  s 
         - ∫_0^t k'_6  X^N,γ_1_4(s)  s.
  
  
  
  
  We assume that light intensity signals can be observed every 2 minutes from a fluorescent microscope 
  and satisfy
  
      Y^N,γ_1(t)= h(X^N,γ_1(t)) + W(t)
  
  where h(x)= (10 × x_4 ) ∧ 10^3 with 10^3 being the measurement range, and W(t) is a sequence of mutually independent standard Gaussian random variables.
  Our goal is to provide estimates of latent dynamic states, including DNA copies, mRNA copies, and protein concentrations, based on readouts.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  A numerical simulation is shown in Fig <ref>, where we first simulated the full dynamic system and observations, which were taken as true signals, and then applied two different types of particle filters to estimate the system. 
  Here, we set the particle population to be 5000.
  From the results, we can observe that both filters can follow the behavior of the true signal, and the particle filter applying the reduced model is within the credible interval (light red region) most of the time.
  In the meanwhile, the particle filter applying the original model took 9.44 seconds in average to output the estimation in each iteration, whereas the particle filter  applying the reduced model took only 0.95 seconds in each iteration, 10 times faster than its opponent.
  These results indicate that the proposed particle filter applying the reduced model is both accurate and computationally efficient in solving filtering problems for multi-scale stochastic chemical reaction network systems. 
  
  
  
  
  § CONCLUSION
  
  
  In this paper, we provide efficient particle filters to solve the filtering problem for multiscale stochastic reaction network systems, based on a time-scale separation technique.
  We first show that the solution of the filtering problem for the original systems can be well approximated by the solution of the filtering problem for a reduced model that represents the dynamics as a hybrid process. 
  Since the reduced model is based on exploiting the time-scale separation in the original network, the reduced model can greatly reduce the computational complexity to simulate the dynamics.
  This enables us to develop efficient particle filters to solve the filtering problem for the original model by applying the particle filter to the reduced model.
  Finally, a numerical example is presented to illustrate our approach.
  
  There are a few topics deserving further investigation in future work.
  First, since the variability of model parameters can greatly influence the performance of a filter, it is worth investigating if we can put parameter inference and latent state estimation into the same framework to improve the performance of the proposed particle filter.
  Second, it is worth extending the convergence result of the particle filter (see Corollary <ref>) to unbounded function scenarios, as many biological states are not necessarily bounded.
  Third, we will also apply the proposed filter to the optimal control problem for an intracellular system.
  
  
  
  
  In this section, we give the outline of the proof of theorem <ref>.
  Here, we utilize a framework proposed in <cit.>, which require us to construct random variables X̃^N,γ_1_t_1:t_i, X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i) and Z̃^γ_1(t_i) on a common probability space (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) at each time point t_i 
  such that 
  
  	
    * random variables ( X̃^N,γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i)) take values in ℝ^n× i×ℝ^m× i×ℝ and have the same law as (  X^N,γ_1_t_1:t_i,  Y^N,γ_1_t_1:t_i, Z^N,γ_1(t_i)) on ( Ω, ℱ,  ℙ^N,γ_1). 
  	
    * the random variable ( X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^γ_1(t_i)) takes values in ℝ^n× i×ℝ^m× i×ℝ and has the same law as (  X^γ_1_t_1:t_i,  Y^γ_1_t_1:t_i, Z^γ_1(t_i)) on ( Ω, ℱ,  ℙ^γ_1). 
  	
    * f̂_ϕ,t^N,γ_1(Ỹ^γ_1_t_1:t_i) →f̂_ϕ,t^γ_1(Ỹ^γ_1_t_1:t_i) in ℚ̃^γ_1-probability. 
  	
    * lim_N→∞𝔼_ℚ̃^γ_1[|Z̃^N,γ_1(t_i)-Z̃^γ_1(t_i)|]=0. 
  
  If we succeed in finding the above random variables, then the convergence result is guaranteed by the following theorem.
   
  	For any i> 0 and any measurable function ϕ(·), if <ref>, <ref>, <ref>, and <ref> are satisfied,
  	then π^N,γ_1_t_i (ϕ) - π̃^N,γ_1_t_i (ϕ)→ 0 in ℙ-probability.
  
  
  
  
  	
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We can construct such random variables and the common probability space as follows.
  First, for any i∈ℕ_>0, we can apply the Skorokhod representation theorem to constructing  a common probability space (Ω̃_1^γ_1, ℱ̃_1^γ_1, ℚ̃_1^γ_1) on which are defined ℝ^n× i-valued random variables X̃^N,γ_1_t_1:t_i and X̃^γ_1_t_1:t_i with the laws of X^N,γ_1_t_1:t_i and X^γ_1_t_1:t_i, respectively, such that lim_N→∞X̃^N,γ_1_t_1:t_i= X̃^γ_1_t_1:t_i ℚ̃_1^γ_1-almost surely.
  Then, we term (Ω̃_2^γ_1, ℱ̃_2^γ_1, ℚ̃_2^γ_1) as a copy of the probability space ( Ω, ℱ_t_i, ℙ^γ_1) and Ỹ^γ_1_t_1:t_i as a copy of Y^γ_1_t_1:t_i in this probability space.
  Moreover, we define Z̃^Nγ_1(t_i)≜∏_j=1^i g(X̃^Nγ_1(t_j),Ỹ^γ_1(t_j)) and Z̃^γ_1(t_i)≜∏_j=1^i g(X̃^γ_1(t_j),Ỹ^γ_1(t_j)) on the product probability space (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) ≜(Ω̃_1^γ_1×Ω̃_2^γ_1, ℱ̃_1^γ_1×ℱ̃_2^γ_1,ℚ̃_1^γ_1×ℚ̃_2^γ_1), 
  where  X̃^N,γ_1(t_j) is the j-th column of X̃^N,γ_1_t_1:t_j, X̃^γ_1(t_j) is the j-th column of X̃^γ_1_t_1:t_j, and Ỹ^γ_1(t_j) is the j-th column of Ỹ^γ_1_t_1:t_j.
  
  Obviously, such defined random variables and the common probability space satisfy <ref> and <ref>.
  The condition <ref> can be proven straightforwardly by the fact that lim_N→∞X̃^N,γ_1_t_1:t_i= X̃^γ_1_t_1:t_i ℚ̃_1^γ_1-almost surely.
  Similarly, we can also show <ref> by using (<ref>) and (<ref>). 
  Consequently, Theorem <ref> is proven.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Let (Ω̃_2^γ_1, ℱ̃_2^γ_1, ℚ̃_2^γ_1) be a copy of the probability space ( Ω, ℱ_t_i, ℙ^γ_1), and Ỹ^γ_1_t_1:t_i be a copy of Y^γ_1_t_1:t_i in this probability space.
  Therefore, columns of Ỹ^γ_1_t_1:t_i are mutually independent standard m-dimensional Gaussian random variables.
  Then, we further define random variables Z̃^Nγ_1(t_i)≜∏_j=1^i g(X̃^Nγ_1(t_j),Ỹ^γ_1(t_j)) and Z̃^γ_1(t_i)≜∏_j=1^i g(X̃^γ_1(t_j),Ỹ^γ_1(t_j)) in the product probability space (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) ≜(Ω̃_1^γ_1×Ω̃_2^γ_1, ℱ̃_1^γ_1×ℱ̃_2^γ_1,ℚ̃_1^γ_1×ℚ̃_2^γ_1), 
  where  X̃^N,γ_1(t_j) is the j-th column of X̃^N,γ_1_t_1:t_j, X̃^γ_1(t_j) is the j-th column of X̃^γ_1_t_1:t_j, and Ỹ^γ_1(t_j) is the j-th column of Ỹ^γ_1_t_1:t_j.
  Finally, we can show these random variables and the product probability space to satisfy <ref> and <ref>.
  
  
  	Assume that (<ref>) and (<ref>) hold, then the above-defined probability space (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) and random variables X̃^N,γ_1_t_1:t_i, X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i) and Z̃^γ_1(t_i)  satisfy  <ref> and <ref>.
  
  
  	It follows from Lemma <ref>, the fact that X^N,γ_1(·) (X^γ_1(·)) and Y^N,γ_1(·) (Y^γ_1(·)) are independent under ℙ^N,γ_1 (ℙ^γ_1), and the fact that the mappings (  X̃^N,γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i)→Z̃^N,γ_1(t_i), (  X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i)→Z̃^γ_1(t_i), (   X^N,γ_1_t_1:t_i, Y^N,γ_1_t_1:t_i) →  Z^N,γ_1(t_i) and (   X^γ_1_t_1:t_i, Y^γ_1_t_1:t_i) →  Z^γ_1(t_i) are measurable and identical.
  
  
  
  
  
  
  
  
  	
  
  
  
  
  
   §.§ Convergence result
  
  This subsection contributes to showing the above constructed random variables to satisfy <ref> and <ref> and, therefore,  π^N,γ_1_t_i (ϕ) - π̃^N,γ_1_t_i (ϕ)→ 0 in ℙ-probability. First, we work on <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	Assume that (<ref>) and (<ref>) are satisfied, then the above-defined probability space (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) and random variables X̃^N,γ_1_t_1:t_i, X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i) and Z̃^γ_1(t_i) satisfy lim_N→∞𝔼_ℚ̃^γ_1[|Z̃^N,γ_1(t_i)- 	Z̃^γ_1(t_i) |^q]= 0 for any q≠ 0
  	and, therefore, <ref>.
  
  
  
  	By the Cauchy-Schwarz inequality, we have
  	
      𝔼_ℚ̃^γ_1[|Z̃^N,γ_1(t_i)- 	Z̃^γ_1(t_i) |^q]≤
         √(𝔼_ℚ̃^γ_1[(Z̃^N,γ_1(t_i)/Z̃^γ_1(t_i) - 	1)^2q]𝔼_ℚ̃^γ_1[(Z̃^γ_1(t_i))^2q]).
  
  	Note that 𝔼_ℙ^γ_1[(Z^γ_1(t_i))^2q] is bounded due to the boundedness of function h(·).
  	Then, we can also conclude that 𝔼_ℚ̃^γ_1[(Z̃^γ_1(t_i))^2q] is finite by the boundedness of  𝔼_ℙ^γ_1[(Z^γ_1(t_i))^2q] and the fact that Z^γ_1(t_i) (under ℙ^γ_1) and Z̃^γ_1(t_i) have the same law (c.f. Lemma <ref>).
  	Therefore, to prove the result, we only need to show 𝔼_ℚ̃^γ_1[(Z̃^N,γ_1(t_i)/Z̃^γ_1(t_i) - 	1)^2q] to converge to zero.
  	By simple calculations, we have that 
  	
      𝔼_ℚ̃^γ_1[(Z̃^N,γ_1(t_i)/Z̃^γ_1(t_i) - 	1)^2q]
         =𝔼_ℚ̃^γ_1[𝔼_ℚ̃^γ_1[(Z̃^N,γ_1(t_i)/Z̃^γ_1(t_i) - 	1)^2q| 𝒳̃^N,γ_t_i⋁𝒳̃^γ_t_i.]] 
         =  ∑_j'=0^2 C_2q^j'(-1)^j'𝔼_ℚ̃^γ_1[ 
      	A_j']
  
  	where 𝒳̃^N, γ_1_t_i and 𝒳̃^γ_1_t_i are  σ-algebras generated by X̃^N,γ_1_t_1:t_i and X̃^γ_1_t_1:t_i, respectively, C_2q^j' are binomial coefficients, and 
  	
      A_j'=
         exp( j'/2∑_j=1^i∑_ℓ^m[h^2_ℓ(X̃^γ_1(t_j))-h^2_ℓ(X̃^N,γ_1(t_j))]  ) 
          ×exp( j'/2∑_j=1^i∑_ℓ^m[h_ℓ(X̃^N,γ_1(t_j))-h_ℓ(X̃^γ_1(t_j))]^2 ).
  
  	The last line in (<ref>) follows from the Gaussianity of Ỹ^γ_1_t_1:t_i under ℚ̃^γ_1.
  	Finally, 𝔼_ℚ̃^γ_1[(Z̃^N,γ_1(t_i)/Z̃^γ_1(t_i) - 	1)^2q] → 0 follows immediately from (<ref>), Lemma  <ref>, and the dominant convergence theorem, which proves the result.
  
  
  Then, we work on <ref> using the following lemma.
  
  
  	Let {A_N}_N, {B_N}_N, A, and B be random variables on a probability space.
  	If A_N→ A and B_N → B in probability, 
  	then A_NB_N→ AB in probability.
  
  
  
  
  
  
  
  
  
  	
  	
  	
  	 
  	
  	
  	
  	
  	
  	
  	
  
  
  
  
  
  
  
  
  
  Specifically, we use the above lemma to prove <ref> by taking A_N, A, B_N and B as ĝ^N,γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i), ĝ^γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i), 1/ĝ^N,γ_1_1,t_i  (Ỹ^γ_1_t_1:t_i) and 1/ĝ^γ_1_1,t_i  (Ỹ^γ_1_t_1:t_i), respectively.
  We first show that A_N→ A and B_N→ B in probability.
  
  
  	Assume that (<ref>) and (<ref>) are satisfied, and X̃^N,γ_1_t_1:t_i, X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i), Z̃^γ_1(t_i), and (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) are defined as in the previous subsection.
  	Then, for any bounded continuous function ϕ(·),
  	
  	there holds ĝ^N,γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i) →ĝ^γ_1_ϕ,t_i  (Ỹ^γ_1_t_1:t_i) in ℚ̃^γ_1 probability.
  	Moreover, by the continuous mapping theorem, there also holds 1/ĝ^N,γ_1_1,t_i  (Ỹ^γ_1_t_1:t_i)→1/ĝ^γ_1_1,t_i  (Ỹ^γ_1_t_1:t_i)  in ℚ̃^γ_1 probability.
  
  
  
  	
  	Note that ( X̃^N,γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i)) have the same laws as (  X^N,γ_1_t_1:t_i,  Y^N,γ_1_t_1:t_i, Z^N,γ_1(t_i)) in ( Ω, ℱ,  ℙ^N,γ_1) (c.f. Lemma <ref> and <ref>). 
  	Then, by the fact that ρ^N,γ_1_t_i(ϕ)=ĝ^N,γ_1_ϕ,t_i  ( Y^N,γ_1_t_1:t_i) holds ℙ^N,γ_1-a.s., we have that ĝ^N,γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i)=𝔼_ℚ̃^γ_1[Z̃^N,γ_1(t_i)ϕ(X̃^N,γ_1(t_i))| 𝒴̃^γ_1_t_i.] holds ℚ̃^γ_1-almost surely, where 𝒴̃^γ_1_t_i is the minimum σ-algebra generated by Ỹ^γ_1_t_1:t_i.
  	Similarly, we can show that ĝ^γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i)=𝔼_ℚ̃^γ_1[Z̃^N,γ_1(t_i)ϕ(X̃^γ_1(t_i))| 𝒴̃^γ_1_t_i.] holds ℚ̃^γ_1-almost surely.
  	Therefore, by Jensen's inequality and the triangle inequality, we can conclude
  	
      𝔼_ℚ̃^γ_1[
      	| ĝ^N,γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i)-ĝ^γ_1_ϕ,t_i (Ỹ^γ_1_t_1:t_i)
      	|
      	]
         ≤𝔼_ℚ̃^γ_1[
      	|(Z̃^N,γ_1(t_i)-Z̃^γ_1(t_i))ϕ(X̃^γ_1(t_i))|
      	]
             +𝔼_ℚ̃^γ_1[
      	|Z̃^N,γ_1(t_i)(ϕ(X̃^N,γ_1(t_i))-ϕ(X̃^γ_1(t_i)))|
      	]
  
  	
  	The first term on the right hand side of the inequality converges to zero by the Cauchy-Schwarz inequality and Lemma <ref>. 
  	The second term also converges to zero by the Cauchy-Schwarz inequality, Lemma <ref> and the dominant convergence theorem. 
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	Therefore, the result is shown.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	
  	
  
  
  
  
  	
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	
  	
  
  
  
  
  
  Finally, we can show <ref> and Theorem <ref>.
  
  
  	  Assume that (<ref>) and (<ref>) are satisfied, and X̃^N,γ_1_t_1:t_i, X̃^γ_1_t_1:t_i, Ỹ^γ_1_t_1:t_i, Z̃^N,γ_1(t_i), Z̃^γ_1(t_i), and (Ω̃^γ_1, ℱ̃^γ_1, ℚ̃^γ_1) are defined as in the previous subsection.
  	Then, for any bounded continuous function ϕ(·),
  	
  	the condition <ref> holds.
  
  
  
      By (<ref>) and the fact that Ỹ^γ_1 and Y^N,γ_1 (under ℙ^N,γ_1) have the same law, we have that 
     f̂^N,γ_1_ϕ,t(Ỹ^γ_1) = ĝ^N,γ_1_ϕ,t(Ỹ^γ_1)/ĝ^N,γ_1_1,t(Ỹ^γ_1)  ℚ̃^γ_1-almost surely.
  	Similarly, we can also conclude f̂^γ_1_ϕ,t(Ỹ^γ_1) = ĝ^γ_1_ϕ,t(Ỹ^γ_1)/ĝ^γ_1_1,t(Ỹ^γ_1) ℚ̃^γ_1-almost surely.
  	Then, by applying <ref> and <ref> to these relations, we show the lemma.
  
  
  The proof of Theorem <ref>):
  	The theorem follows immediately from Lemma <ref>, Lemma <ref>, Lemma <ref>, and Theorem <ref>.
  ▪
  
  
  
  
  
  
  
  10
  	
  	doucet2009tutorial
  	A. Doucet and A. M. Johansen, “A tutorial on particle filtering and smoothing:
  	Fifteen years later,” Handbook of nonlinear filtering, vol. 12,
  	no. 656-704, p. 3, 2009.
  	
  	chuang2013robust
  	C.-H. Chuang, C.-L. Lin, et al., “Robust estimation of stochastic
  	gene-network systems,” Journal of Biomedical Science and Engineering,
  	vol. 6, no. 02, p. 213, 2013.
  	
  	golightly2006bayesian
  	A. Golightly and D. J. Wilkinson, “Bayesian sequential inference for
  	stochastic kinetic biochemical network models,” Journal of
  		Computational Biology, vol. 13, no. 3, pp. 838–851, 2006.
  	
  	calderazzo2019filtering
  	S. Calderazzo, M. Brancaccio, and B. Finkenstädt, “Filtering and inference
  	for stochastic oscillators with distributed delays,” Bioinformatics,
  	vol. 35, no. 8, pp. 1380–1387, 2019.
  	
  	boys2008bayesian
  	R. J. Boys, D. J. Wilkinson, and T. B. Kirkwood, “Bayesian inference for a
  	discretely observed stochastic kinetic model,” Statistics and
  		Computing, vol. 18, no. 2, pp. 125–135, 2008.
  	
  	kang2013separation
  	H.-W. Kang, T. G. Kurtz, et al., “Separation of time-scales and model
  	reduction for stochastic reaction networks,” The Annals of Applied
  		Probability, vol. 23, no. 2, pp. 529–583, 2013.
  	
  	park2008problem
  	J. H. Park, N. S. Namachchivaya, and R. B. Sowers, “A problem in stochastic
  	averaging of nonlinear filters,” Stochastics and Dynamics, vol. 8,
  	no. 03, pp. 543–560, 2008.
  	
  	park2010dimensional
  	J. H. Park, R. B. Sowers, and N. S. Namachchivaya, “Dimensional reduction in
  	nonlinear filtering,” Nonlinearity, vol. 23, no. 2, p. 305, 2010.
  	
  	park2011particle
  	J. H. Park, N. Namachchivaya, and H. C. Yeong, “Particle filters in a
  	multiscale environment: homogenized hybrid particle filter,” Journal of
  		Applied Mechanics, vol. 78, no. 6, 2011.
  	
  	imkeller2013dimensional
  	P. Imkeller, N. S. Namachchivaya, N. Perkowski, H. C. Yeong, et al.,
  	“Dimensional reduction in nonlinear filtering: a homogenization approach,”
  	The Annals of Applied Probability, vol. 23, no. 6, pp. 2290–2326,
  	2013.
  	
  	feinberg2019foundations
  	M. Feinberg, Foundations of Chemical Reaction Network Theory.
  	Springer, 2019.
  	
  	anderson2011continuous
  	D. F. Anderson and T. G. Kurtz, “Continuous time markov chain models for
  	chemical reaction networks,” in Design and analysis of biomolecular
  		circuits, pp. 3–42, Springer, 2011.
  	
  	ethier1986markov
  	S. N. Ethier and T. G. Kurtz, “Markov processes: Characterization and
  	convergence,” 1986.
  	
  	gillespie2000chemical
  	D. T. Gillespie, “The chemical langevin equation,” The Journal of
  		Chemical Physics, vol. 113, no. 1, pp. 297–306, 2000.
  	
  	gupta2014scalable
  	A. Gupta, C. Briat, and M. Khammash, “A scalable computational framework for
  	establishing long-term behavior of stochastic reaction networks,” PLoS
  		Comput Biol, vol. 10, no. 6, p. e1003669, 2014.
  	
  	anderson2018non
  	D. F. Anderson, D. Cappelletti, M. Koyama, and T. G. Kurtz, “Non-explosivity
  	of stochastically modeled reaction networks that are complex balanced,” Bulletin of mathematical biology, vol. 80, no. 10, pp. 2561–2579, 2018.
  	
  	rullan2018optogenetic
  	M. Rullan, D. Benzinger, G. W. Schmidt, A. Milias-Argeitis, and M. Khammash,
  	“An optogenetic platform for real-time, single-cell interrogation of
  	stochastic transcriptional regulation,” Molecular cell, vol. 70,
  	no. 4, pp. 745–756, 2018.
  	
  	kallianpur1968estimation
  	G. Kallianpur and C. Striebel, “Estimation of stochastic systems: Arbitrary
  	system process with additive white noise observation errors,” The
  		Annals of Mathematical Statistics, vol. 39, no. 3, pp. 785–801, 1968.
  	
  	crisan2001particle
  	D. Crisan, “Particle filters—a theoretical perspective,” in Sequential
  		Monte Carlo methods in practice, pp. 17–41, Springer, 2001.
  	
  	anderson2007modified
  	D. F. Anderson, “A modified next reaction method for simulating chemical
  	systems with time dependent propensities and delays,” The Journal of
  		chemical physics, vol. 127, no. 21, p. 214107, 2007.
  	
  	duncan2016hybrid
  	A. Duncan, R. Erban, and K. Zygalakis, “Hybrid framework for the simulation of
  	stochastic chemical kinetics,” Journal of Computational Physics,
  	vol. 326, pp. 398–419, 2016.
  	
  	calzolari2006approximation
  	A. Calzolari, P. Florchinger, and G. Nappo, “Approximation of nonlinear
  	filters for markov systems with delayed observations,” SIAM journal on
  		control and optimization, vol. 45, no. 2, pp. 599–633, 2006.
  	
  
  
  
  
  