Slack Reactants: A State-Space Truncation Framework to Estimate Quantitative Behavior of the Chemical Master Equation

Consider a mass-action 1-dimensional birth-death model,

Since the purpose of this new network is to approximate the qualitative behavior of the original system (5), we minimize the contribution of the slack reactant $Y$ for modeling the associated stochastic system. Hence, we assign $Y$ a special reaction intensity – instead of $\lambda_{Y\to X}(x,y)=\kappa_{1}y$ using mass-action, we choose

and we use the same intensity $\lambda_{X\to Y}(x,y)=\kappa_{2}x$ for the reaction $X\to Y$. By forcing $Y$ to have “zero-order kinetics,” we ensure that the computed rates remain the same throughout the chemical network except for on the imposed boundary. This choice of reaction intensities not only preserves the conservative law $X(t)+Y(t)=X(0)+Y(0)$ but also prevent the slack reactant $Y$ from having negative counts with the characteristic term $\mathbbm{1}_{\{y\geq 0\}}$.

In general, by using slack reactants, any conservative laws can be deduced in a similar fashion to encode the desired truncation directly in the CRN. We have found this perspective to be advantageous with respect to studying complex CRNs. Rather than thinking about the software implementation of the truncation, it is often easier to design the truncation in terms of slack reactants, then implement the already-truncated system exactly.

We now provide an algorithm to automatically determine the slack reactants required for a specified truncation. It is often the case that a “natural” formulation arises (typically by replacing zero-order complexes with slack reactants) but when that is not the case, one can still find a useful approximation by following our algorithm.

Consider any CRN $(\mathcal{S},\mathcal{C},\mathcal{R},\Lambda)$. Suppose we wish to apply a number of conservation bounds to reduce the state-space of the associated chemical master equation, e.g. many equations of the form

for $i=1,2,\dots,m$. Collecting the conservation vectors $\{\mathbf{w}_{i}\}_{i=1}^{m}$ into the rows of the $m\times d$ matrix $W$ and each integer $\{N_{i}\}$ into the $m\times 1$ vector $\mathbf{N}$, we may concisely write

where $\preceq$ denotes an element-wise less-than-or-equal operator and $\mathbf{X}$ is the column vector of reactants ${\{X_{i}\}}$.

We augment the $d\times 1$ vector $\mathbf{X}$ with a $m\times 1$ vector $\mathbf{Y}$ of slack reactants to convert the inequality to an equality:

where $I$ is the $m\times m$ identity matrix. This matrix multiplication implies $m$ conservative laws among $X_{i}$’s and slack reactants $Y_{i}$’s. In other words, we inject the variables $\mathbf{Y}=(Y_{1},Y_{2},\dots,Y_{m})$ into the existing reactions in such a way as to enforce these augmented conservation laws.

To do so, we need to properly define several CRN quantities. Let $S$ be the $|\mathcal{C}|\times|\mathcal{R}|$ connectivity matrix: if the $r$-th reaction establishes an edge from the $i$-th complex to the $j$-th, then $-S_{ir}=S_{jr}=1$ and the other entries in column $r$ are 0. Let $C$ be the $|\mathcal{S}|\times|\mathcal{C}|$ complex matrix $C$. Each column of $C$ corresponds to a unique complex and each row corresponds to a reactant. The entry $C_{ij}$ is the multiplier of $X_{i}$ in complex $j$ (which may be implicitly 0). The stoichiometry matrix is simply ${\Gamma=CS}$, whose $i$-th column indicates the reaction vector of the $i$-th reaction. Note that the CRN admits a conservation law $r_{1}X_{1}(t)+r_{2}X_{2}(t)+\cdots+r_{d}X_{d}(t)=r_{1}X_{1}(0)+r_{2}X_{2}(0)+\cdots+r_{d}X_{d}(0)$ if and only if $\mathbf{r}^{\mkern-2.0mu\top}\Gamma=\mathbf{0}$, where $r=(r_{1},r_{2},\dots,r_{d})^{\mkern-2.0mu\top}$.

Now consider a modified network where we inject $Y_{1},Y_{2},\dots,Y_{K}$ into the complexes. Let $D_{ij}$ denote the (unknown) multiplier of $Y_{i}$ in the $j$-th complex. We require the following to be true:

Then by using the same connectivity matrix $S$ and a new complex matrix $\widetilde{C}=\begin{pmatrix}C\\ D\end{pmatrix}$, we obtain a new CRN that enforces the conservative laws (10). To find $D$ satisfying (11), we choose $D$ such that

for an $m\times|\mathcal{C}|$ matrix $U$ with each row of the form $u_{i}\mathbf{1}^{\top}=u_{i}(1,1,\dots,1)^{\mkern-2.0mu\top}$ for some positive integer $u_{i}$. Since $\mathbf{1}^{T}S=\mathbf{0}$, this guarantees (11).

A new network can be generated with the complex matrix $\widetilde{C}$ and the connectivity matrix $S$. We call this network a regular slack network with species $\widetilde{\mathcal{S}}$ and complexes $\widetilde{\mathcal{C}}$. For the $j$-th reaction $\nu_{j}\to\nu^{\prime}_{j}\in\mathcal{R}$ corresponding to the $j$-th column of the matrix $CS$, we denote by $\tilde{\nu}_{j}\to\tilde{\nu}^{\prime}_{j}\in\widetilde{\mathcal{R}}$ the reaction corresponding to the $j$-th column of the matrix $\widetilde{C}S$. That is, $\tilde{\nu}_{j}\to\tilde{\nu}^{\prime}_{j}$ is the reaction obtained from $\nu_{j}\to\nu^{\prime}_{j}$ by adding slack reactants.

We model the slack network by $X^{\mathbf{N}}(t)=(X^{\mathbf{N}}_{1}(t),X^{\mathbf{N}}_{2}(t),\dots,X^{\mathbf{N}}_{d}(t))$ where each of the entries represents the count of species in the new network. Note that the count of each slack reactant $Y_{i}$ is fully determined by species counts $X^{\mathbf{N}}_{i}$’s because of the conservation law(s). As such, we do not explicitly model the $Y_{i}$’s.

The intensity function $\lambda^{\mathbf{N}}_{\tilde{\nu}\to\tilde{\nu}}$ of $X^{\mathbf{N}}$ for a reaction $\tilde{\nu}\to\tilde{\nu}^{\prime}$ is defined as

where $y_{i}=N_{i}-(w_{i1}x_{1}+w_{i1}x_{1}+\cdots+w_{id}x_{d})$, and $\nu\to\nu^{\prime}$ is the reaction in $\mathcal{R}$. Then we denote by $(\widetilde{\mathcal{S}},\widetilde{\mathcal{C}},\widetilde{\mathcal{R}},\Lambda^{\mathbf{N}})$ a new system with slack reactants obtained from the algorithm, where $\mathcal{K}_{N}$ is the collection of kinetic rates $\{\lambda^{N}_{\tilde{\nu}\to\tilde{\nu}^{\prime}}:\tilde{\nu}\to\tilde{\nu}^{\prime}\in\widetilde{\mathcal{R}}\}$. We refer this system with slack reactants to a slack system.

Here we highlight that the connectivity matrix $S$ and the complex matrix $C$ of the original network are preserved for a new network with slack reactants. Thus the original network and the new network obtained from the algorithm have the same connectivity. This identical connectivity allows the qualitative behavior of the original network, which solely depends on $S$ and $C$, to be inherited to the new network with slack reactants. We particularly exploit the inheritance of accessibility and the inheritance of Poissonian stationary distribution in Section 6.

The algorithm we introduced to construct a network with slack reactants provides is valid and unique up to any user-defined conservation bounds (8), and the outcome is the matrix $D$, shown in (12), that indicates the stoichiometric coefficient of slack reactants at each complex.

As we showed in Example 3.1, to minimize the ‘intrusiveness’ of a slack network, we can simplify a slack network by setting as many ${D_{ij}=0}$ as possible. To do that, we choose the entries of $\mathbf{u}$ such that $u_{i}$ is the maximum entry of the $i$-th row of $AC$. We further optimize the effect of the slack reactants by removing the ‘redundant’ stoichiometric coefficient of slack reactants. For example, for a CRN,

the algorithm generates the following new CRN with a single slack reactant $Y$

However, by breaking up the connectivity, we can also generate another network

The network in (18) is more intrusive than the network in (19) in the sense of accessibility. At any state where $Y=0$, the system associated with (18) will remain unchanged because no reaction can take place. However, the reaction $A\to Y$ in (19) can always occur despite $Y=0$. Hence (19) preserves the accessibility of the original system associated with (17) as any state for $A$ is accessible from any other state in the original reaction system (17). We refer such a system with slack reactants generated by canceling redundant slack reactants to an optimized slack system. In Section 6, we explore the accessibility of an optimized reaction network with slack reactants in a more general setting.

Finally, we can make a network with slack reactants admit a better approximation of a given CRN by choosing an optimized conservation relation in (8). First, we assume that only a single conservation law and a single slack reactant are added to a given CRN. For the purpose of state space truncation onto finitely many states, a single conservation law is enough as all species could be bounded by $N$ as shown in (8). Let this single conservation law be

Then the matrix $W$ in (9) is a vector $(w_{1},w_{2},\dots,w_{d})^{\mkern-2.0mu\top}$, and we denote this by $\mathbf{w}$. By the definition of the intensities (13) for a network with slack reactants, some reactions are turned off when $Y=0$, i.e. $w_{1}X_{1}+\dots+w_{d}X_{d}=N$. Geometrically, a reaction outgoing from the hyperplane $w_{1}X_{1}+\dots+w_{d}X_{d}=N$ is turned off (Figure 2).

Hence we optimize the estimation with slack reactants by minimizing such intrusiveness of turned off reactions. To do that, we choose $\mathbf{v}$ which minimizes the number of the reactions in $\{\nu\to\nu^{\prime}\in\mathcal{R}:(\nu^{\prime}-\nu)\cdot\mathbf{w}>0\}$.

A well known state truncation algorithm is known as the Finite State Projection (FSP) method [32]. For a given continuous-time Markov chain, the associated FSP model is restricted to a finite state space. If the process escapes this truncated space, the state is sent to a designated absorbing state (see Figure 1B). For a fixed time $t$, the probability density function of the original system can be approximated by using the associated FSP model with sufficiently many states. The long-term dynamics of the original system, however, is not well approximated because the probability flow of the FSP model leaks to the designated absorbing state in the long run.

To fix this limitation of FSP, Gupta et al. proposed the stationary Finite State Projection method (sFSP)[16]. This method also projects the original state space onto a finite state space as the FSP method intended to. But sFSP does not create a designated absorbing state, as all outgoing transitions from the projected finite state space are merged to a single state $x^{*}$ ‘inside’ the finite state space (Figure 1C). The sFSP has been frequently used to estimate the long-term distribution of discrete stochastic models. However, if the size of the truncated state space is not sufficiently large, this method could fail to provide accurate estimation for the first passage time. To demonstrate this case, we consider the following simple 2-dimensional model. In the network shown in Figure 3A, two $X_{1}$ proteins are dimerized into protein $X_{2}$ while $X_{1}$ is being produced at relatively high rate. The state space of the original model is the full 2-dimensional positive integer grid. We estimate the time until the system reaches one of the two states indicated in red in Figure 3C, and we use alternative methods to do this.

For the sFSP, we project the original state space onto the rectangle by restricting $X_{1}\leq N$ and $X_{2}\leq N$ for some $N>0$, and we fix the origin $(0,0)$ as the designated state $x^{*}$ (Figure 3C). If the process associated with the sFSP model escapes the rectangle, it transports to the designated state immediately. On the other hand, we also consider a slack network shown in Figure 3B, where we introduce the conservation law $X_{1}+2X_{2}+Y=N$ for some $N>0$. Let $\tau=\inf\{t\geq 0:X_{1}=1\text{ and }X_{2}\in\{1,2\}\}$ be the first passage time we want to estimate.

By using the inverse matrix of the truncated transition matrix for each method, as detailed in [9], we obtain the mean of $\tau$ by using different values of $N$. We also obtain an ‘almost’ true mean of $\tau$ by using $10^{5}$ Gillespie simulations of the original process. As shown in Figure 3B, the slack network model provides a more accurate mean first passage time estimation for the size of truncation in between 100 and 400 if the designated return state of sFSP is $(0,0)$.

The inaccurate estimate from the sFSP is due to the choice of a return state. The sFSP model escapes the confined space often because the production rate of $X_{1}$ is relatively high. When it returns back to the origin, it is more likely to visit the target states $\{X_{1}=1\text{ and }X_{2}\in\{1,2\}\}$ than the original process.

Figure 3 shows that the mean first passage time of this system using sFSP depends significantly on the location of the chosen designated state. One of the two states is a particularly poor choice for sFSP, but it illustrates the idea that without previous knowledge of the system it can be difficult to know which states will perform well.

We display the behavior of individual solutions of the original model, the slack network, and the sFSP model in Figure 3D. The trajectory plots show that within the time interval $[0,500]$, almost half of the 100 samples from both the original model and the slack network model stay far away from the target states, while all the 100 sample trajectories from the sFSP model stay close to the target states. We also illustrate this point in Figure 3E with heat maps of the three models at $t=500$. Note that only in the case of the sFSP, the probability densities are concentrated at the target states.

The finite buffer method was proposed to estimate the stationary probability landscape with state space truncation and a novel state space enumeration algorithm [7, 8]. For a given stochastic model associated with a CRN, the finite buffer method sets inequalities among species such as (8), so-called buffer capacities. Then at each state $\mathbf{x}$, the transition rate of a reaction $\nu\to\nu^{\prime}$ is set to be zero if at least one of the inequality does not hold at $\mathbf{x}+\nu^{\prime}-\nu$. We note the algorithm, described in Section 3.2, for generating a slack network uses the same inequalities. Thus, the finite buffer method and the slack reactant method truncate the state space in the same way. We have shown, in Section 3.3 this type of truncation can create ‘additional’ absorbing states. These additional absorbing states change the accessibility between the states which means the mean first passage times cannot be accurately estimated. However, the regular slack systems preserve the network structure of the original network. Hence we are able to prove, as we already noted, that regular slack networks inherit the accessibility of the original network as long as the original network is ‘weakly reversible’ as we will define below.

We demonstrate this disparity between the the finite buffer and slack methods with the following network. Consider the mass-action system (14) with a fixed initial state $\mathbf{x}_{0}=(a_{0},b_{0})$. We are interested in estimating the mean first passage time to a target state $\mathbf{x}_{T}=(10,10)$. Note that the state space is irreducible (i.e., every state is accessible from any other state) as the network consists of unions of cycles. This condition, the union of cycles, is precisely what is meant by weakly reversible.[13, 5] Thus, the original stochastic system has no absorbing state and is accessible to the target state $\mathbf{x}_{T}$.

To use the finite buffer method on this network, we set $2X_{A}+X_{B}\leq N$ as the buffer capacity, where $X=(X_{A},X_{B})$ is the associated stochastic process. (Here we choose $N>30$ so the state space contains the target state $\mathbf{x}_{T}$.) Hence when $X$ satisfies $2X_{A}+X_{B}=N$, the reactions $\emptyset\to A$, $B\to A$ and $\emptyset\to B$ cannot be fired as $2X_{A}+X_{B}$ exceeds the buffer capacity. We now demonstrate the system has a new absorbing state. By first using reaction $A\to B$, to deplete all $A$, and then $\emptyset\to B$, every state can reach the state $(0,N)$ in finite time with positive probability. The state $(0,N)$ is absorbing state because no other reactions can occur. Reactions $A\to\emptyset$ and $A\to B$ require at least one $A$ species, and any other reactions lead to states exceeding the buffer capacity. Therefore, the finite buffer method has introduced a new absorbing state not present in the original model so that it is not accessible to $\mathbf{x}_{T}$ with positive probability.

Now, we show the explicit network structure of our slack network formulation will preserve the accessibility of the original system. We consider the same inequality $2X_{A}+X_{B}\leq N$ as above with $N>30$. We generate the slack network by using the algorithm shown in Section 3.2,

Note that the associated stochastic process $X=(X_{A},X_{B},Y)$ admits the conservation relation $2X_{A}+X_{B}+Y=N$ implying that $2X_{A}+X_{B}\leq N$. The state $(0,N,0)$ can not be reached as the only state that is accessible to $(0,N,0)$ is $(0,N-1,2)$, but it violates the conservation law.

As we highlighted in Section 3.3, slack networks preserve the connectivity of the original network (14), hence the network (20) is also weakly reversible. Thus the state space of the stochastic process associated with the slack network is irreducible by Corollary 6.1. This implies that there is no absorbing state and the system is accessible to $\mathbf{x}_{T}$, unlike the stochastic process associated with the finite buffer relation.

Consider a Lotka-Volterra model with migration shown in Figure 4A. In this model, species $A$ is the prey, and species $B$ is the predator. Clearly, the state space this model is infinite $(A,B)$ such that $A\geq 0$, $B\geq 0$. We will use slack reactants to determine the expected time to extinction of either species. More specifically, let $K=\{(A,B):A=0\text{ or }B=0\}$. We will calculate the mean first arrival time to $K$ from an initial condition $(A(0),B(0))$. (In our simulations in Figure 4, we chose $(A(0),B(0))=(3,3)$.)

To generate our slack network, we choose a conservative bound $w\cdot(A,B)^{\top}\leq N$ with $w=(1,1)$. As we discussed in Section 3.3, this $w$ minimizes the intrusiveness of slack reactants because the number of reactions $\nu\to\nu^{\prime}$ such that $(\nu^{\prime}-\nu)\cdot w>0$ is minimized. By using the algorithm introduced in Section 3.2, we generate a regular slack network (LABEL:eq:regular_slack) with a slack reactant $Y$,

As the slack reactant $Y$ in reactions $B+Y\rightleftharpoons 2Y\rightleftharpoons A+Y$ can be canceled, we further generate the optimized slack network shown in Figure 4A. We let $A(0)+B(0)+Y(0)=N$, which is the conservation quantity of the new network.

Let $\tau$ be the first passage time from our initial condition to $K$. First, we examine the accessibility of the set $K$. Because our reaction network contains creation and destruction of all species (i.e., $B\rightleftharpoons\emptyset\rightleftharpoons A$) the original model is irreducible and any state is accessible to $K$. Furthermore, Theorem 6.1 guarantees that the stochastic model associated with the optimized slack network is also accessible to $K$ from any state.

Next, by showing there exists a Lyapunov function satisfying the condition of Theorem 5.4 for the original model, we are guaranteed the first passages times from our slack network will converge to the true first passage times. (See Appendix E.1 for more detail.) Therefore, as the plot shows in Figure 4B, the mean first extinction time of the slack network converges to that of the original model, as $N$ increases. The mean first passage time of the original model was obtained by averaging $10^{9}$ sample trajectories. These trajectories were computed serially on a commodity machine and took $4.6$ hours to run. In contrast, the mean first passage times of the slack systems were computed analytically on the same computer and took at most 13 seconds. Figure 4 also shows that using only $10^{3}$ samples is misleading as the simulation average has not yet converged to the true mean first passage time. Finally, as expected from Theorem 5.1 the probability density of the slack network converges to that of the original network (see Figure 4C).

We now consider a protein synthesis model with a toggle switch, see Figure 5A. Protein species $X$ and $Z$ may be created, but only when their respective genes $D^{X}$ or $D^{Z}$ are in the active (unoccupied) state, $D^{X}_{0}$ and $D^{Z}_{0}$. Each protein acts as a repressor for the other by binding at the promoter of the opposite gene and forcing it into the inactive (occupied) state ($D^{X}_{1}$ and $D^{Z}_{1}$). In this system we consider only one copy of each gene, so that, $D^{X}_{0}+D^{X}_{1}=D^{Z}_{0}+D^{Z}_{1}=1$ for all time. Thus, we focus primarily on the state space of protein numbers only $(X,Z)$.

The deterministic form of such systems are often referred to as a “bi-stable switch” as it is characterized by steady states $(X^{*},0)$ and ($X$ “on” and $Z$ “off”) and $(0,Z^{*})$ ($X$ “off” and $Z$ “on”). This stochastic form of toggle switch has been shown to exhibit a multi-modal long-term probability density due to switches between these two deterministic states due to rapid significant changes in the numbers of proteins $X$ and $Z$ by synthesis or degradation (depending on the state of promoters) [1]. Figure 5C shows that the associated stochastic system admits a tri-modal long-term probability density. Thus the system transitions from one mode to other modes and, for the kinetic parameters chosen in Figure 5, rarely leaves the region $R=\{(X,Z)|0\leq X\leq 30\text{ or }0\leq Z\leq 30\}$. Significant departures from a stable region of a genetic switch may be associated with irregular and diseased cell fates. As such, the first passage time of this system outside of $R$ may indicate the appearance of an unexpected phenotype in a population. Because this event is rare, estimating first passage times with direct stochastic simulations, such as with the Gillespie algorithm [15], will be complicated by the long time taken to exit the region.

As in the previous example, slack systems provide a valuable tool for direct calculation of mean first passage times. In this example, we consider the time a trajectory beginning at state $(X,Z,D^{X}_{0},D^{Z}_{0})=(0,0,1,1)$ enters the target set $K=\{(X,Z)|X>30\text{ and }Z>30\}=R^{c}$ and compute $\tau$, the first passage time to $K$.

Since the species corresponding to the status of promoters ($D^{X}_{0},D^{X}_{1},D^{Z}_{0}$ and $D^{Z}_{1}$) are already bounded, we use the conservative bound $X+Z\leq N$ to generate a regular slack network in Figure 5B with the algorithm introduced in Section 3.2. The original toggle switch model is irreducible (because of the degradation $X\to 0$, $Z\to 0$ and protein synthesis $Z+D^{X}_{0}\leftarrow D^{X}_{1}$, $X+D^{Z}_{0}\leftarrow D^{Z}_{1}$ reactions). Moreover by Theorem 6.1, the degradation reactions guarantee that the slack system is also accessible to $K$ from any state.