Acceleration of Moment Bound Optimization for Stochastic Chemical Reactions Using Reaction-wise Sparsity of Moment Equations

1. Problem Statement

Stochastic chemical reaction systems, common in synthetic biology and cellular processes, are governed by the Chemical Master Equation (CME). Analytical solutions to the CME are generally intractable due to its infinite dimensionality. Consequently, researchers rely on moment-based approaches to characterize system dynamics (e.g., mean and variance).

However, the moment equations form an infinite hierarchy where lower-order moments depend on higher-order ones. To make this tractable, Semidefinite Programming (SDP) has been used to compute mathematically guaranteed upper and lower bounds on stationary moments without resorting to unverified approximations (like moment closure).

The Core Challenge:
The size of the SDP constraints (specifically the moment and localizing matrices) grows combinatorially with the number of molecular species ($n$) and the truncation order ($\mu$). This rapid growth makes the computational cost prohibitive for systems with more than a few species, limiting the practical application of rigorous moment bounding.

2. Methodology

The authors propose a sparsity-exploiting matrix decomposition method to reduce the dimensionality of the SDP constraints. The methodology relies on the observation that in chemical kinetics, each reaction typically involves only a small subset of molecular species (reactants), leading to a "reaction-wise" sparsity structure in the moment equations.

Key Steps:

1. Reaction-wise Decomposition of Moment Equations:
   The standard moment equation (a linear ODE for moments) is rewritten as a sum of terms, where each term corresponds to a specific chemical reaction $R_i$.
   $$0 = \sum_{i=1}^r \left( \langle A_i(\alpha), M \rangle + \sum_{k \in K_i} \langle A_i^{(k)}(\alpha), N^{(k)} \rangle \right)$$
   Here, $M$ is the global moment matrix, $N^{(k)}$ are localizing matrices, and $K_i$ represents the set of reactant species for reaction $i$.

2. Identification of Structural Sparsity:
   The authors prove (Theorem 2) that for any specific reaction $i$, only a specific subset of monomials (and thus moments) appears in its corresponding term. Many entries in the coefficient matrices $A_i$ and $A_i^{(k)}$ are structurally zero because the reaction does not involve certain species.

3. Matrix Permutation and Partitioning:
   Based on the identified non-zero monomials, the global moment matrix $M$ is permuted and partitioned into:

   • A core set of monomials ($\Phi_i$) relevant to the reaction.
   • A "unused" set ($\bar{\Phi}_i$) that does not appear in the reaction term.
     The matrix is then decomposed into smaller sub-matrices ($\hat{M}_{i,j}$) by partitioning the unused sets further.

4. Relaxed Optimization Formulation:
   Instead of enforcing positive semidefiniteness on the single, massive global matrix $M \succeq 0$, the authors reformulate the problem to enforce positive semidefiniteness on the smaller, decomposed sub-matrices ($\hat{M}_{i,j} \succeq 0$).

   • Theoretical Basis: This approach is a relaxation of the original problem. However, the authors link it to chordal decomposition, a known technique in sparse SDP. While the decomposition introduces a relaxation (potentially widening the bounds), it is structurally similar to chordal decomposition, which typically yields tight bounds.

3. Key Contributions

• Characterization of Reaction-wise Sparsity: The paper formally defines the sparsity structure of moment equations induced by the reactant sets of individual chemical reactions.
• Decomposition Algorithm: A constructive method is provided to decompose the large global moment and localizing matrices into smaller, reaction-specific matrices based on this sparsity.
• Computational Acceleration: The proposed reformulation significantly reduces the dimension of the largest semidefinite constraint, which is the dominant factor in SDP computational cost (scaling cubically with matrix dimension).
• Theoretical Justification: The method is shown to be a relaxation related to chordal decomposition, ensuring that while bounds may be slightly looser, they remain practically useful and mathematically valid.

4. Results

The authors validated their method using a negative feedback gene regulatory system with a DNA sponge (7 molecular species, 14 reactions).

• Accuracy: For truncation orders $\mu = 5$ and $\mu = 6$, the bounds obtained by the proposed method were comparable to the original (non-decomposed) method. The gap between upper and lower bounds remained small (order of $10^{-2}$).
• Performance:
  • Matrix Size: The maximum dimension of the semidefinite matrices was significantly reduced. For example, at $\mu=6$, the max dimension dropped from 120 (original) to 48 (proposed).
  • Computation Time: The proposed method reduced computation time by approximately 20% for $\mu=5$ and $\mu=6$.
  • Scalability: The reduction in matrix size suggests that the method enables the analysis of larger systems that were previously computationally intractable using standard SDP approaches.

5. Significance

This work addresses a critical bottleneck in the rigorous analysis of stochastic biological systems. By leveraging the inherent sparsity of chemical reaction networks, the authors enable the application of guaranteed moment bounding to larger, more complex systems.

• Practical Impact: It allows synthetic biologists and systems engineers to obtain rigorous performance guarantees (bounds on protein copy numbers, noise levels, etc.) for complex gene circuits without relying on potentially inaccurate moment closure approximations.
• Theoretical Contribution: It bridges the gap between chemical kinetics and sparse semidefinite programming, demonstrating that reaction-wise sparsity can be effectively exploited similarly to graph sparsity in other optimization domains.

In summary, the paper provides a scalable framework for bounding stochastic chemical kinetics, making rigorous analysis feasible for systems with higher molecular complexity.