Noise in gene expression can lead to reversible phenotypic switching. Several experimental studies have shown that the abundance distributions of proteins in a population of isogenic cells may display multiple distinct maxima. Each of these maxima may be associated with a subpopulation of a particular phenotype, the quantification of which is important for understanding cellular decision-making. Here, we devise a methodology which allows us to quantify multimodal gene expression distributions and single-cell power spectra in gene regulatory networks. Extending the commonly used linear noise approximation, we rigorously show that, in the limit of slow promoter dynamics, these distributions can be systematically approximated as a mixture of Gaussian components in a wide class of networks. The resulting closed-form approximation provides a practical tool for studying complex nonlinear gene regulatory networks that have thus far been amenable only to stochastic simulation. We demonstrate the applicability of our approach in a number of genetic networks, uncovering previously unidentified dynamical characteristics associated with phenotypic switching. Specifically, we elucidate how the interplay of transcriptional and translational regulation can be exploited to control the multimodality of gene expression distributions in two-promoter networks. We demonstrate how phenotypic switching leads to birhythmical expression in a genetic oscillator, and to hysteresis in phenotypic induction, thus highlighting the ability of regulatory networks to retain memory. (5), and cancer cells (6). It has been argued that such stochastic transitions in gene activity can affect stem cell lineage decisions (7,8). Similarly, they may present advantageous strategies when cells make decisions in changing environments (9). Here, we develop a quantitative methodology which allows us to explore the implications of phenotypic switching, and the phenomena associated with it.It is known that gene regulatory networks involving slow promoter switching may lead to distinct expression levels having significant lifetimes; hence, overall expression levels are characterized by bimodal distributions (10-12) or, more generally, by mixture distributions. However, it remains to be resolved how modeling can generally describe and parameterize these distributions. A positive resolution is crucial for the development of testable quantitative and predictive models, e.g., when investigating the sensitivity of bimodality against variation of model parameters, for estimating rate constants from experimentally measured distributions, in the design of synthetic circuitry with tunable gene expression profiles, but, most importantly, when determining the implications of phenotypic decision-making.A class of theoretical models based on the Chemical Master Equation (CME) predicts bimodal protein distributions in the absence of bistability in the corresponding deterministic model (12, 13), some of which have been verified experimentally (1,4,14). Recent efforts to quant...