Reachability in Biochemical Dynamical Systems by Quantitative Discrete Approximation

Technical Summary: Reachability in Biochemical Dynamical Systems by Quantitative Discrete Approximation

Problem Statement
Systems biology aims to reconstruct executable models of living organisms to investigate emergent behaviors arising from underlying biochemistry. While Ordinary Differential Equations (ODEs) are the standard for modeling these biochemical dynamical systems, they present significant challenges for global analysis. Specifically, the reachability problem—determining whether a dynamical phenomenon (e.g., oscillation or specific concentration levels) occurs for all considered initial conditions—is generally undecidable due to the unbounded and uncountable nature of continuous state spaces.

Existing analysis methods often rely on numerical simulation, which is limited to "local" analysis starting from a single initial point. To achieve global analysis, researchers employ finite discrete approximations, such as rectangular abstractions. However, traditional rectangular abstraction techniques suffer from high levels of "spurious behavior." Because these methods over-approximate the vector field on the boundaries between rectangles (facets) without tracking the specific history of trajectories, they generate many state sequences that do not correspond to any actual trajectory in the continuous system. Increasing the partition density does not always eliminate this issue effectively.

Methodology
The authors propose a novel computational technique for finite discrete approximation of continuous dynamical systems, specifically targeting biochemical systems modeled by multi-affine vector fields (affine in each variable). The core of the method is the Quantitative Discrete Approximation Automaton (QDAA).

1. State Definition: Unlike standard rectangular abstractions where a state is merely a rectangle, a state in the QDAA is defined as a pair $\langle H, E \rangle$, where $H$ is a rectangle (a cell in the partitioned state space) and $E$ is an "entry set." The entry set $E$ represents a subset of a specific facet of $H$ through which trajectories enter the rectangle. This local memory encodes the history of the system's evolution.
2. Quantification via Focal Subsets: For each state $\langle H, E \rangle$, the method identifies "focal subsets." These are subsets of the entry set from which all trajectories evolve toward a specific exit facet (or remain within the rectangle indefinitely).
3. Weighted Transitions: Transitions between states are assigned weights based on the relative volume (Lebesgue measure) of these focal subsets. Specifically, the weight of a transition from state $A$ to $B$ corresponds to the volume of the part of $A$'s entry set that leads to $B$, divided by the total volume of $A$'s entry set. This transforms the QDAA into a discrete-time Markov chain.
4. Discretization and Approximation: To make the computation feasible, continuous volumes are approximated by discretizing facets into a uniform grid of $\kappa$-tiles. Numerical simulations are employed to identify entry regions and focal subsets. The parameter $\kappa$ controls the granularity of the discretization.
5. Convergence: The method is designed such that as the discretization parameter $\kappa$ increases, the approximate volume measure converges to the continuous volume measure, thereby reducing the discrepancy between the original continuous behavior and its approximation.

Key Contributions
The paper presents four primary contributions:

1. A Novel Computational Technique: The introduction of the QDAA for the finite discrete approximation of multi-affine dynamical systems. This technique augments standard rectangular abstraction with local memory (entry sets) and quantitative weights.
2. Convergence Proof: A theoretical demonstration (Theorem 3.3) that the QDAA converges to the original continuous system's behavior as the discretization parameter $\kappa$ approaches infinity.
3. Reachability Algorithm: The development of specific algorithms (Algorithms 1 and 2) to compute the reachable state space of the QDAA. These algorithms utilize breadth-first search and numerical simulations (including backward simulations to refine entry sets) to determine successors and transition weights.
4. Empirical Evaluation: The application and evaluation of the method on several biological models, including an oscillatory Lotka-Volterra variant, an enzyme kinetics model, and a real-world case study of ammonium assimilation in E. coli.

Results
The authors evaluated the method on models of dimensions two, four, and seven.

• Oscillatory and Enzyme Models: The results demonstrated that increasing $\kappa$ improves the precision of the reachability analysis. For the oscillatory model, the QDAA provided a more accurate representation of reachable states compared to standard rectangular abstraction, which over-approximated the reachable space to the entire phase space.
• E. coli Case Study: In the ammonium transport model, the QDAA results for the concentration of intracellular ammonium ($NH_4in$) were significantly closer to numerical simulation data than those obtained via standard rectangular abstraction. The QDAA results were within one order of magnitude of the simulation, whereas the rectangular abstraction over-approximated the bounds significantly.
• Spurious Behavior: The evaluation confirmed that the QDAA effectively reduces the extent of spurious behavior (false positives in reachability) inherent in conservative rectangular abstractions.

Significance and Claims
The paper claims that the proposed method offers a practical shift in applying formal methods to complex systems. By combining numerical simulation with rectangular abstraction, the QDAA avoids the extremes of being purely an over-approximation or under-approximation. Instead, it provides a parameterized approximation where the level of precision can be tuned.

The authors assert that the method is sufficient for approximating reachability in complex biochemical dynamical systems. It serves two potential roles:

1. As a parameterized simulation technique to obtain global views of dynamics.
2. As a tool to quantify and eliminate spurious counterexamples when used in conjunction with rectangular abstraction for model checking.

The work is presented as a full version of a conference paper, offering a theoretical foundation and a prototype implementation (in C++) that demonstrates satisfactory results for computing reachability in low-to-mid dimensional biological models. The authors note that future work includes integrating the method into existing model checking software and developing parallel algorithms to handle higher-dimensional systems.