  Phenotypical variability in the absence of genetic variation often reflects
complex energetic landscapes associated with underlying gene regulatory
networks (GRNs). In this view, different phenotypes are associated with
alternative states of complex nonlinear systems: stable attractors in
deterministic models or modes of stationary distributions in stochastic
descriptions. We provide theoretical and practical characterizations of these
landscapes, specifically focusing on stochastic slow promoter kinetics, a time
scale relevant when transcription factor binding and unbinding are affected by
epigenetic processes like DNA methylation and chromatin remodeling. In this
case, largely unexplored except for numerical simulations, adiabatic
approximations of promoter kinetics are not appropriate. In contrast to the
existing literature, we provide rigorous analytic characterizations of multiple
modes. A general formal approach gives insight into the influence of parameters
and the prediction of how changes in GRN wiring, for example through mutations
or artificial interventions, impact the possible number, location, and
likelihood of alternative states. We adapt tools from the mathematical field of
singular perturbation theory to represent stationary distributions of Chemical
Master Equations for GRNs as mixtures of Poisson distributions and obtain
explicit formulas for the locations and probabilities of metastable states as a
function of the parameters describing the system. As illustrations, the theory
is used to tease out the role of cooperative binding in stochastic models in
comparison to deterministic models, and applications are given to various model
systems, such as toggle switches in isolation or in communicating populations,
a synthetic oscillator and a trans-differentiation network.
