Quantifying epigenetic stability with minimum action paths

Abstract

Chromatin can adopt multiple stable, heritable states with distinct histone modifications and varying levels of gene expression. Insight on the stability and maintenance of such epigenetic states can be gained by mathematical modeling of stochastic reaction networks for histone modifications. Analytical results for the kinetic networks are particularly valuable. Compared to computationally demanding numerical simulations, they often are more convenient at evaluating the robustness of conclusions with respect to model parameters. In this communication, we developed a second-quantization-based approach that can be used to analyze discrete stochastic models with a fixed, finite number of particles using a representation of the SU(2) algebra. We applied the approach to a kinetic model of chromatin states that captures the feedback between nucleosomes and the enzymes conferring histone modifications. Using a path-integral expression for the transition probability, we computed the epigenetic landscape that helps to identify the emergence of bistability and the most probable path connecting the two steady states. We anticipate the generalizability of the approach will make it useful for studying more complicated models that couple epigenetic modifications with transcription factors and chromatin structure.

FIG. 1.

Illustration of the kinetic model for the interconversion between modified (green, X) and unmodified nucleosomes (grey, Y). (a) Recruited conversion defined in Eq. 1 that requires a pair of (un-)modified nucleosomes to alter the state of a nucleosome. (b) Noisy conversion (Eq. 2) with first order kinetics.

FIG. 2.

Phase portrait determined using Eq. (20) with kinetic parameters c₁/c₂ = 3 (a) and 12 (b). The red dashed lines are zero-energy paths and green dots are steady state solutions. The blue paths represent deterministic trajectories. The number of nucleosomes was held fixed at N = 60.

FIG. 3.

Comparison between the steady state distribution (−logP_eq, red dots) and the quasi potential (Φ) computed using Eq. 21 (black solid line) and the Fokker-Planck equation (blue dashed line) for c₁/c₂ = 3 (a) and 12 (b). The number of nucleosomes was held fixed at N = 60.

FIG. 4.

Correlation between the exact transition rates (k) computed from diagonalizing the transition matrix and the barrier height of the quasi-potential (a) or the mean first passage time (τ) estimated using the Fokker-Planck equation (b). Each data point corresponds to an independent calculation for integer values of the parameter c₁/c₂ between 5 and 120. The total nucleosome number was fixed as N = 60.