On the Properties of the Reaction Counts Chemical Master Equation

Abstract

The reaction counts chemical master equation (CME) is a high-dimensional variant of the classical population counts CME. In the reaction counts CME setting, we count the reactions which have fired over time rather than monitoring the population state over time. Since a reaction either fires or not, the reaction counts CME transitions are only forward stepping. Typically there are more reactions in a system than species, this results in the reaction counts CME being higher in dimension, but simpler in dynamics. In this work, we revisit the reaction counts CME framework and its key theoretical results. Then we will extend the theory by exploiting the reactions counts’ forward stepping feature, by decomposing the state space into independent continuous-time Markov chains (CTMC). We extend the reaction counts CME theory to derive analytical forms and estimates for the CTMC decomposition of the CME. This new theory gives new insights into solving hitting times-, rare events-, and a priori domain construction problems.

1. Introduction

Continuous-time Markov chains (CTMC) accurately capture the dynamics of a broad range of biochemical reaction systems. The CTMCs come in two flavours: The discrete state space, the chain transitions from one state to another, or the continuous state space, the state transitions are smooth yet non-differentiable. Each type has an intuitive interpretation. Discrete state spaces describe the population or counts of chemicals in a system, while continuous states spaces describe concentrations of chemicals in a system [1,2]. It is tempting to generalise concentrations to be simply scaled populations, which would mean that all systems could be mapped into the continuous setting and solved there. However, this is not the case mathematically. It was shown that for large populations, the discrete CTMC and continuous CTMC were equivalent [3,4]. If the large population assumption is not satisfied, mathematically, the discrete state space system has to be studied in its own right.

The probability distribution of a CTMC over the discrete state space in the biochemical context is found by solving the chemical master equation (CME). We first present the CME and then describe its components. The formula of the CME is written as follows,

$$\frac{\partial p\left(Z\right(t)=\mathbf{z})}{\partial t}=\sum _{n=1}^{N_r}{a}_n(\mathbf{z}-{\nu }_n)p(Z\left(t\right)=\mathbf{z}-{\nu }_n)-\sum _{n=1}^{N_r}{a}_n\left(\mathbf{z}\right)p(Z\left(t\right)=\mathbf{z}).$$

(1)

In the equation above, $Z\left(t\right)$ is a CTMC over the state space $\mathrm{\Omega }\subset {\mathbb{N}}_0^{N_s}.$ Each state is a vector of $N_s\in \mathbb{N}$ non-negative integers describing the population of the species at time $t\ge 0.$ Starting from the left-hand side, the equation states that the change in probability of being in a state $z\in \mathrm{\Omega }$ at time t is given by the sum of two terms. The first term describes the probability flowing into the state $\mathbf{z}$ and the second term describes the probability flowing out of the state $\mathbf{z}.$ There are $N_r\in \mathbb{N}$ reactions perturbing the population; the terms ${\nu }_n$ are referred to as stoichiometric vectors, they capture the net change in population if the nth reaction fires. Hence, subtracting the stoichiometric vectors ${\nu }_n$ from $\mathbf{z}$ gives us the set of all previous states which might transition into $\mathbf{z}$ and “contribute” probability. The last piece of information in the equation are the rates at which reactions fire, these are encapsulated in the propensity functions $\{a_n:\mathrm{\Omega }\to {\mathbb{R}}_{+},\mathrm{for}n=1,\dots ,N_r\}.$ In summary, the CME describes the change in probability of observing a CTMC in a certain state, which is equal to the sum of probability transitioning into the state, minus the sum of probability leaving the state.

Even though the dynamics of the derivative are fairly simple, solving the CME is a hard problem [5,6,7,8,9,10,11,12]. A special case of the CME which is truly unyielding to approximation is when a system is close to its population boundary (for example, close to zero). In this scenario, none of the state transitions around a state near the boundary can be summarised by the average transitions out of the state. This impedes most approximation methods, and even if a method was successful, its computational effort would be equivalent to that of solving the Finite State Projection for the same error. An approach which was proposed in various different contexts to simplify dynamics, was to study the number of reactions fired rather than the population counts [6,13,14,15,16,17]. Because reactions cannot “un-fire”, the state transitions are only moving forward. Furthermore, given that the system’s starting population is known, the species count in any state can be reconstructed from the reaction counts.

In this work, we will reintroduce the reaction counts variant of the CME. We will then explore the theoretical and structural results which emerge from its forward stepping process. Then, we will show how the reaction counts CME can be used to construct solutions to the classical CME. Lastly, we decompose the state space of the CME into independent CTMCs and give analytical forms and estimates to calculate their probabilities. The purpose of this work is to gain intuition and explore structural properties of the CME, hence, the motivation is more theory rather than the application of the CME. In light of that, we omit a discussion of numerical methods in this paper and envisage this aspect in future research.

2. Reaction Counts CME

2.1. Formulation

We denote the state space of the reaction counts CME by $\mathrm{\Lambda }\subset {\mathbb{N}}_0^{N_r}.$ Each state is a vector of non-negative integers of length $N_r.$ Each element of the vector represents the number of times its corresponding reaction has fired. In the CME setting, the change in state after a reaction fires is given by the stoichiometric vector, which we denoted earlier by $v_n,$ for $n=\{1,\dots ,N_r\}.$ This vector quantifies the net change in populations after the reaction has fired. Analogously, in the reaction counts CME setting, a change in state indicates that a reaction has fired, and its corresponding reaction count is incremented by one. Therefore, the stoichiometric vector for the nth reaction is simply the identity vector ${\mathbf{1}}_n;$ the vector is zero in all but the nth position, where it is one.

We bridge the species counts CME and the reaction counts CME with the mapping $\mathrm{\Gamma }:\mathrm{\Omega }\times \mathrm{\Lambda }\to \mathrm{\Omega },$ that is, for $(x_0,r)\in \mathrm{\Omega }\times \mathrm{\Lambda },$

$$\begin{array}{ccc}\hfill \mathrm{\Gamma }(x_0,r)& :=& x_0+[v_1,v_2,\dots ,v_{N_r}]r,\hfill \\ & =& x_0+\sum _{n=1}^{N_r}{v}_n{r}_n.\hfill \end{array}$$

(2)

The mapping above links the reaction counts CME to the species counts CME by stating that, given the starting state $x_0$ and the reactions which have fired are known, then the current state in the species counts is the starting species state plus the sum of stochiometries of all the reactions which have fired. It is important to note that the map $\mathrm{\Gamma }$ is injective into the species state space $\mathrm{\Omega }.$ In most cases, $\mathrm{\Lambda }$ is of higher dimension than $\mathrm{\Omega },$ hence, $\mathrm{\Gamma }$ is seldom bijective. For our purpose, we only need the pull back map of $\mathrm{\Gamma }.$ We denote the pull back map as ${\mathrm{\Gamma }}^{-1}:\mathrm{\Omega }\times \mathrm{\Omega }\to \mathcal{P}\left(\mathrm{\Lambda }\right),$ where for $(x_0,x)\in \mathrm{\Omega }\times \mathrm{\Omega },$

$${\mathrm{\Gamma }}^{-1}(x_0,x):=\{r\in \mathrm{\Lambda }:\mathrm{\Gamma }(x_0,r)=x\}.$$

(3)

With the mapping between $\mathrm{\Lambda }$ and $\mathrm{\Omega }$ established, we can now inherit the propensity function from the species counts over to the reaction counts. For $n=1,\dots ,N_r,$ the reaction counts propensity of the nth reaction is given by,

$${\alpha }_n(x_0,r):=\left\{\begin{array}{cc}{a}_n\left(\mathrm{\Gamma }(x_0,r)\right)\hfill & \mathrm{if}\mathrm{\Gamma }(x_0,r)\in \mathrm{\Omega },\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(4)

With the state space, stoichiometry, and the propensities established in the reaction counts setting, we are ready to formulate the reaction counts based Kurtz process. For $x_0\in \mathrm{\Omega },$ let ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0}$ be the reaction counts Kurtz process,

$$R_{x_0}\left(t\right):=\sum _{n=1}^{N_r}{\mathbf{1}}_n\mathcal{P}\left({\int }_0^t{\alpha }_n(x_0,R_{x_0}\left(s\right))ds\right).$$

(5)

The corresponding reaction counts CME for the process above is given by,

$$\frac{dp(R_{x_0}\left(t\right)=r)}{dt}=\sum _{n=1}^{N_r}{\alpha }_n(x_0,r-{\mathbf{1}}_n)p(R_{x_0}\left(t\right)=r-{\mathbf{1}}_n)-\sum _{n=1}^{N_r}{\alpha }_n(x_0,r)p(R_{x_0}\left(t\right)=r),$$

(6)

for all $r\in \mathrm{\Lambda }.$ The initial condition for the reaction counts CME is a points mass on the origin,

$$p(R_{x_0}\left(0\right)=r)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}r=\mathbf{0},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(7)

Upon first glance, the CME of the reaction counts is simpler in its complexity compared to its species counter part. In the reaction case, the processes only move forward, which gives rise to simple forward propagating dynamics. Furthermore, given solutions to the reaction counts CME exist, with the simple application of the push forward measure, we would obtain:

$$p(X\left(t\right)=x)=\sum _{r\in {\mathrm{\Gamma }}^{-1}(X\left(0\right),x)}p(R_{X\left(0\right)}\left(t\right)=r).$$

(8)

The relationship given above is the critical motivation for studying the reaction counts CME. In principle, if structures and results could be attained in the reaction counts setting, then by using $\mathrm{\Gamma },$ these results can be mapped into species setting. The simplest example of this is the proof for the existence of solutions of the CME for finite time. We now show that the reaction counts CME has analytical solutions, then using the push forward measure in Equation (8), we prove that the solutions of the CME exist for a broad range of propensity functions.

2.2. Analytical Solutions of the Reaction Counts CME

Firstly, we will establish the notations and assumptions needed to present the results. For brevity, we denote the sum of all propensities ${\alpha }_n$ as:

$$\alpha (·,r):=\sum _{n=1}^{N_r}{\alpha }_n(·,r).$$

(9)

The reaction counts CME in Equation (6) is a difference differential equation, hence, by arranging the probability states in a vector $u\left(t\right):={\left(p(R_{x_0}\left(t\right)=r)\right)}_{r\in \mathrm{\Lambda }},$ we reduce solving the reaction counts CME to solving the ODE,

$$\frac{du\left(t\right)}{dt}=A^{*}u\left(t\right).$$

For $r,l\in \mathrm{\Lambda },$ the matrix $A^{*}$ has the properties: ${\left[A^{*}\right]}_{r,r}\le 0,$${\left[A^{*}\right]}_{r,l}\ge 0,$ with ${\sum }_{l\in \mathrm{\Lambda }}{\left[A^{*}\right]}_{l,r}=0.$ We will prove that $A^{*}$ is a generator of a one-parameter semigroup. Before this, let us consider a simple two reaction system to gain some intuition into the different components of the reaction counts setting.

Example 1.

Let us consider the birth-death process in the context of the reaction counts CME. The birth-death process in the species count is given by,

$$Z\left(t\right)=x_0+\mathcal{P}\left({\int }_0^t{c}_{b}ds\right)-\mathcal{P}\left({\int }_0^t{c}_{d}Z\left(s\right)ds\right),$$

where $x_0\in {\mathbb{N}}_0$ is the starting population, and the birth and death rates are denoted by $c_b$ and $c_d,$ respectively. We can translate this process into the reaction counts setting using the mapping,

$${\mathrm{\Gamma }}_{x_0}(r_1,r_2)=x_0+r_1-r_2.$$

We restrict the reaction counts state space $\mathrm{\Lambda }\subset {\mathbb{N}}_0\times {\mathbb{N}}_0$ to only contain states which yield non-negative values by applying ${\mathrm{\Gamma }}_{x_0}.$ Substituting the mapping ${\mathrm{\Gamma }}_{x_0}$ and the process $Z\left(t\right)$ into Equations (2) to (5) leads to the reaction counts birth-death process formulation:

$$R_{x_0}\left(t\right)={\left(\begin{array}{c}1\\ 0\end{array}\right)}^T\mathcal{P}\left({\int }_0^t{c}_{b}ds\right)+{\left(\begin{array}{c}0\\ 1\end{array}\right)}^T\mathcal{P}\left({\int }_0^t{c}_d{\mathrm{\Gamma }}_{x_0}\left({\left(R_{x_0}\left(s\right)\right)}_1,{\left(R_{x_0}\left(s\right)\right)}_2\right)ds\right).$$

(10)

The relationship between the reaction count state space and the species count state space is visualised in Figure 1A. The colours in the figure show how the species count state space partitions the reaction count state space. To demonstrate the forward moving structure of the reaction counts CME, we derive the generator for the first nine states of its state space:

$$\frac{du\left(t\right)}{dt}=A^{*}u\left(t\right),$$

where,

$$u\left(t\right)={\left(\frac{dp(R_{x_0}\left(t\right)=(0,0))}{dt},\frac{dp(R_{x_0}\left(t\right)=(1,0))}{dt},\frac{dp(R_{x_0}\left(t\right)=(0,1))}{dt},\dots \right)}^T,$$

$$A^{*}=\left[\begin{array}{cccccc}\frac{(0,0)}{-c_b-c_d{x}_0}& \frac{(1,0)}{0}& \frac{(0,1)}{0}& \frac{(2,0)}{0}& \frac{(1,1)}{0}& \frac{(0,2)}{0}\\ c_b& -c_b-c_d(x_0+1)& 0& 0& 0& 0& \dots \\ c_d{x}_0& 0& -c_b-c_d(x_0-1)& 0& 0& 0& \dots \\ 0& c_b& 0& -c_b-c_d(x_0+2)& 0& 0& \dots \\ 0& c_d(x_0+1)& c_b& 0& -c_b-c_d{x}_0& 0& \dots \\ 0& 0& c_d(x_0-1)& 0& 0& -c_b-c_d(x_0-2)& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right].$$

The initial condition is the identity vector with the total probability mass on state $(0,0).$ We can observe that due to the forward moving nature of the reaction counts setting, $A^{*}$ is a lower diagonal matrix. In Figure 1B,C, we present a birth-death process with initial state $x_0=0,$ birth rate $c_b=1.0,$ and death rate $c_d=0.1.$ Specifically, we show the probability distributions of the first three states of that process which correspond to state 1 in the species count setting. In Figure 1B, it is shown that the probability distribution of being in state 1 in the species count CME setting bounds the corresponding reaction counts CME distributions. In Figure 1C, we progressively add the first three reaction counts states corresponding to the state 1, showing that the reaction counts CME approximates the species counts CME from below.

Theorem 1.

(Sunkara 09 ([15])) Given a system with $N_s$ species and $N_r$ reactions at a starting state $x_0\in \mathrm{\Omega }.$ If $A_r$ is the generator for the reaction counts CME (6), then:

1.

There exists a permutation matrix P such that $P^T{A}^{*}P$ is lower triangular;

2.

The spectrum of $A^{*}$ is the set:

$$spec\left(A^{*}\right)=\left\{-\sum _{n=1}^{N_r}{\alpha }_n(x_0,r):r\in \mathrm{\Lambda }\right\}.$$

(11)

The property of the reaction counts setting, that the process only moves forward, aids in showing that the generator could be rewritten as a lower triangular matrix. We know that the spectrum of a lower triangular matrix are the diagonal elements of the matrix. In our case, the diagonal elements are simply the negative sum of the outgoing propensities. This then reveals that the spectrum of the generator $A^{*}$ is in the negative real numbers. Then by the spectral mapping theorem, solutions to the reaction counts CME exist. We can further exploit this “forward stepping” structure to write the analytical solutions for the probability at each state in $\mathrm{\Lambda }.$

Proposition 1.

(Sunkara 09 ([15])) Given a system with $N_s$ species and $N_r$ reactions at a starting state $x_0\in \mathrm{\Omega }.$ Let $C\in (0,1].$ If ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0}$ is the reaction counts Kurtz process (5) with $\mathrm{\Lambda }\ne \varnothing$, then the solution to (6) at $t>0$ is given by:

$$p(R_{x_0}\left(t\right)=\mathbf{0})=e^{-\alpha (x_0,r)t}C,$$

(12)

for $r=\mathbf{0},$ and,

$$p(R_{x_0}\left(t\right)=r)=\sum _{n=1}^{N_r}\left({\alpha }_n(x_0,r-{\mathbf{1}}_n){\int }_0^t{e}^{-\alpha (x_0,r)(t-s)}p(R_{x_0}\left(s\right)=r-{\mathbf{1}}_n)ds\right),$$

(13)

for $r\in \mathrm{\Lambda }\backslash \mathbf{0}.$

Combining Proposition 1 with the push forward measure in Equation (8), we see that at finite time t the solutions to the CME (1) exist. Unfortunately, since the reaction counts CME cannot reach a stationary state (in cases where the state space $\mathrm{\Lambda }$ is not finite), the proof for the existence of stationary solutions for a generalised CME is still an open problem. Now that the solutions of the reaction counts CME have been established, we can probe further into its properties. Firstly, we explore how the state space of the reaction count CME has a naturally hierarchical partitioning, and how this partitioning can be used to build a sequence of approximate sub-processes.

3. Partitioning the State Space

Definition 1.

For $m\in \mathbb{N},$ we define,

$${\nabla }_m:=\{r\in \mathrm{\Lambda }:{∥r∥}_1=m\},$$

to be the set of all states which are reachable by m steps from the origin $\mathbf{0}.$ The state space Λ naturally partitions into non-intersecting subsets,

$$\mathrm{\Lambda }={\nabla }_0\cup {\nabla }_1\cup {\nabla }_2\cup \dots .$$

If we consider truncating the state space progressively, we find that this is equivalent to the Finite State Projection with an N-step domain expander [9,18].

Lemma 1.

For $x_0\in \mathrm{\Omega },$ let ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0}$ be the reaction counts Kurtz process (5). Then for any $m\in \mathbb{N}$,

$$\sum _{r_1\in {\nabla }_m}\sum _{n=1}^{N_r}{\alpha }_n(x_0,r_1){\int }_0^{t}p(R_{x_0}\left(s\right)=r_1)ds=\sum _{\tilde{m}=m+1}^{\infty }\sum _{r_2\in {\nabla }_{\tilde{m}}}p\left(R_{x_0}\left(t\right)\right)=r_2).$$

(14)

Proof. 

Fix $m\in \mathbb{N}.$ Let ${\left(R_{x_0,{\nabla }_m}\left(t\right)\right)}_{t\ge 0}$ be the same process as ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0},$ but with the restriction that all propensities of states $r\in {\cup }_{\tilde{m}=m+1}^{\infty }{\nabla }_{\widehat{m}}$ are zero. Since in the reaction counts case, states only transition forward, we have that,

$$p(R_{x_0}\left(t\right)=r)=p(R_{x_0,{\nabla }_m}\left(t\right)=r),\mathrm{for}\mathrm{all}r\in {\cup }_{\tilde{m}=0}^m{\nabla }_{\tilde{m}}.$$

(15)

That is, the probability of both processes are the same for states leading up to the states in ${\nabla }_m.$ However, since the process $R_{x_0,{\nabla }_m}\left(t\right)$ does not evolve past ${\nabla }_m$, by the conservation of probability we have that,

$$\begin{array}{cc}\hfill \sum _{\tilde{m}=0}^m\sum _{r\in {\nabla }_{\tilde{m}}}p(R_{x_0,{\nabla }_m}\left(t\right)=r)+\sum _{r_1\in {\nabla }_m}\sum _{n=1}^{N_r}{\alpha }_n(x_0,r_1){\int }_0^{t}p(R_{x_0,{\nabla }_m}\left(s\right)=r_1)ds& =1.\hfill \end{array}$$

(16)

Substituting in Equation (15) for states which appear before ${\nabla }_{m+1}$ reduces the above expression to,

$$\begin{array}{cc}\hfill \sum _{\tilde{m}=0}^m\sum _{r\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r)+\sum _{r_1\in {\nabla }_m}\sum _{n=1}^{N_r}{\alpha }_n(x_0,r_1){\int }_0^{t}p(R_{x_0}\left(s\right)=r_1)ds& =1.\hfill \end{array}$$

(17)

By the conservation of probability,

$$\sum _{\tilde{m}=0}^{\infty }\sum _{r_2\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r_2)=1,$$

hence, substituting this into the right-hand side of Equation (17) and subtracting the like terms gives,

$$\sum _{r_1\in {\nabla }_m}\sum _{n=1}^{N_r}{\alpha }_n(x_0,r_1){\int }_0^{t}p(R_{x_0}\left(s\right)=r_1)ds=\sum _{\tilde{m}=m+1}^{\infty }\sum _{r_2\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r_2).$$

 ☐

Lemma 1 is an alternative formulation of the principle behind the Finite State Projection method. We extend on this result and can show that for a desired error $\epsilon \in (0,1)$, there exists a subset of the state space $\mathrm{\Lambda }$ which will produce an approximation with the desired error.

Theorem 2.

For $x_0\in \mathrm{\Omega },$ let ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0}$ be the reaction counts Kurtz process (5). For $\epsilon >0$ there exists $m\in \mathbb{N}$ such that,

$$\sum _{\tilde{m}=m+1}^{\infty }\sum _{r\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r)<\epsilon .$$

Proof. 

The proof is an extension of Lemma 1. We define,

$${\varphi }_m:=\sum _{r_1\in {\nabla }_m}\sum _{n=1}^{N_r}{\alpha }_n(x_0,r_1){\int }_0^{t}p(R_{x_0}\left(s\right)=r_1)ds,$$

and,

$${\psi }_{m+1}:=\sum _{\tilde{m}=m+1}^{\infty }\sum _{r_2\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r_2),$$

the left- and right-hand terms of Equation (14). Then for $m\in \mathbb{N},$ Lemma 1 can be reformulated to state,

$$\begin{array}{cc}\hfill {\varphi }_m& ={\psi }_{m+1};\hfill \end{array}$$

separating the ${\nabla }_{m+1}$ states from ${\psi }_{m+1}$ gives us,

$$\begin{array}{c}=\sum _{r\in {\nabla }_{m+1}}p(R_{x_0}\left(t\right)=r)+\underset{:={\psi }_{m+2}}{\underbrace{\sum _{\tilde{m}=m+2}^{\infty }\sum _{r_2\in {\nabla }_{\tilde{m}}}p(R_{x_0}\left(t\right)=r_2)}},\hfill \\ =\sum _{r\in {\nabla }_{m+1}}p(R_{x_0}\left(t\right)=r)+{\psi }_{m+2}.\hfill \end{array}$$

Combining the last and first term in the equalities gives us,

$$\begin{array}{cc}\hfill {\psi }_{m+1}& =\sum _{r\in {\nabla }_{m+1}}p(R_{x_0}\left(t\right)=r)+{\psi }_{m+2}.\hfill \end{array}$$

Rearranging the equality reduces to, ${\psi }_{m+2}-{\psi }_{m+1}<0.$ We have shown that ${\psi }_{•}$ is monotonically decreasing. By definition, ${\psi }_{•}$ has the properties that ${\psi }_1=1.0-p(R_{x_0}\left(t\right)=\mathbf{0})$ and ${lim}_{m\to \infty }{\psi }_m=0.$ Hence, there exists an m such that ${\psi }_m<\epsilon .$ ☐

Theorem 2 is an alternative proof for the fact that, when studying transient dynamics (finite time) of the CME for an arbitrary precision, one can always find a finite state space to project the CME onto. The key structure that we used to prove this is that the state space of reaction counts has a natural partitioning and that the probability only flows forward. This gave us the monotonicity needed to prove that state space truncation of the CME is well-posed. In the next section, we further decompose each state in the reaction counts state space into paths over the state space.

4. Paths

Definition 2.

For $m\in \{2,3,\dots \}$, a vector $g=(g_1,\dots ,g_m)\in {\mathrm{\Lambda }}^m$ is said to be an admissible path if for every index $i\in \{2,\dots ,m\}$ there exists an $n\in \{1,\dots ,N_r\}$ such that $g_i=g_{i-1}+{\mathbf{1}}_n,$ where ${\mathbf{1}}_n$ is the identity vector with one in the nth position.

Definition 3.

For $m\in \{2,3,\dots \},$ we denote the set of all admissible paths of length m by,

$${\gamma }_m:=\{g\in {\mathrm{\Lambda }}^m:g\mathrm{is}\mathrm{admissible}\}$$

Definition 4.

For a point $r\in \mathrm{\Lambda },$ we denote the set of all admissible paths of length ${∥r∥}_1+1$ which start at the origin $\mathbf{0}$ and end at the state r by,

$$G\left(r\right):=\{g\in {\gamma }_{{∥r∥}_1+1}:g_1=\mathbf{0}\mathrm{and}{g}_{{∥r∥}_1+1}=r\}.$$

4.1. Path Chains

We define a Markov chain over a path in the reaction counts state space. The chain must be such that it mimics a “realisation” over the reaction state space. Using Proposition 1, we define a path chain over an admissible path in the reaction state space. It is important to note that a path chain is a continuous-time Markov chain (CTMC) [19], where all reactions which transition the state off the path of the chain are accrued into a sink state. For brevity, we simply remove this component and only describe the states of interest (see Figure 2). Also, given all chains are from the same reaction counts setting, we omit explicitly stating $x_0$ in the propensity functions. We now define a path chain.

Definition 5.

For $m\in \mathbb{N},$$g\in {\gamma }_m,$ and $C\in (0,1],$ we define a path chain to be the stochastic process ${\left(X_{g,C}\left(t\right)\right)}_{t\ge 0},$ with state space $\{g\cup \{sink\left\}\right\}$ and the probability distribution given by,

$$p(X_{g,C}\left(t\right)=g_1):=e^{-\alpha \left(g_1\right)t}C,$$

(18)

$$p(X_{g,C}\left(t\right)=g_k):=\beta (g_{k-1},g_k){\int }_0^t{e}^{-\alpha \left(g_k\right)(t-s)}p(X_{g,C}\left(s\right)=g_{k-1})ds,$$

(19)

for $k\in \{2,\dots ,m\}.$ Then to conserve probability,

$$p(X_{g,C}\left(t\right)=g_k)=sink):=1.0-\sum _{g_k\in g}p(X_{g,C}\left(t\right)=g_k).$$

The propensity function α is inherited from Equation (9) and $\beta (g_k,g_{k+1}):={\alpha }_n\left(g_k\right),$ where n is the reaction index which transitions the state $g_i$ to $g_{i+1}.$

Before stating the new proposition, let us recall the decompositions performed so far. Firstly, we showed that every state in the species state space can be decomposed into multiple reaction states in the reaction count setting. That is, a population state can be decomposed into all the different ways in which it can be visited. We then ventured further and showed that a path on the species counts state space corresponds to a path with only forward stepping transitions in the reaction counts state space. Now, we extend this further by showing that any state in the reaction counts state space can be decomposed into the sum of all independent path chains which start at the origin and end in that state.

Proposition 2.

For $x_0\in \mathrm{\Omega },$ let ${\left(R_{x_0}\left(t\right)\right)}_{t\ge 0}$ be the Reaction counts Kurtz process (5).

For $r\in \mathrm{\Lambda }\backslash \mathbf{0},$ let

$$G\left(r\right):=\{g\in {\gamma }_{{∥r∥}_1+1}:g_1=\mathbf{0}and{g}_{{∥r∥}_1+1}=r\},$$

be the set of all admissible paths of length ${∥r∥}_1+1$ which start at $\mathbf{0}$ and ending at r (Definition 4).

Then,

$$p(R_{x_0}\left(t\right)=r)=\sum _{g\in G\left(r\right)}p(X_{g,1}\left(t\right)=r).$$

(20)

Proof. 

The Proposition will be proved using mathematical induction. Firstly, we show the base case; that every state which is reachable by one step from the origin satisfies Equation (20). Let $r=I_n,$ with $n\in \{1,\dots ,N_r\}.$ Applying Proposition 1 we know that,

$$p(R_{x_0}\left(t\right)=I_n)={\alpha }_n\left(\mathbf{0}\right){\int }_0^t{e}^{-a\left(I_n\right)(t-s)}p(R_{x_0}\left(s\right)=\mathbf{0})ds.$$

(21)

Now we consider the path chains leading to $I_n.$ The state can only be reached in one way, hence, we have that $G\left(I_n\right)=\{g=(\mathbf{0},I_n)\}.$ Now calculating the probability of the evolution of a path chain (Definition 5) over $g\in G\left(I_n\right)$ gives us,

$$p(X_{g,1.0}\left(t\right)=I_n):=\beta (\mathbf{0},I_n){\int }_0^t{e}^{-\alpha \left(I_n\right)(t-s)}p(X_{g,1}\left(s\right)=\mathbf{0})ds.$$

(22)

Given the two equations above match, we can conclude that for $r\in {\nabla }_1,$

$$p(R\left(t\right)=r)=\sum _{g\in G\left(r\right)}p(X_{g,1}\left(t\right)=r).$$

Now we perform the inductive step. Fix $m>1$. Assume that for $\check{r}\in {\nabla }_m,$

$$p(R_{x_0}\left(t\right)=\check{r})=\sum _{g\in G\left(\check{r}\right)}p(X_{g,1}\left(t\right)=\check{r}).$$

(23)

Then for $r\in {\nabla }_{m+1},$ Proposition 1 states that,

$$\begin{array}{cc}\hfill p\left(R\right(t)=r)=& \sum _{n=1}^{N_r}{\alpha }_n(r-I_n){\int }_0^t{e}^{-\alpha \left(r\right)(t-s)}p(R_{x_0}\left(s\right)=r-I_n)ds,\hfill \end{array}$$

using the inductive assumption Equation (23) on the right-hand side inside the integral,

$$\begin{array}{cc}\hfill =& \sum _{n=1}^{N_r}{\alpha }_n(r-I_n){\int }_0^t{e}^{-\alpha \left(r\right)(t-s)}\sum _{g\in G(r-I_n)}p(X_{g,1}\left(s\right)=r-I_n)ds.\hfill \end{array}$$

Rewriting the propensities in the path propensity notation gives,

$$\begin{array}{cc}\hfill =& \sum _{n=1}^{N_r}\beta (r-I_n,r){\int }_0^t{e}^{-\alpha \left(r\right)(t-s)}\sum _{g\in G(r-I_n)}p(X_{g,1}\left(s\right)=r-I_n)ds,\hfill \\ \hfill =& \sum _{n=1}^{N_r}\sum _{g\in G(r-I_n)}\beta (r-I_n,r){\int }_0^t{e}^{-\alpha \left(r\right)(t-s)}p(X_{g,1}\left(s\right)=r-I_n)ds.\hfill \end{array}$$

Since all paths in $G\left(r\right)$ are paths from the set ${\cup }_{n=1}^{N_r}G(r-I_n),$ with r as the last state, the two summations reduce to give,

$$\begin{array}{cc}\hfill =& \sum _{\widehat{g}\in G\left(r\right)}\beta ({\widehat{g}}_m,r){\int }_0^t{e}^{-\alpha \left(r\right)(t-s)}p(X_{\widehat{g},1}\left(s\right)={\widehat{g}}_m)ds.\hfill \end{array}$$

By Definition 5, the term above captures all the paths in $G\left(r\right)$, therefore,

$$\begin{array}{cc}\hfill =& \sum _{\widehat{g}\in G\left(r\right)}p(X_{\widehat{g},1}\left(s\right)=r)ds.\hfill \end{array}$$

We have shown that given Equation (23) holds on ${\nabla }_m,$ then Equation (23) also holds for ${\nabla }_{m+1}.$ Hence, by mathematical induction, for any $r\in \mathrm{\Lambda },$

$$p(R_{x_0}\left(t\right)=r)=\sum _{g\in G\left(r\right)}p(X_{g,1}\left(t\right)=r)ds.$$

 ☐

We have shown that each state in the reaction counts CME can be decomposed into the framework of “path chains” that we have defined. In essence, we have shown that the path chains are trajectories of the reaction counts Kurtz process. Even though at first glance the result is not a startling revelation, the real novelty lies in how we derived the probability distributions for these trajectories. By carefully decomposing the state space and concurrently deriving the corresponding probability distributions, we were able to give an analytical form for the probability of individual trajectories of the Kurtz process. Using the probability distribution given over the path chains in Definition 5, we can derive the classical stochastic simulation algorithm [20]. Using the path chains framework we have introduced so far, we can build realisations for applications such as: Hitting time problems, rare event problems, observable problems, etc. In the following section, we further build on the properties of path chains by proposing a method for estimating the path probabilities of trajectories.

4.2. Gated- and Un-Gated Path Chains

Given a path or trajectory, we define the path probability to be the probability distribution of the last state in the chain. We can use Definition 5 to analytically compute the probability of any path. However, we can exploit the properties of the path chains further and derive some simple estimates for the path probability. We borrow the notion of gating and un-gating introduced by Sunkara (2013) [21], and show how these can be used to calculate estimates for the path probabilities.

The concept of gating is mathematically tedious to prove, however, the concept is fairly simple. Given a chain g of length $m,$ gating this chain at position $j<m$ means that we set the propensity of leaving state $g_j$ to zero. Hence, the probability will simply flow into state $g_j$ and remain there (hence, the term “gating”). Then we can “un-gate” the state $g_j,$ which means that we reset the time to zero, set the initial probability at state $g_j$ to the accumulated probability from the gating, and then allow the process to continue. The reason we are interested in this notion of gating and un-gating is that for path chains, it happens so that gating and un-gating the chains gives an upper bound for the original path probability. We will now prove that this is indeed the case.

Definition 6.

For an admissible path $g\in {\gamma }_m,$ we define a truncated path $g^j:=(g_j,g_{j+1},\dots ,g_m),$ where the first $j-1$ terms of the path are ignored.

Definition 7.

Fix $\mathrm{\Delta }t<\infty$ and $x_0\in \mathrm{\Omega }.$ Let $g\in {\gamma }_m$ be a path of length $m+1.$ We define ${\left(X_{g,C}\left(t\right)\right)}_{t\in [0,\mathrm{\Delta }t]}$ to be a path chain over g in the time interval $[0,\mathrm{\Delta }t].$ Then for $j\in \{2,\dots ,m\}$ we define an un-gated path chain over $g^j$ to be given by:

$$p(X_{g^j,C}\left(t\right)=g_k):=\beta (g_{k-1},g_k){\int }_0^t{e}^{-\alpha \left(g_k\right)(t-s)}p(X_{g,C_j}\left(s\right)=g_{k-1})ds,$$

(24)

for $k\in \{j+1,\dots ,m\}$ and,

$$p(X_{g^j,C}\left(t\right)=g_j):=e^{-\alpha \left(g_j\right)t}{C}_j$$

where,

$$C_j:=\beta (g_{j-1},g_j){\int }_0^{\mathrm{\Delta }t}p(X_{g,C}\left(s\right)=g_{j-1}).$$

To conserve probability,

$$p(X_{g^j,C}\left(t\right)=g_k)=sink):=1.0-\sum _{g_k\in g^j}p(X_{g,C}\left(t\right)=g_k).$$

We see that the un-gated path chain is a normal path chain with the initial probability reset to the total probability which has flown into the gated state (see Figure 3). We now show that gating and un-gating over a path gives an upper bound for the path probability. First we consider the case of a single gating and un-gating at the position $g_2.$

Lemma 2.

Fix $g\in {\gamma }_m$ and $\mathrm{\Delta }t>0.$ Given two path chains ${\left(X_{g,C}\left(t\right)\right)}_{t\in [0,\mathrm{\Delta }t]}$ and ${\left(X_{g^2,C_2}\left(t\right)\right)}_{t\in [0,\mathrm{\Delta }t]},$ it follows that for all $t\in [0,\mathrm{\Delta }t],$

$${\int }_0^{t}p(X_{g,C}\left(s\right)=2)dt\le {\int }_0^{t}p(X_{g^2,C_2}\left(s\right)=2)dt,$$

and for $j\in \{3,\dots ,m\}$

$$p(X_{g,C}\left(t\right)=g_j)\le p(X_{g^2,C_2}\left(t\right)=g_j).$$

Proof. 

We break down the proof into two steps. First we investigate the bounds formed by gating at the second step, and then we prove the bound over the remaining states in the chain. We begin by showing that for $t\in [0,\mathrm{\Delta }t$],

$${\int }_0^{t}p(X_{g,C}\left(s\right)=2)dt\le {\int }_0^{t}p(X_{g^2,C_2}\left(s\right)=2)dt.$$

(25)

From the definition of the path chain (Definition 5), the solution for the probability of being in state $g_2$ at time t is given by,

$$p(X_{g,C}\left(t\right)=g_2)=\frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)-\alpha \left(g_2\right)}\left(e^{-\alpha \left(g_2\right)t}-e^{-\alpha \left(g_1\right)t}\right).$$

(26)

Using the solution above, we can calculate the probability which would accumulate in $g_2$ if the chain was gated at $g_2,$

$$C_2:=\frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)}\left(1-e^{-\alpha \left(g_1\right)\mathrm{\Delta }t}\right).$$

(27)

Substituting the above terms into the definition of the probability of being is state $g_2$ at time t in the un-gated path chain ${\left(X_{g^2,C_2}\left(t\right)\right)}_{t\in [0,\mathrm{\Delta }t]}$ gives us,

$$p(X_{g^2,C_2}\left(t\right)=g_2)=\frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)}\left(1-e^{-\alpha \left(g_1\right)\mathrm{\Delta }t}\right)e^{-\alpha \left(g_2\right)t}.$$

(28)

Taking the integral of Equations (26) and (28) over the interval $[0,t]\subseteq [0,\mathrm{\Delta }t]$ gives us,

$$\begin{array}{cc}\hfill {\int }_0^{t}p(X_{g,C}\left(s\right)=g_2)ds& =\frac{C\beta (g_1,g_2)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_1\right)}\left(\alpha \left(g_2\right)(e^{-\alpha \left(g_1\right)t}-1)-\alpha \left(g_1\right)(e^{-\alpha \left(g_2\right)t}-1)\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\int }_0^{t}p(X_{g^2,C_2}\left(s\right)=g_2)ds& =\frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)\alpha \left(g_2\right)}\left(1-e^{-\alpha \left(g_1\right)\mathrm{\Delta }t}\right)\left(1-e^{-\alpha \left(g_2\right)t}\right).\hfill \end{array}$$

(30)

To simplify notation, let us define,

$$\widehat{\varphi }\left(t\right):={\int }_0^{t}p(X_{g^2,C_2}\left(s\right)=g_2)ds\mathrm{and}\varphi \left(t\right):={\int }_0^{t}p(X_{g,C}\left(s\right)=g_2)ds.$$

Then the difference between $\widehat{\varphi }\left(t\right)$ and $\varphi \left(t\right)$ reduces to,

$$\begin{array}{cc}\hfill \widehat{\varphi }\left(t\right)-\varphi \left(t\right)& =\frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)\alpha \left(g_2\right)}\left(1-e^{-\alpha \left(g_1\right)\mathrm{\Delta }t}\right)\left(1-e^{-\alpha \left(g_2\right)t}\right)\hfill \\ & -\frac{C\beta (g_1,g_2)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_1\right)}\left(\alpha \left(g_2\right)(e^{-\alpha \left(g_1\right)t}-1)-\alpha \left(g_1\right)(e^{-\alpha \left(g_2\right)t}-1)\right).\hfill \end{array}$$

(31)

Since $e^{-\alpha \left(g_1\right)\mathrm{\Delta }t}\le e^{-\alpha \left(g_1\right)t}$ for all $t\in [0,\mathrm{\Delta }t],$ we attain the lower bound,

$$\begin{array}{c}\ge \frac{C\beta (g_1,g_2)}{\alpha \left(g_1\right)\alpha \left(g_2\right)}\left(1-e^{-\alpha \left(g_1\right)t}\right)\left(1-e^{-\alpha \left(g_2\right)t}\right)\hfill \\ -\frac{C\beta (g_1,g_2)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_1\right)}\left(\alpha \left(g_2\right)(e^{-\alpha \left(g_1\right)t}-1)-\alpha \left(g_1\right)(e^{-\alpha \left(g_2\right)t}-1)\right).\hfill \end{array}$$

(32)

We show that the right-hand side of the equation above is always positive. We define the right hand-side of Equation (32) as a function,

$$\begin{array}{c}\hfill \mathrm{\Psi }(g_1,g_2,t):=\frac{C\beta (g_1,g_2)\left(e^{-\alpha \left(g_2\right)t}\alpha \left(g_2\right)-e^{-\alpha \left(g_1\right)t}\alpha \left(g_1\right)+e^{-(\alpha \left(g_1\right)+\alpha \left(g_2\right))t}(\alpha \left(g_1\right)-\alpha \left(g_2\right))\right)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_2\right)}.\end{array}$$

(33)

We notice that the function $\mathrm{\Psi }$ has the property,

$$\mathrm{\Psi }(g_1,g_2,t)=\mathrm{\Psi }(g_2,g_1,t).$$

Firstly, if we consider the case $\alpha \left(g_1\right)=\alpha \left(g_2\right),$ we can see that $\psi (g_1,g_2,t)=0.$ Using the symmetry of $\mathrm{\Psi },$ we only need to consider the case $\alpha \left(g_1\right)<\alpha \left(g_2\right).$ Then it follows that,

$$\begin{array}{cc}\hfill \mathrm{\Psi }(g_1,g_2,t)& =\frac{C\beta (g_1,g_2)\left(e^{-\alpha \left(g_2\right)\mathrm{\Delta }t}\alpha \left(g_2\right)-e^{-\alpha \left(g_1\right)t}\alpha \left(g_1\right)+e^{-(\alpha \left(g_1\right)+\alpha \left(g_2\right))t}(\alpha \left(g_1\right)-\alpha \left(g_2\right))\right)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_2\right)},\hfill \end{array}$$

replacing $e^{-\alpha \left(g_1\right)t}\alpha \left(g_1\right)$ with $e^{-\alpha \left(g_1\right)t}\alpha \left(g_2\right)$ gives us,

$$\begin{array}{c}\ge \frac{C\beta (g_1,g_2)\left((e^{-\alpha \left(g_2\right)t}-e^{-\alpha \left(g_1\right)t})\alpha \left(g_2\right)+e^{-(\alpha \left(g_1\right)+\alpha \left(g_2\right))t}(\alpha \left(g_1\right)-\alpha \left(g_2\right))\right)}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)\alpha \left(g_2\right)},\hfill \\ =C\beta (g_1,g_2)\left(\underset{>0}{\underbrace{\frac{(e^{-\alpha \left(g_2\right)t}-e^{-\alpha \left(g_1\right)t})}{(\alpha \left(g_1\right)-\alpha \left(g_2\right))\alpha \left(g_1\right)}}}+\underset{>0}{\underbrace{\frac{e^{-(\alpha \left(g_1\right)+\alpha \left(g_2\right))\mathrm{\Delta }t})}{\alpha \left(g_1\right)\alpha \left(g_2\right)}}}\right),\hfill \\ \ge 0.\hfill \end{array}$$

Hence, the right-hand side of Equation (33) is positive. So it follows that for $t\in [0,\mathrm{\Delta }t],$

$${\int }_0^{t}p(X_{g,C}\left(s\right)=2)ds\le {\int }_0^{t}p(X_{g^2,C_2}\left(s\right)=2)ds.$$

We have shown that the outflow of probability in the un-gated chain is more than in the original path chain. We now show that this phenomenon continues all the way down the chain. We prove this using mathematical induction. First, let us consider the base case $j=3.$ By the definition of path chains,

$$\begin{array}{cc}\hfill p(X_{g,C}\left(t\right)=g_3)-p(X_{g^2,C_2}\left(t\right)=g_3)& =\beta (g_2,g_3){\int }_0^t{e}^{-\alpha \left(g_3\right)(t-s)}\left(p(X_{g,C}\left(s\right)=g_2)-p(X_{g^2,C_2}\left(s\right)=g_2)\right)ds,\hfill \\ & \le \beta (g_2,g_3){\int }_0^t\left(p(X_{g,C}\left(s\right)=g_2)-p(X_{g^2,C_2}\left(s\right)=g_2)\right)ds,\hfill \end{array}$$

applying condition (25) gives us that,

$$\begin{array}{c}\le 0.\hfill \end{array}$$

Hence, for $t\in [0,\mathrm{\Delta }],$

$$p(X_{g,C}\left(t\right)=g_3)\le p(X_{g^2,C_2}\left(t\right)=g_3).$$

We now build the inductive step. Fix $k\in \{4,\dots ,m\},$ assume that, for all $t\in [0,\mathrm{\Delta }t],$

$$p(X_{g,C}\left(t\right)=g_k)\le p(X_{g^2,C_2}\left(t\right)=g_k).$$

(34)

Then by the definition of path chains,

$$\begin{array}{cc}\hfill p(X_{g,C}\left(t\right)=g_{k+1})-p(X_{g^2,C_2}\left(t\right)=g_{k+1})& =\beta (g_k,g_{k+1}){\int }_0^t{e}^{-\alpha \left(g_{k+1}\right)(t-s)}\left(p(X_{g,C}\left(s\right)=g_k)-p(X_{g^2,C_2}\left(s\right)=g_k)\right)ds,\hfill \\ & \le \beta (g_k,g_{k+1}){\int }_0^t\left(p(X_{g,C}\left(s\right)=g_k)-p(X_{g^2,C_2}\left(s\right)=g_k)\right)ds,\hfill \end{array}$$

applying the inductive assumption Equation (34) gives us that,

$$\begin{array}{c}\le 0.\hfill \end{array}$$

Hence, by mathematical induction, for $j\in \{3,\dots ,m\}$ and $t\in [0,\mathrm{\Delta }t],$ the probabilities of the un-gated chain bound the probability on the original path chain, that is,

$$p(X_{g,C}\left(t\right)=g_j)\le p(X_{g^2,C_2}\left(t\right)=g_j).$$

 ☐

The intuition behind Lemma 2 is that once we gate and un-gate a state, the un-gated chain has more probability in it than the regular path chain beyond the gated state. Since the dynamics of the path chain beyond the gated state are not changed, we would expect that the un-gated chain (having more probability) would be an upper bound for the original path chain. The proof, as seen above, is sadly tedious. We now consider a cascade of gating and un-gating, that is, we gate and un-gate as we transition through the path chain. We prove that this gives an upper bound for the path probabilities.

4.3. Cascade of Gating and Un-Gating

We now set up a cascade of gating and un-gating at ever state of the path chain. We achieve this by defining the starting probability of the un-gating. Previously, an un-gated path chain at position j had a initial probability of,

$$C_j:=\beta (g_{j-1},g_j){\int }_0^{\mathrm{\Delta }t}p(X_{g,C}\left(s\right)=g_{j-1})ds.$$

We now introduce a recurrent initial probability,

$${\tilde{C}}_j:=\beta (g_{j-1},g_j){\int }_0^{\mathrm{\Delta }t}p(X_{g^{j-1},{\tilde{C}}_{j-1}}\left(s\right)=g_{j-1})ds,$$

(35)

for $j\in \{2,\dots ,m-1\}.$ Using this new recursive initial probability, we will build a sequence of path chains and using the idea of gating and un-gating at each transition of the path chain, we can construct an upper bound for the path probability.

Theorem 3.

Fix $g\in {\gamma }_m$ and $\mathrm{\Delta }t>0.$ Let ${\left(X_{g,C}\left(t\right)\right)}_{t\in [0,\mathrm{\Delta }t]}$ be a path chain over $g.$ We recursively define a set of un-gated chains, $\{{(X_{g^j,{\tilde{C}}_j}\left(t\right))}_{t\in [0,\mathrm{\Delta }t]}:j\in \{2,\dots ,m-2\}\}.$ Then for $k\in \{3,\dots ,m-1\},$

$${\int }_0^{\mathrm{\Delta }t}p(X_{g,C}\left(s\right)=g_k)ds\le {\int }_0^{\mathrm{\Delta }t}p(X_{g^{k-1},{\tilde{C}}_{k-1}}\left(s\right)=g_k)ds.$$

Proof. 

We prove this result by repeatedly applying Lemma 2. Fix $k\in \{3,\dots ,m-1\},$ then by the first clause of Lemma 2 we have that,

$$\begin{array}{cc}\hfill {\int }_0^{\mathrm{\Delta }t}p(X_{g^{k-1},{\tilde{C}}_{k-1}}\left(s\right)=g_k)ds& \ge {\int }_0^{\mathrm{\Delta }t}p(X_{g^{k-2},{\tilde{C}}_{k-2}}\left(s\right)=g_k)ds,\hfill \end{array}$$

since $g_k$ is the third state of the chain $g^{k-2},$ we can apply the second clause of Lemma 2 to get,

$$\begin{array}{c}\ge {\int }_0^{\mathrm{\Delta }t}p(X_{g^{k-3},{\tilde{C}}_{k-3}}\left(s\right)=g_k)ds,\hfill \end{array}$$

stepping backward in the chain using Lemma 2, we reduce down to,

$$\begin{array}{c}⋮\hfill \\ \ge {\int }_0^{\mathrm{\Delta }t}p(X_{g,C}\left(s\right)=g_k)ds.\hfill \end{array}$$

 ☐

In summary, by gating and un-gating along a path chain we can easily construct estimates for the probability of the trajectories of interest. We now consider an example where we build a path chain of the birth-death process introduced in Example 1, to show how to compute the gating and un-gating steps.

Example 2.

Let us consider the birth-death process introduced in Example 1, with the parameters $x_0=2,c_b=1.0,$ and $c_d=0.1.$ We will now construct the path probability of the chain:

$$g=\left((0,0),(1,0),(1,1),(1,2),(1,3)\right).$$

The generator for this path chain is given by,

$$A^{*}=\left[\begin{array}{cccccc}\frac{(0,0)}{-c_b-2c_d}& \frac{(1,0)}{0}& \frac{(1,1)}{0}& \frac{(1,2)}{0}& \frac{(1,3)}{0}& \frac{sink}{0}\\ c_b& -c_b-3c_d& 0& 0& 0& 0\\ 0& -3c_d& -c_b-2c_d& 0& 0& 0\\ 0& 0& 2c_d& -c_b-c_d& 0& 0\\ 0& 0& 0& c_d& 0& 0\\ 2c_d& c_b& c_b& c_b& 0& 0\end{array}\right].$$

We can solve for $p(X_g,1\left(t\right)=(1,3))$ by simply taking the matrix exponential at time t and applying it to the respective initial probability. Now we will construct the gating and un-gating approximation. This approximation is calculated using the following recursive equation: Let m be the number of states in the chain $g,$ then,

$$u_n\left(t\right)=\frac{\beta (g_{n-1},g_n)}{\alpha \left(g_{n-1}\right)}\left(1.0-e^{-\alpha \left(g_{n-1}\right)t}\right)u_{n-1}\left(t\right),$$

for $n=2,\dots ,m$ with $u_1=1.0.$ The functions α and β are as given in Definition 5. Hence, $u_5\left(t\right)$ is the gated and un-gated estimate for the probability $p(X_{g,1}\left(t\right)=(1,3)).$ In Figure 4A, we plot the true path probability, $p(X_{g,1}\left(t\right)=(1,3)),$ and the upper bound for the path probability, $u_5\left(t\right),$ for time interval $t\in [0,10]$ seconds. The graph shows that $u_5\left(t\right)$ is an upper bound for $p(X_g,1\left(t\right)=(1,3))$ at any time point. It is interesting to observe that the error between the two functions decreases over time, however it is not clear whether the difference will converge to zero as time goes to infinity. To test whether this behaviour can be observed in other parameter settings, we plot the analytical path probability and its upper bound for three different parameter settings in Figure 4B. We observe the same effect in all three presented cases. This raises the question whether this observation is a structural property or only an artefact of the specific model which is considered here.

5. Conclusions

We began by formulating the reaction counts CME and showing that its forward stepping characteristic yields analytical solutions. Then, with a simple application of the push forward measure, we could prove the existence of solutions of the CME for finite time. We further decomposed the reaction counts state space into independent CTMCs and gave analytical forms and estimates for their path probabilities. Hence, we have derived analytical theory for trajectories which arise from the solutions of the CME. This theory can be used for computing observable estimates of the underlying CTMCs. We can also use the paths to estimate the hitting times and rare event probabilities. A natural future direction would be to investigate if the current framework of the time independent propensities can be translated to the time dependent propensities. An even more challenging task would be, in the context of parameter inference, to reverse time on the trajectories and to cascade probability backwards.