Accelerated Bayesian inference of gene expression models from snapshots of single-cell transcripts

TL;DR

This paper introduces a fast Bayesian inference method for gene expression models from single-cell transcript data, enabling efficient model selection and parameter estimation with a 1000-fold speedup over existing techniques.

Contribution

The authors developed a hybrid numerical algorithm combining analytical results and kinetic Monte Carlo to accelerate Bayesian inference for gene expression models.

Findings

Accelerated inference enables practical Bayesian model comparison.

Bayesian evidence outperforms BIC and AIC in model selection.

Parameters can be reasonably constrained with modest data samples.

Abstract

Understanding how stochastic gene expression is regulated in biological systems using snapshots of single-cell transcripts requires state-of-the-art methods of computational analysis and statistical inference. A Bayesian approach to statistical inference is the most complete method for model selection and uncertainty quantification of kinetic parameters from single-cell data. This approach is impractical because current numerical algorithms are too slow to handle typical models of gene expression. To solve this problem, we first show that time-dependent mRNA distributions of discrete-state models of gene expression are dynamic Poisson mixtures, whose mixing kernels are characterized by a piece-wise deterministic Markov process. We combined this analytical result with a kinetic Monte Carlo algorithm to create a hybrid numerical method that accelerates the calculation of time-dependent mRNA distributions by 1000-fold compared to current methods. We then integrated the hybrid algorithm into an existing Monte Carlo sampler to estimate the Bayesian posterior distribution of many different, competing models in a reasonable amount of time. We validated our method of accelerated Bayesian inference on several synthetic data sets. Our results show that kinetic parameters can be reasonably constrained for modestly sampled data sets, if the model is known \textit{a priori}. If the model is unknown,the Bayesian evidence can be used to rigorously quantify the likelihood of a model relative to other models from the data. We demonstrate that Bayesian evidence selects the true model and outperforms approximate metrics, e.g., Bayesian Information Criterion (BIC) or Akaike Information Criterion (AIC), often used for model selection.

Equations

$$\mathrm{Gene}\mathrm{state}i$$

$$\mathrm{mRNA}$$

$$\mathcal{L}=$$

$$\frac{\text{d}}{\text{d}t}\lambda (t)={\beta }_{s(t)}-\delta \lambda (t)$$

$$\hat{\mathbb{P}}\left(N_{\text{mRNA}}\left(t\right)=m,s\left(t\right)=i\right)=\sum _{k=1}^{N_s}\frac{{\delta }_{i,s_k\left(t\right)}}{N_s}\frac{{\lambda }_k^m\left(t\right)e^{-{\lambda }_k\left(t\right)}}{m!},$$

$$\mathbb{P}(\vec{\theta }|h,\mathcal{M})=\frac{\mathbb{P}(h|\vec{\theta },\mathcal{M})\mathbb{P}(\vec{\theta }|\mathcal{M})}{\mathbb{P}(h|\mathcal{M})}=\frac{\mathcal{L}\cdot \mathbb{P}(\vec{\theta }|\mathcal{M})}{\mathbb{P}(h|\mathcal{M})}$$

$$\mathbb{P}(h|\mathcal{M})=\int \underset{\mathcal{L}}{\underbrace{\mathbb{P}(h|\vec{\theta },\mathcal{M})}}\underset{\text{Prior}}{\underbrace{\mathbb{P}(\vec{\theta }|\mathcal{M})}}\text{d}\vec{\theta }$$

$$\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t\right)=n\right)=$$

$$e^{-n\delta t}\Theta \left(m-n+\frac{1}{2}\right).$$

$$\mathbb{P}\left(N_{\text{mRNA}}^{\text{new}}\left(t\right)=m\right)=\frac{{\lambda }^m\left(t\right)e^{-\lambda \left(t\right)}}{m!},$$

$$\frac{\text{d}}{\text{d}t}{\mathbb{P}}_m\left(t\right)=$$

$$+{\beta }_{s\left(t\right)}{\mathbb{P}}_{m-1}\left(t\right)+\delta \left(m+1\right){\mathbb{P}}_{m+1}\left(t\right),$$

$$\mathcal{G}\left(z,t\right):=\sum _{m=0}^{\infty }{z}^m{\mathbb{P}}_m\left(t\right),$$

$${\partial }_t\mathcal{G}\left(z,t\right)=\delta \left(1-z\right){\partial }_z\mathcal{G}\left(z,t\right)+{\beta }_s\left(z-1\right)\mathcal{G}\left(z,t\right).$$

$$\mathcal{G}\left(z,t\right)=$$

$$\exp \left(\frac{{\beta }_s}{\delta }\left(z-1\right)\left(1-e^{-\delta t}\right)\right),$$

$${\mathcal{G}}^{\text{init}}\left(z,t\right)=$$

$${\mathcal{G}}^{\text{new}}\left(z,t\right)=$$

$$\mathcal{G}\left(z,t_1\right)={\mathcal{G}}^{\text{init}}\left(z,t_1\right)\times {\mathcal{G}}^{\text{new}}\left(z,t_1\right).$$

$$\mathcal{G}\left(z,t|t_1\le t\le t_2\right)={\mathcal{G}}^{\text{init}}\left(z,t\right)\cdot \exp \{\frac{z-1}{\delta }[{\beta }_{s_0}\times$$

$$\left(1-e^{-\delta t_1}\right)e^{-\delta \left(t-t_1\right)}+{\beta }_{s_1}\left(1-e^{-\delta \left(t-t_1\right)}\right)]\}.$$

$$\mathcal{G}$$

$$\exp \left[\frac{z-1}{\delta }\sum _{n=1}^N{\beta }_{s_{n-1}}\left(1-e^{-\delta \left(t_n-t_{n-1}\right)}\right)e^{-\delta \left(t_N-t_n\right)}\right],$$

$$\mathcal{G}\left(z,t|t\ge t_N\right)=$$

$${\mathcal{G}}^{\text{new}}\left(z,t\right):=$$

$$\psi \left({\left\{t_n\right\}}_{n=1}^N\right):=$$

$$e^{-\delta \left(t-t_N\right)}+{\beta }_{s_N}\left(1-e^{-\delta \left(t-t_N\right)}\right)]$$

$$\mathcal{G}^{\text{init}}\left(z,t\right)\equiv\sum_{m=0}^{\infty}\mathbb{P}_{m,0}\sum_{n=0}^{\infty}\left(\begin{array}[]{c}m\\ n\end{array}\right)\left(1-e^{-\delta t}\right)^{m-n}z^{n}e^{-n\delta t}.$$

$$\left[1-p+pz\right]^{m}=\sum_{n=0}^{\infty}\left(\begin{array}[]{c}m\\ n\end{array}\right)\left(1-p\right)^{m-n}z^{n}p^{n}.$$

$$\frac{\text{d}}{\text{d}t}\lambda \left(t\right)={\beta }_{s\left(t\right)}-\delta \lambda \left(t\right),\text{and }\lambda \left(0\right)=0.$$

$${\mathcal{G}}^{\text{new}}\left(z,t\right)=$$

$$=$$

$$\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t\right)>0\right)=$$

$$1-\sum _{m=1}^{\infty }{\mathbb{P}}_{m,0}{\left(1-e^{-\delta t}\right)}^m.$$

$$\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t\right)>0\right)$$

$$=$$

Full text

stochastic gene expression $|$ Bayesian statistical inference $|$ model selection $|$ piecewise-deterministic Markov process $|$ single-cell experiments

I Introduction

Gene expression is a biochemical process driven by the chance collisions of molecules, which can result in strong stochastic signatures and cell-to-cell variability in gene dynamics. Advances in single-cell and single-molecule technologies have provided unprecedented resolution on the stochastic dynamics of gene expression Pichon et al. (2018). Dynamic assays measure gene expression in living cells either directly via transcript tagging Golding et al. (2005); Larson et al. (2011); Corrigan et al. (2016); Fukaya et al. (2016) or indirectly via fluorescent or luminescent proteins Ozbudak et al. (2002); Elowitz et al. (2002); Suter et al. (2011); Nicolas et al. (2018). Static assays measure transcript levels in fixed cells either using a cocktail of fluorescently-labelled DNA oligos that bind specific transcripts Raj et al. (2008); Zenklusen et al. (2008) or via single-cell RNA sequencing Tang et al. (2009); Hwang et al. (2018). Static assays are popular because they do not require genetic modifications and are easily multiplexed. The disadvantage is that static assays only provide population snapshots of transcripts levels and cannot follow the dynamics of transcription in a single cell through time.

To this end, static assays have relied upon mathematical models to infer dynamic properties of gene expression in single cells from the measured snapshot of transcript levels; see Munsky et al. (2015) for a review. Inference requires (1) appropriate models of stochastic gene expression, (2) numerical methods to calculate the time-dependent mRNA distribution in a population of cells given any underlying model and associated parameters, and (3) calculating the likelihood that measured data were sampled from the calculated distribution. We recently developed a Bayesian approach (BayFISH) that uses this likelihood to infer best-fitting parameters from single cell data and quantify their uncertainty using the posterior distribution Gómez-Schiavon et al. (2017). Although Bayesian inference is the most complete and rigorous approach, it requires significantly more computation than other approximate methods, e.g. maximum likelihood.

We have developed a hybrid numerical method that accelerates the calculation of time-dependent mRNA distributions by 1000-fold compared to standard methods. We integrated this method into BayFISH and, for the first time, one can estimate the Bayesian posterior distribution of many competing models in a reasonable amount of time. The Bayesian evidence rigorously quantifies the likelihood of a model relative to other models given the data, and we show that HME selects the true model and outperforms approximate metrics, e.g., BIC or AIC, typically used for model selection. Our accelerated Bayesian inference represents a significant advance over existing methods used for inferring gene expression models from snapshots of single cell transcripts.

II Results

Our inference method uses data from single-molecule RNA Fluorescence In Situ Hybridization (smFISH), but could include single cell data from other static assays. The smFISH technique labels transcripts with fluorescent DNA oligos and measures both the number of mature mRNAs ($m$) and the number of gene loci with high-activity transcription sites ($TS$); see Figure 1a. A typical smFISH data set is a histogram $h=h(\vec{\omega})$ where $\vec{\omega}\in\Omega$ is the set of all possible states $(m,TS)$ that can be measured in a cell.

II.1 Connecting models of gene expression to single cell data

A broad spectrum of measured gene expression profiles in bacteria and eukaryotes is well-explained by discrete state gene expression models Raj and van Oudenaarden (2008); Nicolas et al. (2017), summarized by the following reactions:

[TABLE]

In this article, we adopt a two-allele, 3-state model (Figure 1b) as a case study for modelling stochastic gene expression in eukaryotes and for testing our method of accelerated Bayesian inference. We further focus on dynamic smFISH experiments that perturb gene expression (e.g. induction) and then measure mRNA distributions at different times after induction to infer dynamics and kinetic parameters. Induction can change one or more of the model parameters (Figure 1c). The smFISH data from an induction experiment consist of a joint histogram $h=h(\vec{\omega},t_{\ell})$, where $t_{\ell}$ are independent observations made at different times before and after induction. If the changed parameters are unknown a priori, then one should evaluate all possible induction models, which leads to a combinatorial explosion in model space. For example, there are $2^{8}=256$ candidate induction models for the 3-state model shown in Figure 1c.

A likelihood approach is used to connect mathematical models of stochastic gene expression to single cell data. Formally, the likelihood $\cal{L}$ is the probability that a candidate model $\mathcal{M}$ and its associated parameter set $\vec{\theta}$ would generate a given set of data ($h$). The number of parameters (i.e., dimension of $\vec{\theta}$) is determined by the model structure $\mathcal{M}$. Mathematically, the likelihood $\mathcal{L}$ is a function of the joint probability distribution $\mathbb{P}(\vec{\omega},t_{\ell}|\vec{\theta},\cal{M})$ of a candidate model $\mathcal{M}$ and its associated parameters $\vec{\theta}$ at discrete observation times:

[TABLE]

where $\Phi$ is the set of observation times and $\mathbb{M}_{\ell}$ is the Multinomial coefficient associated with each $h(\vec{\omega},t_{\ell})$ that arises because the data were not ordered.

In our Bayesian inference work flow (Figure 1d), each candidate model $\cal{M}$ in the class of possible models $\left\{\cal{M}\right\}$ will require a large number ($\geq 10^{6}$) of Monte Carlo steps where, at each step, numerical simulations calculate the time-dependent mRNA distributions and evaluate the likelihood that different parameter sets $\vec{\theta}$ for that model generated the observed data. Our previous software Gómez-Schiavon et al. (2017) took days to perform the likelihood calculations for one model, which highlights the challenge of using Bayesian inference to evaluate hundreds of models and perform model selection. Below, we develop a hybrid method that both accelerates numerical simulation and likelihood calculations, and (in contrast to standard methods) scales with the number of multi-core processors, thus, allowing for efficient parallelization.

II.2 A novel hybrid method to calculate the time evolution of discrete-state models

While exact time-dependent solutions exist for two-state models Shahrezaei and Swain (2008); Kumar et al. (2014); Tiberi et al. (2018), it is hard to generalize this analysis to models with more than two states. It is therefore necessary to solve the general time-dependent problem using numerical simulations. There are two classes of numerical procedures to solve the time evolution of a discrete-state model for a given set of parameters. The first class forward-evolves the chemical master equations (CME), which are a system of infinitely many coupled ordinary differential equations (ODEs) that describe the joint probability distributions $\mathbb{P}(\vec{\omega},t)$ as a function of time Munsky and Khammash (2006); Cao et al. (2016). To be numerically feasible, a truncation scheme (e.g. only consider mRNA levels below a maximum $M$) is used to reduce the infinite size of the dynamical system. While this approach delivers accurate estimates of the temporal evolution of the truncated system, there are two shortcomings. First, the number of ODEs scale as $S^{2}M$ where $S$ is the number of genetic states for each allele. The ODE system becomes unwieldy for mammalian cells where the number of observed mRNAs per cell can be $\mathcal{O}(10^{3})$ Raj et al. (2006); Little et al. (2013); Bahar Halpern et al. (2015); Skinner et al. (2016). Second, the forward integration of the CME requires stiff ODE solvers, which can place demands on memory resources and hinder parallel processing. The second class of numerical procedures utilize kinetic Monte Carlo methods (e.g., continuous time Markov chain simulation Gillespie (1976, 1977)) to sample the temporal evolution of the joint probability distribution $\mathbb{P}(\vec{\omega},t)$. While this approach is computationally less expensive, it comes at the cost of having to sample over many runs to achieve equivalent accuracy to CME.

In this article, we propose a hybrid simulation method (Poisson Mixture with Piecewise Deterministic Markov Switching Rates, or PM-PDMSR) which leverages analytical results and the efficiency of the kinetic Monte Carlo method. The key result is that the mRNA distribution can be exactly calculated for any realization (trajectory) of the genetic state, $s(t)$; see Methods. Once transient, initial conditions have burned off ($t\gg\delta^{-1}$), where $\delta$ is the mRNA degradation rate, the mRNA ($N_{\rm mRNA}$) distribution is always Poisson, $\mathbb{P}(N_{\rm mRNA}=m)=\lambda^{m}(t)e^{-\lambda(t)}/m!$ with a dynamic rate $\lambda(t)$ satisfying the following piecewise ODE:

[TABLE]

with an initial condition $\lambda(0)=0$. Given any trajectory $s(t)$, we can exactly compute the mRNA distribution $\mathbb{P}(m|s(t))$; see Figure 2a-b. Our goal, however, is to determine the joint distribution $\mathbb{P}(\vec{\omega},t)$, which requires us to generate $N_{s}$ sample paths of $s(t)$ that cover $\mathbb{P}(s,t)$. The sample paths in the small genetic state space ($S^{2}$-dimensional) are efficiently generated using standard kinetic Monte Carlo methods. After accumulating a large number of sample paths $N_{s}$ generated by the underlying model, the mixture of the Poisson distributions recovers the mRNA distribution:

[TABLE]

where where $\lambda_{k}\left(t\right)$ is the solution of (3) subject to the $k^{\text{th}}$ sample path of genetic switching trajectory $s_{k}\left(t\right)$ and $\delta_{i,j}$ is the Kronecker delta (see Figure 2c-e).

A detailed description of the hybrid simulator is given in the Methods section. We evaluated the efficiency of the hybrid simulator relative to CME in performing a single step of the Bayesian inference work flow, i.e. simulate the joint distribution $\mathbb{P}(\vec{\omega},t)$ and calculate the likelihood $\cal{L}$ that this joint distribution produced a given data set ($h$). We bench-marked the simulators on diverse classes of discrete-state models, parameter sets, and data sets; see Figure 3. The hybrid simulator is up to $10^{3}$ more efficient for models with increased genetic states, $S=3$ and $4$. The efficiency gain of the hybrid simulator originates from the fact that $\mathbb{P}(m|s(t))$ is solved exactly in mRNA space (and is independent of the size of $M$) and that $\mathbb{P}(s,t)$ is sampled efficiently in genetic-state space via kinetic Monte Carlo techniques. Last, we tested the scaling of efficiency of different simulators on a modern machine, which can execute 32 parallel threads. The hybrid method scaled linearly with the number of processors, i.e. 32 processors in parallel ran just as fast as one processor. Surprisingly, the simultaneous execution of 32 stiff ODE solvers for the CME simulator took 16 times longer per processor than one processor alone; see Methods.

II.3 Bayesian inference and uncertainty quantification of model parameters

Equipped with an efficient simulator of the time-dependent joint probability distribution and likelihood calculation for any model and parameter set, we first turned our attention to uncertainty quantification of model parameters $\vec{\theta}$ for a fixed model $\cal{M}$. Given a likelihood, Bayesian inference uses Bayes formula to update any prior beliefs $\mathbb{P}(\vec{\theta}|\cal{M}$) and calculate the posterior distribution $\mathbb{P}(\vec{\theta}|h,\cal{M})$ of parameters $\vec{\theta}$ given the data $h$ and a fixed model $\cal{M}$:

[TABLE]

As done previously, we resorted to Markov chain Monte Carlo (MCMC) with a Metropolis–Hastings (MH) sampler to estimate the posterior distribution $\mathbb{P}(\vec{\theta}|h,\cal{M})$; see Methods and Gómez-Schiavon et al. (2017). We assumed that the prior $\mathbb{P}(\vec{\theta}|\cal{M})$ is log-uniform. At each MCMC step, the MH sampler randomly proposes a nearby parameter set and computes the ratio of the posterior probability $\mathbb{P}(\vec{\theta}|h,\mathcal{M})$ relative to that of the current parameter set, and probabilistically accepts or rejects the proposal with a prescribed criterion that only depends on the ratio of the likelihood values. The denominator $\mathbb{P}(h|\mathcal{M})$ in (5) is identical for any parameter set $\vec{\theta}$ and cancels during the calculation of the ratio.

We bench-marked our approach on two synthetic data sets that were generated by sampling ($N=100$ or $1000$ cells per time point for a total of 4 time points) from a two-allele, 3-state induction model, where the induction stimulus decreased the downward transition rates; see Methods. Here, the model was known a priori and our goal was to infer the kinetic parameters and perform uncertainty quantification by comparing their posterior distributions (Figure 4a-b) to the ground truth (GT) parameters used to generate the sampled synthetic data set (Figure 4c-d). Our method constrained the posterior parameter distribution around the ground truth, and a 10-fold increase in the number of sampled cells dramatically reduced uncertainty in the inferred parameters. This observation holds true for a synthetic data set generated by a different two-allele, 3-state induction model; see Supplemental Materials.

II.4 Model selection using the full Bayesian framework

Knowing that our method of accelerated Bayesian inference can reliably constrain the kinetic parameters for a given model, we turned our attention to the harder problem of model selection. The goal was to identify the correct model from 64 possible types of two-allele, 3-state induction models given the same synthetic data set in Fig. 4, which was sampled from a ground-truth model and its parameters. Bayesian analysis naturally provides a quantitative measure of the likelihood of any model $\cal{M}$, i.e., the probability of the model to reproduce the experimentally observed data $h$. The measure, referred to as the marginalized likelihood or evidence MacKay (2005); Sivia and Skilling (2006), is the denominator of (5):

[TABLE]

The evidence is simply the probability that a model $\cal{M}$ produced data $h$ and is equal to the sum of the probabilities of the model (i.e. likelihood) over all sets of parameters that could have produced the data. The evidence is a convolution of the likelihood with the prior $\mathbb{P}(\vec{\theta}|\mathcal{M})$, which quantifies the belief regarding the initial parameter distributions. The dimensionality of $\vec{\theta}$ does not have to be identical for two different models and this prior inherently penalizes models with too many parameters; see Discussion.

The evidence for a model $\mathcal{M}$ is not calculated during the MCMC sampling of the posterior distribution and has to be computed separately. Computing the evidence is a sophisticated problem Green (1995); Kass and Raftery (1995); Weinberg (2012); Pajor (2017) and we adopted an Importance Sampler of the Harmonic Mean Estimator (IS-HME) proposed Robert and Wraith Robert and Wraith (2009), which re-samples the posterior distribution estimated by the MCMC to compute the evidence of each model; see Methods. We first carried out the MCMC calculations of posterior distributions for each of the 64 possible types of two-allele, 3-state induction models for the data sets described in Fig. 4. We then used IS-HME to compute the evidence of each model given the underlying data set. We compared the IS-HME evidence to maximum likelihood metrics used for model selection, such as the Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC) Gómez-Schiavon et al. (2017). Both BIC and AIC are approximations to the Bayesian evidence and become equivalent in the limit of large sample sizes; see Discussion.

Our results demonstrate that the IS-HME evidence $\mathbb{P}(h|\mathcal{M})$ of the ground truth model dominates over other models ($\geq 95\%$) when using Bayesian inference on the larger data set ($N=1000$ cells sampled per time point); see Figure 5. The BIC approximation also selected the ground-truth model (although incorrect models exhibited significant probabilities, e.g. $>5\%$) whereas the AIC failed to select the correct model. When the sample size dropped to $N=100$ cells per time point, even IS-HME evidence could not reliably select the ground-truth model with this under-powered data set.

III Discussion

Piecewise-deterministic Markov processes (PDMP) have become a useful, coarse-grained description of stochastic gene dynamics, where the underlying discrete variable $s(t)$ captures the stochastic dynamics of gene states and the continuous variable $\lambda(t)$ captures the first moment of downstream gene products Bokes et al. (2013); Lin and Doering (2016); Lin and Galla (2016); Hufton et al. (2016); Bressloff (2017); Lin et al. (2018); Lin and Buchler (2018). The key insight of our manuscript was proving that the time-dependent mRNA distribution of any underlying $s(t)$ is asymptotically a Poisson distribution with a rate $\lambda(t)$, and that the time-dependent joint probability distributions of discrete-state models are dynamic Poisson mixtures, whose mixing kernels are characterized by a PDMP. This significantly expands upon a related framework, which only considered the stationary distribution of discrete-state models Herbach (2018). More generally, our analysis helps bridge a gap between mechanistic discrete-state models and statistical models used in single cell analysis. For example, Wang et al. recently proposed an statistical model of gene expression, which postulated that mRNA distributions are Poisson mixtures Wang et al. (2018), and our work justifies this assumption.

We used our insight to develop a hybrid method that calculates the time-dependent joint distribution more efficiently than standard numerical methods that forward-integrate the Chemical Master Equation (CME). The efficiency arises because our method analytically solves the mRNA distribution and rapidly samples many path $s(t)$ of discrete-switching events using a kinetic Monte Carlo algorithm. We benchmarked the hybrid method and showed that it is $\mathcal{O}\left(10^{3}\right)$ more efficient than previous methods that directly integrate the CME. Furthermore, the hybrid method runs as efficiently in parallel on a multi-core processor than it does on a single processor. The stiff CME integrators ran more slowly in parallel and this sublinear scaling persisted for different integrators. We suspect that the slow-down arises from the competing memory demands of stiff CME integrators running on a multi-core processor. While there is room to improve stiff integration and parallelization, current approaches are limited compared to the hybrid method because they must integrate the CME for a large number of mRNA states, e.g. 0 to 1000 mRNAs per cell.

We incorporated the hybrid algorithm into BayFISH and were able, for the first time, to use a full Bayesian framework for model selection and uncertainty quantification of parameters from single-cell smFISH data. We adopted the Bayesian framework for model selection because it naturally quantifies “Occam’s factor” MacKay (2005); Kass and Raftery (1995) and, thus, avoids over-fitting. For example, the top models based on Bayesian evidence are not the most complex models with the largest number of parameters that change upon induction, e.g. $\mathcal{M}=11110110$; see Figure 5. The evidence resists over-fitting because when the dimensionality of parameter space is high, the value of a uniform probability density of the prior parameter distribution $\mathbb{P}\left(\vec{\theta}|\mathcal{M}\right)$ in (6) is small due to normalization. Thus, Bayesian evidence will favor a model that is complex enough to have a large likelihood but not too complex to decrease the prior parameter density. We note that when the data sample size is large, such that the posterior distributions $\mathbb{P}\left(\vec{\theta}|h,\mathcal{M}\right)$ can be approximated by a multivariate normal distribution, the logarithm of the evidence converges asymptotically to the Schwartz index (commonly known as the Bayesian Information Criterion, BIC) Kass and Raftery (1995); Sivia and Skilling (2006). Here, we bench-marked the ability of Bayesian evidence and the BIC metric to select the correct model from synthetic data sets generated by a ground-truth model and parameters. Figure 5 shows that while the BIC (but not AIC) ranked models similarly to Bayesian evidence for the larger data set ($N$=1000 cells per time point), BIC requires an even larger sample size to confidently converge to the correct model. This is an important issue because most biology labs are ill-equipped to generate and analyze large smFISH data sets, and their sample sizes are typically $N$ = 100 - 1000 cells per time point. Our work demonstrates why Bayesian inference should be used for modestly sampled data sets. We show that $N$=100 cells per time point is sufficient for parameter inference and uncertainty quantification if one has high confidence in the underlying model; see Figure 4. However, if the goal of the smFISH experiments is model selection, then these smaller data sets are under-powered and the experimentalist needs to increase data sampling by atleast 10-fold; see Figure 5a and 5c. Here, we only considered one round of experiments followed by Bayesian inference, but multiple cycles of data collection and analysis are the norm. Our framework quantifies certainty in both models and parameters using Bayesian evidence and posterior distributions. Future work can complete the data collection and analysis cycle by using the evidence and posterior distributions to rationally dictate the next round of experiments, i.e. different sampling times and densities, that are most informative for constraining models and/or parameters.

In this article, we adopted a Markov chain Monte Carlo algorithm with Metropolis-Hastings sampling (i.e., MCMC-MH) to compute the posterior distributions of the model parameters Gómez-Schiavon et al. (2017); Munsky et al. (2018). However, there is room to further improve the speed of Bayesian inference. First, Hamiltonian Monte Carlo (HMC) algorithms are more efficient at sampling posterior distributions in high dimensional parameter spaces because they uses local sensitivity, i.e., the partial derivatives of the likelihoods with respect to the model parameters Duane et al. (1987); Neal et al. (2011); Betancourt (2017). Second, although PM-PDMRS is efficient at generating sample paths in the space $\left(s,\lambda\right)$, evaluating the convolution to calculate the joint distribution (4) is the rate-limiting step in the likelihood calculation. Thus, transforming the experimental data $h$ into the mixing kernel of the Poisson mixtures $\rho_{s}\left(\lambda\right)$ would accelerate Bayesian inference. Last, one could use low-order moments of PM-PDMSR and experimental data to formulate a sufficient statistics for likelihood-free Approximate Bayesian Computation Toni et al. (2009), thus, replacing the explicit calculation of the likelihood $\mathcal{L}$.

IV Methods

IV.1 Poisson mixture with piecewise deterministic Markov switching rates

We illustrate our central theoretical results using a single-allele model. However, the results generalize to multiple-allele models because the states of a multiple-allele model can be relabeled as internal states of a single-allele model.

Central theoretical result I: Given a trajectory of genetic state $s(t)$, the total number of mRNA, $N_{\text{mRNA}}\left(t\right)$, is the sum of two variables $N_{\text{mRNA}}^{\text{initial}}\left(t\right)$ and $N_{\text{mRNA}}^{\text{new}}\left(t\right)$. $N_{\text{mRNA}}^{\text{initial}}\left(t\right)$ describes the number of initial mRNAs that remain at time $t$. The probability distribution of $N_{\text{mRNA}}^{\text{initial}}\left(t\right)$ is a Binomial mixture weighted by the initial mRNA distribution $\mathbb{P}_{m,0}:=\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t=0\right)=n\right)$:

[TABLE]

Here, $\Theta\left(\cdot\right)$ is the Heaviside step function. $N_{\text{mRNA}}^{\text{new}}\left(t\right)$ describes the number of new mRNAs that were synthesized after $t>0$ and still remain at time $t$. The probability distribution of $N_{\text{mRNA}}^{\text{new}}\left(t\right)$ is a Poisson distribution with rate $\lambda\left(t\right)$

[TABLE]

where $\lambda(t)$ satisfies equation $\dot{\lambda}\left(t\right)=\beta_{s\left(t\right)}-\delta\lambda\left(t\right)$ with an initial condition $\lambda\left(0\right)=0$.

Proof. We denote the probability of the total number of mRNA $N_{\text{mRNA}}$ at time $t$ by $\mathbb{P}_{m}\left(t\right)$. For a given trajectory of the genetic state $s(t)$, the temporal evolution of satisfies the chemical master equation (CME)

[TABLE]

for all $m\in\mathbb{Z}_{\geq 0}$ and with a boundary condition $\mathbb{P}_{-1}=0$. We prove the central theoretical result I by using the probability generating function defined by

[TABLE]

where $\mathbb{P}_{m}\left(t\right)=\partial_{z}^{m}\mathcal{G}\left(z,t\right)/m!$. We first solve for $\mathcal{G}\left(z,t\right)$ over an interval of time when $s\left(t\right)$ is constant. We will then extend our analysis to include piecewise intervals of time with different value of constant $s$, similar to the $s\left(t\right)$ generated by a genetic state model. To begin, we apply the operator $\partial_{t}$ to $\mathcal{G}\left(z,t\right)$ and use (12) to obtain the partial differential equation (PDE):

[TABLE]

This linear PDE can be solved using the method of characteristics Kampen (2007) and the general solution is

[TABLE]

where $\mathbb{P}_{m,0}:=\mathbb{P}_{m,0}\left(t=0\right)$ is the initial mRNA distribution of the system. For reasons that will become apparent below, we label the initial-distribution-dependent part by $\mathcal{G}^{\text{init}}\left(z,t\right)$ and the rest of the terms by $\mathcal{G}^{\text{new}}\left(z,t\right)$:

[TABLE]

The above solution applies to a constant $s\left(t\right)$. We now consider a piecewise-constant trajectory for any genetic state $s\left(t\right)$. To specify the discrete state and the switching events, we label $s\left(t\right)$ by the ordered pairs $\left(t_{i},s_{i}\right)$ for $i=0,1,\ldots,N$ before an observation time $t$. The gene starts with a state $s_{0}$ at $t_{0}:=0$, switches to $s_{1}$ at $t_{1}$,…until the final switching event to $s_{N}$ at time $t_{N}$. Our aim is to compute the generating function $\mathcal{G}\left(z,t\right)$ after $N$ switching events ($t\geq t_{N}$).

The solution $\mathcal{G}\left(z,t,|t\leq t_{1}\right)$ is identical to (15) with $s=s_{0}$ before the first switching event. At $t=t_{1}$, the generating function is

[TABLE]

Note that after $t_{1}$ and before $t_{2}$, the genetic state changes to $s_{1}$ and only the transcription rate in (12) changes from $\beta_{s_{0}}$ to $\beta_{s_{1}}$. The initial condition of the generating function of this period ($t_{1}\leq t\leq t_{2}$ ) is precisely $\mathcal{G}\left(z,t_{1}\right)$ in the above equation. Matching the “initial condition” for $\mathcal{G}\left(z,t_{1}\right)$, we arrive at

[TABLE]

We iterate this procedure for each piecewise “episode” until $t=t_{N}$,

[TABLE]

and for $t\geq t_{N}$,

[TABLE]

The total solution $\mathcal{G}\left(z,t|t\geq t_{N}\right)$ is factorized into two terms, $\mathcal{G}^{\text{init}}\left(z,t\right)$ and $\mathcal{G}^{\text{new}}\left(z,t\right)$, for any $N$ and $t$. The probability generating function of the sum of two independent random variables is the product of the generating functions of the random variables. This hints that we can define two random variables, $X_{1}$ and $X_{2}$, which have generating functions $\mathcal{G}^{\text{init}}\left(z,t\right)$ and $\mathcal{G}^{\text{new}}\left(z,t\right)$ respectively.

Our next task is to identify the variables $X_{1}$ and $X_{2}$ and their probability distributions. For $X_{1}$, we expand $\mathcal{G}^{\text{init}}\left(z,t\right)$ to arrive at

[TABLE]

Recall that the generating function of a $\text{Binomial}\left(m,p\right)$ distribution is

[TABLE]

The probability distribution of $X_{1}$ is therefore identified to be a binomial mixure with a temporally decaying parameter $p=\exp\left(-\delta t\right)$ and a mixing kernel defined by the initial distribution $\mathbb{P}_{m,0}$. The physical meaning of $X_{1}\left(t\right)$ is the number of initial mRNA molecules that remain at time $t$, i.e., $N_{\text{mRNA}}^{\text{initial}}\left(t\right)$. These mRNA can only degrade with the decay rate $\delta$. Each of the mRNA decays independently and, at time $t$, there is a probability $\exp\left(-\delta t\right)$ that a specific mRNA has not degraded. Importantly, when $t\gg 1/\delta$, this distribution will be concentrated at $m=0$ (see Corollary I below).

The total mRNA $N\left(t\right)=X_{1}\left(t\right)+X_{2}\left(t\right)=N_{\text{mRNA}}^{\text{initial}}\left(t\right)+X_{2}\left(t\right)$, so $X_{2}\left(t\right)$ is identified to be the number of new mRNA molecules that were synthesized after $t=0$ but which have not degraded at time $t$. We refer to this variable as $N_{\text{mRNA}}^{\text{new}}\left(t\right)$. The squarer bracket of $\mathcal{G}^{\text{new}}\left(z,t\right)$ in (20) is the piecewise solution $\lambda\left(t\right)$ of the following ODE for a given genetic trajectory $s\left(t\right)$:

[TABLE]

We now expand $\mathcal{G}^{\text{new}}\left(z,t\right)$:

[TABLE]

where $q_{m}\left(\lambda\left(t\right)\right)$ is the probability density function of a Poisson distribution with rate $\lambda\left(t\right)$, as in (11).$\ \square$

Corollary I. The transient timescale for the initial distribution is $\mathcal{O}\left(1/\delta\right)$. When $t\gg\mathcal{O}\left(1/\delta\right)$, the mRNA distribution converges to a Poisson with a dynamic rate parameter $\lambda\left(t\right)$.

Proof. Physically, the timescale of degradation of each initially populated mRNA is $1/\delta$, so for a timescale which is much longer than this, the initial distribution will be fully degraded. Mathematically, the probability that the initial mRNA molecules have not fully decayed is

[TABLE]

In the asymptotic limit $t\gg 1/\delta$, $\exp\left(-\delta t\right)\ll 1$ so

[TABLE]

Here, $\left\langle N_{\text{mRNA}}^{\text{initial}}\left(0\right)\right\rangle$ is the first moment of the initial distribution. $\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t\right)>0\right)$ decays exponentially fast, and we can bound this probability to be smaller than $\varepsilon$ when $t>\delta^{-1}\log\left(\left\langle N_{\text{mRNA}}^{\text{initial}}\left(0\right)\right\rangle/\varepsilon\right)$. $\ \square$

For our numerical simulations detailed below, our total time was $20/\delta$. The contribution of the initial distribution $\mathbb{P}\left(N_{\text{mRNA}}^{\text{initial}}\left(t\right)>0\right)=\left\langle N_{\text{mRNA}}^{\text{initial}}\left(0\right)\right\rangle e^{-20}\approx 2\times 10^{-9}\times\left\langle N_{\text{mRNA}}^{\text{initial}}\left(0\right)\right\rangle\approx 10^{-6}$, because typical numbers of measured mRNAs are bounded by $\mathcal{O}\left(500\right)$.

Central theoretical result II: At long times $t\gg 1/\delta$, the mRNA distribution asymptotically converges to a Poisson mixture regardless of the initial mRNA distribution and genetic switching trajectory $s\left(t\right)$

[TABLE]

where the joint probability density $\rho_{i}\left(\lambda,t\right)$ satisfies the forward Kolmogorov equation

[TABLE]

The initial condition for $\rho_{i}\left(\lambda,t=0\right)$ are defined by

[TABLE]

where $\delta\left(\lambda\right)$ is the Dirac delta distribution at $\lambda=0$.

Proof. The solution of (23) subject to random switching events of $s\left(t\right)$ is a random process. Formally, $\lambda\left(t\right)$ and the discrete switching states $s\left(t\right)$ jointly comprise a piecewise deterministic Markov process (PDMP) Davis (1984); Bena (2006); Faggionato et al. (2009). The forward Kolmogorov equation describing the temporal evolution of the joint probability distribution is (28) Hufton et al. (2016); Lin et al. (2018). Therefore,

[TABLE]

We then use the central theoretical result I and corollary I to show that $\mathbb{P}\left(N_{\text{mRNA}}\left(t\right)=m|\lambda\left(t\right)=\ell,s\left(t\right)=i\right)=\ell^{m}e^{-\ell}/m!$ asymptotically when $t\gg 1/\delta$ to complete the proof. $\ \square$

Efficient numerical method for sampling $\rho_{i}\left(\lambda,t\right)$. Because $\lambda$ is continuous, solving the forward Kolmogorov equation (28) is as complex as solving the full CME, both of which are infinite dimensional systems. Instead, we used an efficient kinetic Monte Carlo simulation Lin et al. (2018); Bokes et al. (2013) to generate a large number of sample paths to estimate the asymptotic joint distribution $\mathbb{P}\left(N_{\text{mRNA}}\left(t\right)=m|\lambda\left(t\right)=\ell,s\left(t\right)=i\right)$ when $t\gg 1/\delta$ using (27). The pseudo code of this algorithm is shown in Algorithm 1.

Computing the joint distribution from sample paths. To compute the time-dependent joint distribution of genetic states and mRNAs, we generated $N_{s}$ sample paths $\left\{\lambda_{k}\left(t\right),s_{k}\left(t\right)\right\}_{k=1}^{N_{s}}$ with Algorithm I. We then used (27) to estimate $\rho_{i}\left(\lambda\right)$

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta function which is equal to $1$ if $i=j$, and [math] otherwise. We determined that $N_{s}\equiv 10^{5}$ is a sufficient number of sample paths to estimate the same joint distribution obtained by forward-integrating the CME. We refer to our method as the Poisson Mixture with a Piecewise Deterministic Markov Switching Rate (PM-PDMSR).

The goal was to compute the joint distribution of genetic states and mRNAs before and after induction. Similar to the situation in many experimental systems Neuert et al. (2013); Senecal et al. (2014), the joint distribution is at stationarity before induction. Upon induction, we assume some model parameters are changed, which results in the time-evolution of the joint probability distribution $\mathbb{P}\left(N_{\text{mRNA}}\left(t\right)=m,s\left(t\right)=i\right)$ towards a new stationary state. We label the kinetic rates ($\kappa_{ij}$ and $\beta_{k}$) before and after induction by $U$ (Unstimulated) and $S$ (Stimulated). To use PM-PDMSR to estimate the stationary distribution before induction, we first solved for the marginal stationary distribution of the genetic state $p_{i}^{\ast}$, where $0=\sum_{j}\left(\kappa_{ij}^{U}-\kappa_{ji}^{U}\right)p_{i}^{\ast}$, $i=1,2,\ldots,S$ for an $S$-state model. We initiated $N_{s}p^{\ast}_{i}$ sample paths in PM-PDMSR at $\lambda=\beta_{i}/\delta$ and state $s=i$ and ran for $t=10/\delta$ so that the Poisson mixture relaxes to stationarity. Upon induction at $t=10/\delta$, we changed model parameters $\kappa_{ij}^{U}\rightarrow\kappa_{ij}^{S}$ and $\beta_{k}^{U}\rightarrow\beta_{k}^{S}$, for $i,j,k\in\left\{1,2,\ldots S\right\}$ and continued simulating the temporal evolution of the joint probability distribution after induction using PM-PDMSR. This is valid because our previous proofs and arguments for (27) apply even when the kinetic rate constants change upon induction.

IV.2 Testing the speed and accuracy of the simulators

We bench-marked the efficiency of PM-PDMSR versus the “gold-standard” simulator, i.e., forward-integration of a truncated CME Munsky and Khammash (2006); Cao et al. (2016). Both simulators were embedded into BayFISH and evaluated on their ability to perform a single Monte Carlo step, i.e., simulate the time-dependent joint distribution of a model and its parameters and to calculate the likelihood of a synthetic data set generated by the same model and parameters. This single-step bench-marking was performed for 1024 diverse models and parameters across two-allele, discrete-state models of increasing complexity (2-state, 3-state, and 4-state induction models).

Generating diverse models, parameters, and synthetic data sets. Each model had a randomly chosen subset of parameters that were one value $\left(\kappa_{ij}^{U},\beta_{k}^{U}\right)$ before induction from $t=0$ to $t=20$ and with different parameters $\left(\kappa_{ij}^{S},\beta_{k}^{S}\right)$ after induction from $t=20$ to $t=22$. The genetic switching rates of parameters $\kappa_{ij}$ ranged between $10^{-2}$ to $10^{2}$ and the transcription rates $\beta_{k}$ range between [math] and $200$. For each instance of a model, we randomly generated the switching rate constants in the logarithmic space ($log_{10}\kappa_{ij}\sim\text{Unif}\left(-2,2\right)$ if $|i-j|=1$) and transcription rate constants in linear space ($\beta_{k}\sim\text{Unif}\left(0,200\right)$). Note that not all parameters were random or changed upon induction: We fixed $\beta_{0}=0$ and the mRNA degradation rate $\delta=1$, and we constrained $\beta_{i}\leq\beta_{j}$ if $i<j$. For the purpose of bench-marking algorithm speed and efficiency, we considered the most complex induction model $\mathcal{M}$ for each discrete model class. The corresponding induction model is $\mathcal{M}=11010$ (2-state), $\mathcal{M}=11110110$ (3-state), and $\mathcal{M}=11111101110$ (4-state).

The synthetic data of each model with its parameters were generated by running 1000 independent trajectories of standard continuous-time Markov chain simulation. We collected the statistics of the trajectories at four discrete times $t_{\ell}=20$, $20.5$, $21$, and $22$ ($N=1000$ cells per time point). The measured allele activity state $TS$ was marginalized: we define $TS=1$ when its internal state $s>0$. The synthetic data set therefore consists of a histogram at discrete times $t_{\ell}$, $h\left(m,TS,t_{\ell}\right)$, which is the number of trajectories with $m$ mRNAs and $TS$ active transcription sites (which can be [math], $1$, or $2$ for a two-allele system) at time $t_{\ell}$. For each model class (2-, 3-, and 4-state models), we repeated the process 1024 times to test diverse parameter combinations.

CME simulators. Given an induction model and associated parameters, we forward propagated the CME using the same parameter values used to create the synthetic data sets in the previous section. We truncated the number of mRNA at $500$ (i.e., there are no transcription event once the system reaches $N_{\text{mRNA}}=500$) with an absolute error tolerance of $10^{-5}$. The truncation number $M$ was motivated by data sets in animals, showing that mRNA populations can be as large as $\mathcal{O}\left(500\right)$ Raj et al. (2006); Little et al. (2013); Bahar Halpern et al. (2015); Skinner et al. (2016). We tested different software platforms, including Matlab, Python (SciPy), and a research software ACME Cao et al. (2016) to forward integrate the same stiff CME. Python’s stiff integrator (using backward differentiation formula, BDF) outperformed other integrators and software platforms. Thus, Python (with SciPy) was chosen to be the platform for direct integration of the CME in the following analysis.

Comparison of the PM-PDMSR and CME simulators. We incorporated these PM-PDMSR and CME simulators into modified BayFISH software to evaluate their speed and accuracy for one Monte Carlo step. Both algorithms were implemented in c++ and compiled using Intel’s icc compiler. All PM-PDMSR and CME simulations were carried out on the same machine with Intel(R) Xeon(R) CPU E5-2695 v3 @ 2.30GHz. We computed the joint distributions of the models and parameter sets used to create the synthetic data sets. These joint distributions and corresponding synthetic data ($h\left(m,TS,t_{\ell}\right)$) were then used to compute the likelihood $\mathcal{L}$ of the generated data $h\left(m,n_{TS},t_{\ell}\right)$ using (2). The execution time of a Monte Carlo step for each simulator for each model class of the generated parameter set was recorded and presented in Fig. 3. We also compared the accuracy of the calculated likelihoods of each synthetic data set. We compared the average error of the PM-PDMSR likelihood relative to the more accurate CME likelihood. We define the average error by

[TABLE]

The relative accuracy of CME versus PM-PDMSR is presented in Supplemental Materials.

To test the parallelization of each simulator, we simultaneously ran 32 simulations on a 32-core machine (Intel(R) Xeon(R) CPU E5-2695 v4 @ 2.10GHz) and recorded the execution time; see Fig. 3 PM-PDMRS on 32 parallel threads takes the same amount of time as running a single thread on a 32-core machine. In contrast, CME on 32 parallel threads takes 16 times longer than a single thread on a 32-core machine. This sub-optimal scaling holds true on the multiple machines that we tested and our algorithms leverage Python’s sub-processes functionality (Pool). We suspect that the slow down in CME is due to the high memory demand of the BDF integrator.

IV.3 Bayesian statistical inference for model parameters and structure

Synthetic data. We synthesized data sets to test if Bayesian statistical inference could identify the ground truth (of the model parameter values and the model structure.) We chose two 3-state ground-truth models for data synthesis: (1) a ON-rate induction model, $\kappa_{01}^{U}=1\rightarrow\kappa_{01}^{S}=12$, $\kappa_{12}^{U}=0.25\rightarrow\kappa_{12}^{S}=20$, $\kappa_{10}^{U}=\kappa_{10}^{S}=3$, $\kappa_{21}^{U}=\kappa_{21}^{S}=10$, $\beta_{0}=0$, $\beta_{1}=25$, and $\beta_{2}=300$, and (2) an OFF-rate induction model, $\kappa_{01}^{U}=\kappa_{01}^{S}=1$, $\kappa_{12}^{U}=\kappa_{12}^{S}=1$, $\kappa_{10}^{U}=10\rightarrow\kappa_{10}^{S}=0.1$, $\kappa_{21}^{U}=10\rightarrow\kappa_{21}^{S}=0.1$, $\beta_{0}=0$, $\beta_{1}=100$, and $\beta_{2}=200$. The protein degradation rate constant is $\delta=1$ by choosing the timescale of the model. We relaxed the models from $t=0$ to $t=20$, and sampled the system at $t=20$, $20.5$, $21$, and $22$. The observation timescale was motivated by our recent experimental procedure Gómez-Schiavon et al. (2017). For each model, we synthesized by sampling $100$ and $1000$ synthetic data at each of the descrete sampled times from the joint probability distribution. The data consists of the sampled and discrete number of mRNA, and whether the gene is active. Again, the genetic space $s$ is marginalized that we defined $s>0$ is an active allele ($TS=1$), otherwise inactive ($TS=0$).

Bayesian analysis. The model class we considered for Bayesian inference is the set of two-allele, 3-state models with $\beta_{0}^{U}=\beta_{0}^{S}=0$. The rest of the parameters are free parameters, but depending on the model structure, some of the perturbed ($S$) parameters may be constrained to the unperturbed ($U$) value. Combinatorially, there are in total $2^{6}=64$ models we considered, as there are six biophysical parameters $\vec{\theta}:=\left(\kappa_{01},\kappa_{12},\kappa_{21},\kappa_{10},\beta_{1},\beta_{2}\right)$.

We adopted a plain Markov chain Monte Carlo algorithm to sample the posterior distribution $\mathbb{P}\left(\vec{\theta}|h,\mathcal{M}\right)$, using the Metropolis-Hastings sampler Metropolis et al. (1953); Hastings (1970); Newman and Barkema (1999). Specifically, we perform random jumps in the logarithm space of the parameters ($\log\kappa_{ij}$ and $\log\beta_{k}$). The jump kernel was chosen to be uniformly distributed $\text{Unif}\left(-D,D\right)$, where the meta-parameter $D$ (diffusivity) globally regulates how wide the isotropic diffusion kernel is. For each model structure, we adjusted the meta-parameter $D$ such that the acceptance rate of the Metropolis-Hastings sampler was between $0.2$ and $0.3$ Roberts et al. (1997). We assumed that our prior is uniform in the logarithm space $\log\theta_{i}$.

Before running the MCMC samplers, we randomized $400$ initial guesses of the model parameters and forward evolved the MCMC for $5\times 10^{4}$ iterations each chain. Most of the chains converged to a unique parameter region and the likelihood value in this region was significantly higher than the few chains trapped in (presumably) local maxima. We independently initiated $32$ MCMC chains with different random initial seeds from the parameter values that maximized the likelihood in the previous test runs. We collected a total length (sum of the length of all $32$ chains, $>10^{7}$) to accurately approximate the posterior distribution $\mathbb{P}\left(\vec{\theta}|h,\mathcal{M}\right)$.

Computing the evidence $\mathbb{P}\left(h|\mathcal{M}\right)$. As described in Robert and Wraith (2009), the evidence $\mathbb{P}\left(h|\mathcal{M}\right)$ is computed by the algebraic identity

[TABLE]

with an importance sampler $\phi\left(\vec{\theta^{\prime}}\right)$ satisfying the normalization condition

[TABLE]

In this work, the posterior distributions were exclusively unimodal. Therefore, we chose the importance sampler to be proportional to an indicator function on an ellipsoid located at the high posterior density region. We first ranked the posterior chain by their likelihood and selected the top 20% parameter sets to construct the importance sampler. We performed a Principle Component Analysis on the selected samples to compute the mean $\bar{\theta_{k}}$, variance $\sigma_{k}^{2}$ and eigenvector $\hat{e}_{k}$, $k=1,2,\ldots 6$, in the six-dimensional parameter space. We used the eigenvalues and eigenvectors to construct an ellipsoid centering at the mean and with the axis length proportional to the eigenvalues along with the eigenvectors:

[TABLE]

where $R_{i,j}:=\left(\hat{e}_{i}\right)_{j}$ is the linear transformation onto the eigenbasis. We tuned the meta-parameter $\alpha$ such that there were precisely 20% of the points of the posterior chains inside the ellipsoid $\mathcal{E}$. These samples were then used to compute the evidence.

One must also specify the prior distribution $\mathbb{P}\left(\vec{\theta}|\mathcal{M}\right)$ to compute the evidence. For simplicity, we imposed a uniform prior in the logarithm space of the parameters, bounded by $\left(10^{-16},10^{4}\right)$. $\mathbb{P}\left(\vec{\theta}|\mathcal{M}\right)$is $1/V$, where $V$ is the bounded volume in the parameter space. Let the posterior chains to be $\left\{\vec{\theta}_{k}\right\}_{k=1}^{N_{P}}$, where $N_{P}$ is the total number of samples in the posterior chains. Then, the marginalized likelihood is estimated computed by the Monte Carlo sampler:

[TABLE]

where $\mathbf{1}_{\vec{\theta}_{k}\in\mathcal{E}}$ is the characteristic function which is equal to $1$ if $\vec{\theta}_{k}$ is in the ellipsoid $\mathcal{E}$, and [math] otherwise. In Figure 5, we presented a normalized probability among all the 64 linear 3-state models we considered:

[TABLE]

which is reported in Fig. 5.

The Schwarz Index is an asymptotic result of the Bayesian evidence $\mathbb{P}\left(h|\mathcal{M}\right)$ when the sample size is large Schwarz (1978). Given a model $\mathcal{M}_{i}$, the Bayesian Information Criterion (BIC) is defined to be twice of the its Schwarz Index

[TABLE]

where $\mathcal{L}_{i}^{\max}$ and $m_{i}$ are maximum likelihood and the number of free parameters of model $\mathcal{M}_{i}$ respectively, and $N$ is the sample size (number of data). Thus, to estimate the normalized probability $\bar{\mathbb{P}}$ using BIC, we use

[TABLE]

We remark that the calcaulation of BIC only involves estimating the maximum likelihood $\mathcal{L}_{i}^{\max}$ of each model $\mathcal{M}_{i}$ and not the full posterior distribution $\mathbb{P}\left(\vec{\theta}|\mathcal{M}_{i}\right)$.

Akaike Information Criterion (AIC) is another commonly adopted matrix information criterion. The motivation of AIC is to minimize the information loss, measured by the Kullback–Leibler divergence (KL) of the prediction from the data. It is derived Akaike (1998) as

[TABLE]

Thus, the (normalized) evidence calculated by the AIC is

[TABLE]

We remark that the negative logarithm of the likelihood function (2) converges to $N\times\text{KL}\left(h\left(\omega\right)/N||\mathbb{P}\left(\omega|\mathcal{M}_{i},\vec{\theta}\right)\right)$ only when $N\gg 1$ such that the Multinomial coefficient $\mathbb{M}_{\ell}$ can be expanded by the Stirling approximation. In most biological cases, the sample size is far from this regime, and the Kullback–Leibler divergence is a poor choice to approximate the likelihood function (2).