Markov State Models of gene regulatory networks

Background Gene regulatory networks with dynamics characterized by multiple stable states underlie cell fate-decisions. Quantitative models that can link molecular-level knowledge of gene regulation to a global understanding of network dynamics have the potential to guide cell-reprogramming strategies. Networks are often modeled by the stochastic Chemical Master Equation, but methods for systematic identification of key properties of the global dynamics are currently lacking. Results The method identifies the number, phenotypes, and lifetimes of long-lived states for a set of common gene regulatory network models. Application of transition path theory to the constructed Markov State Model decomposes global dynamics into a set of dominant transition paths and associated relative probabilities for stochastic state-switching. Conclusions In this proof-of-concept study, we found that the Markov State Model provides a general framework for analyzing and visualizing stochastic multistability and state-transitions in gene networks. Our results suggest that this framework—adopted from the field of atomistic Molecular Dynamics—can be a useful tool for quantitative Systems Biology at the network scale. Electronic supplementary material The online version of this article (doi:10.1186/s12918-017-0394-4) contains supplementary material, which is available to authorized users.


Connection between the Master Equation and the Transition Matrix
Herein, we summarize the main previously derived, theoretical results and refer the reader to cited references for further detail. The discrete, Markovian Chemical Master Equation (CME) describes the time-evolution of the probability distribution over the state-space for a well-mixed system of S reacting species with M possible reactions. The state of the system is given by the molecular population vector x ∈ N S . Let p(x, t) be the probability that the system is found in the state x at time t, and let a µ (x)dt and ν µ be the propensity function and stoichiometric transition vector, respectively, of the possible reactions µ = (1, 2, ..., M ). Then the CME is given by In vector-matrix form, where the off-diagonal element K ij gives the time-independent rate of transitioning from state x i to state x j , and K ii = − j =i K ji (i.e., columns sum to 0). Given the probability at time t, p(X, t), the probability density at a later time (so-called lag-time) τ may be computed by: where exp(Kτ ) is the matrix-exponential of Kτ . Many biochemical reaction networks are modeled as open systems, such that the molecular state space is technically infinite. That is, a given molecule in the system may occur with with any positive integer value for its copy number, though practically only a finite number of states are reachable in finite time. We assume here that the CME is approximated over a finite subspace, e.g., via the Finite State Projection Algorithm [1].
The linear system may be recast using the row-stochastic transition matrix [2] T(τ ) = exp(K T τ ) where the off-diagonal elements T ij (τ ) give the probability that, given the system is in state i (corresponding to molecular population vector x i ), it will then be found in state j after lag time τ . Probability is preserved, and thus rows sum to 1. The relationship between the state reaction matrix K and the stochastic transition matrix T leads to the relationship for the implied timescales: .

Building Markov State Models with PCCA+
The goal of the clustering approach is to derive a new coarse-grained stochastic transition matrix with drastically reduced dimensionality relative to to the original matrix T, and which preserves as accurately as possible the slow system dynamics. Consider a reactive system in which T is defined over a total number of reachable states N (i.e., if a system of S molecular species has a maximum copy number per species of n, then N ∼ n S ). Then the Perron eigenvalue λ 1 = 1 corresponds to the left-eigenvector π according to that is, π T = (π 1 , ..., π N ) gives the stationary probability over N states. If the reactive system has the property of metastability, then the dynamics can be approximately decomposed into fast and slow processes, with fast transitions occuring within metastable states, and slow transitions carrying the system between metastable states. For sufficient separation of fast and slow timescales, Markovian hopping between metastable states is a good model of global dynamics, and the full dynamics may be projected onto this slow subspace. In this case, for a sufficient ordering of the N states, T has a dominant block-diagonal structure with some number C of weakly coupled blocks.
In such systems, a cluster of C eigenvalues can be found near the Perron root, with λ 1 ≥ λ i ≥ λ C corresponding to the slow processes in the system, and all fast processes corresponding to rapidly decaying processes with eigenvalues λ i < λ C . The PCCA+ algorithm [3] determines membership vectors χ which assign states i ∈ {1, ..., N } to clusters of states j ∈ {1, ..., C} (with C < N ) where grades of membership are given by χ j (i) ∈ [0, 1]. (The original PCCA algorithm was based on strict "crisp" assignments χ j (i) ∈ {0, 1}, but the algorithm was not robust numerically). The membership vectors satisfy the linear transformation where ψ = [ψ 1 , ..., ψ C ] is the N × C matrix constructed from the C dominant right-eigenvectors, and B is a non-singular matrix whose elements are determined by an unconstrained optimization procedure. The membership vectors determine a projection of the full N × N transition matrix T onto a reduced (coarse) subspace, and we denote this reduced C × C transition matrix asT. As discussed previously, the optimization procedure can be be based on maximizing the metastability, taken to be measured by trace(T). The coarse-grained stationary probability is given by the projectioñ and the transition probabilities in the reduced space are given bỹ (for the probability to remain in state i within time τ ) and (for the probabilities to jump from state i to state j in τ , where u, v π is the π-weighted inner product of u and v). Therefore we may writẽ 6 whereD = diag(π 1 , ...,π C ), or equivalently, where D = diag(π 1 , ..., π N ). One may compare, on the coarse sub-space, the true system dynamics given by T with the dynamics given by the reduced Markov State Model,T. Given an initial probability vector p(x, t 0 ), the true dynamics would propagate the system according to p T (x, kτ +t 0 ) = p T (x, t 0 )T(τ ) k . By projection, The error after k timesteps is thus given by: As discussed previously in [2] the maximum error associated with the projection of T ontoT in the reduced space (independent of initial state) is given by:

Transition Path Theory
Transition path theory [4] is used in this study to extract information about the kinetics and mechanism of a transition between a starting set of states A and ending set of states B. First, we calculate the flux between states during the transition between A and B according to (17). We obtain a net flux from (18) that extracts only the flux that contributes to the transition. Finally in (19), flux of each pathway can be used to compute the probability of traversing the pathway. Equation (20) and (21) show how to reduce the full transition matrix to a reduced representation if it is desirable to perform the pathway probability calculation on the metastable states of the system.
To compute the probabilities of the transition paths, from a starting state A to an ending state B, one must calculate the effective flux f ij for all pairs of states. The transition probability relevant to the A to B state transition is found by taking the product of forward commitor q − i , transition probability T ij , and backward commitor q + j . The forwards and backwards commitors are probablities of transitioning to B or returning to A respectively. Finally, the flux is obtained by weighting it by its starting state stationary probability.
To get the net flow from A to B, only positive fluxes are considered by taking differences of the fluxes between pairs of states.
After extracting out the pathways by removing their fluxes sequentially from the system, the probability of each pathway is calculated by taking the ratio of the flux of the individual path to the total flux of all paths.
The total stationary probabilityπ i of a cluster of states is simply the sum of its individual states π i in the cluster.π The transition probablities between clustersT IJ is found by summing transition probabilities between states T ij weighted by their stationary probabilities of the states π i and normalized by the stationary probabilities of the clusters π I .  is shown, and each microstate is colored according to its macrostate assignment, as determined by the PCCA+ crisp partitioning. Each of the sixteen panels corresponds to a particular DNA promoter binding configuration A 00 , etc., as defined in SI I, MISA reactions, and Fig. 2B. For this parameter set, the promoter configuration determined the macrostate assignment exactly.