Abstract
BACKGROUND: The dynamics of many gene regulatory networks (GRNs) is characterized by the occurrence of metastable phenotypes and stochastic phenotype switches. The chemical master equation (CME) is the most accurate description to model such stochastic dynamics, whereby the long-time dynamics of the system is encoded in the spectral properties of the CME operator. Markov State Models (MSMs) provide a general framework for analyzing and visualizing stochastic multistability and state transitions based on these spectral properties. Until now, however, this approach is either limited to low-dimensional systems or requires the use of high-performance computing facilities, thus limiting its usability. RESULTS: We present a domain decomposition approach (DDA) that approximates the CME by a stochastic rate matrix on a discretized state space and projects the multistable dynamics to a lower dimensional MSM. To approximate the CME, we decompose the state space via a Voronoi tessellation and estimate transition probabilities by using adaptive sampling strategies. We apply the robust Perron cluster analysis (PCCA+) to construct the final MSM. Measures for uncertainty quantification are incorporated. As a proof of concept, we run the algorithm on a single PC and apply it to two GRN models, one for the genetic toggle switch and one describing macrophage polarization. By comparing the results with reference solutions, we demonstrate that our approach correctly identifies the number and location of metastable phenotypes with adequate accuracy and uncertainty bounds. We show that accuracy mainly depends on the total number of Voronoi cells, whereas uncertainty is determined by the number of sampling points. CONCLUSIONS: A DDA enables the efficient computation of MSMs with quantified uncertainty. Since the algorithm is trivially parallelizable, it can be applied to larger systems, which will inevitably lead to new insights into cell-regulatory dynamics.