Determining noisy attractors of delayed stochastic Gene Regulatory Networks from multiple data sources

. ABSTRACT Motivation: Gene regulatory networks (GRN) are stochastic, thus, do not have attractors, but can remain in conﬁned regions of the state space, the ‘noisy attractors’, which deﬁne the cell type and phenotype. Results: We propose a gamma-bernoulli mixture model clustering algorithm ( Γ BMM), tailored for quantizing states from gamma and bernoulli distributed data, to determine the noisy attractors of stochastic GRN. Γ BMM uses multiple data sources, naturally selects the number of states and can be extended to other parametric distributions according to the number and type of data sources available. We apply it to protein and RNA levels, and promoter occupancy state of a toggle switch and show that it can be bistable, tristable, or monostable depending on its internal noise level. We show that these results are in agreement with the patterns of differentiation of model cells whose pathway choice is driven by the switch. We further apply Γ BMM to a model of the MeKS module of Bacillus subtilis , and the results match experimental data, demonstrating the usability of Γ BMM. Availability: Implementation software is available upon request.


INTRODUCTION
In the deterministic framework, it was hypothesized that cell types are attractors of gene regulatory networks (GRN) dynamics (Kauffman, 1969).To be biologically plausible, the attractors need to be confined patterns of gene expression since even in the simple case of a Boolean network with 20 000 nodes, its state space size is 2 20000 (Kauffman, 1993).An important criticism to this hypothesis (Aldana et al., 2002) is that noise renders attractors a poor model of cell types.This is important given the stochasticity of gene expression (Arkin et al., 1998).To address this, the concept of 'ergodic set' was proposed as a set of states from which the GRN, once entering, does not leave for a time longer than the cell life time in the absence of external perturbations (Ribeiro et al., 2007).An ergodic set can be composed of one or several 'noisy attractors', the set of states where the GRN spends most time when on that ergodic set.Stochastic models of GRN (Ribeiro et al., 2006) can remain in confined regions of the state space for long periods of time, in agreement with the hypothesis that ergodic sets corresponds to cell types (Ribeiro et al., 2007).One example of noise-driven transitions between noisy attractors might be the observed transient and probabilistically transition between two phenotypes in Bacillus subtilis, controlled by the MeKS module (Suel et al., 2006).
While techniques to find noisy attractors and ergodic sets in deterministic Boolean networks are known (Ribeiro et al., 2007), in stochastic GRN this is not a simple task.In (Ribeiro et al., 2007), the 'states' composing the noisy attractors of a stochastic GRN were found by binarizing the proteins' temporal expression levels with Kmeans clustering algorithm (MacQueen, 1967).This is not the most appropriate approach since the number of clusters was predefined, possibly causing information loss.Stochastic GRN's 'logic' can be more complex than binary (e.g., the two genes of a toggle switch (TS) can express at half its maximum rate, simultaneously, for long periods of time (Ribeiro, 2008)).Also, observing the proteins' level alone might not be sufficient to fully capture the behavior of GRN (e.g., miRNA's are known to have a regulatory role (Zhang et al., 2007)).Hence, one problem that needs to be addressed is how to define noisy attractors of a stochastic GRN in a realistic biological setting.For this the key questions are what features of the GRN dynamics should be taken into account to determine whether a GRN is in a noisy attractor, and how to reliably identify noisy attractors from the GRN dynamics given the choice of features to consider.
In the delayed stochastic model of GRN (Ribeiro et al., 2006), able to match gene expression at the single molecule level (Zhu et al., 2007;Yu et al., 2006), what changes with time are the number of proteins and RNA, and genes' promoter occupancy state if transcription factors can bind to their target genes' promoters and affect transcription.Here we use these three sources of information to determine noisy attractors.While proteins (Elowitz et al., 2002) and RNA levels (Golding et al., 2005) can be measured in real time, promoters occupancy states are not directly measurable currently.However, it can be detected indirectly in vivo when a transcription factor is bound to a gene's promoter region (King et al., 2007;Thompson et al., 2004), which is one assumption made here.
Inspired by the joint mixture model clustering framework used in gene clustering (Dai et al., 2008), we built a gamma-bernoulli mixture model, ΓBMM, to quantize and detect noisy attractors of stochastic GRNs from multiple data sources.The algorithm accounts gamma and bernoulli distributed data simultaneously and quantizes their joint effort into state(s), assuming that proteins and RNAs' levels are the summation (including subtraction) of several exponential processes, i.e., gamma distributions, and that the promoter occupancy state, which is composed of binary values, is of bernoulli distribution.Each distribution for a gene is modeled as an independent event, and the quantization is done based on the multiple levels to make the detection of changes in the dynamics of the GRN more robust than when observing protein levels alone.The number of noisy attractors are determined by model selection, where four approximation-based criteria are compared: the Akaike information criterion (AIC) (Akaike, 1974), a modified AIC (AIC3) (Biernacki et al., 1999;Bozdogan, 1987), Bayesian information criterion (BIC) (Pan, 2006;Schwarz, 1978), and integrated classification likelihood-BIC (ICL-BIC) (Ji et al., 2005).BIC and ICL performed better and similarly, hence are employed here.
To define noisy attractors as biologically meaningful entities in cell's dynamics, we assume that by being in a certain noisy attractor the cell will opt for a specific differentiation pathway, or adopt a phenotype out of several possible choices (Kauffman, 1969).We improve previous work (Ribeiro et al., 2007) by introducing a more robust clustering algorithm that uses data fusion to detect noisy attractors from the GRN dynamics and by validating its results by adding a stochastic model of cell differentiation.
In this study, we apply ΓBMM to the TS and to a GRN that drives differentiation in Bacillus subtilis, the MeKS module (Suel et al., 2006), and validate the results by confronting them with the cell types distributions resultant from differentiation.The results show that the clustering algorithm is robust to high levels of noise in GRN dynamics by using multiple data sources to determine the noisy attractors, and that the modeling strategy provides a good match with experimentally measured differentiation patterns.The results provide a better understanding of how stochastic GRNs can regulate their dynamics to determine cell populations differentiation patterns and how multi-modal distributions of phenotypic expression can emerge from the stochastic nature of gene expression and regulation.

Mixture model based clustering and EM algorithm
In model-based clustering methods, each observation o j , where j = 1, . . ., n and n is the number of genes, is drawn from a finite mixture distribution with the prior probability π i , component-specific distribution f (g) i and its parameters θ i .The formula is (Mclachlan et al., 2000) f where θ = {(π i , θ i ) : i = 1, . . ., g} denotes all unknown parameters, with the restriction that 0 < π i ≤ 1 for any i and g i=1 π i = 1.Note that g is the number of components in this model.From here onwards, we ignore the superscript (g) from f (g) i for simplicity.In ΓBMM, we define θ = [π, θ 1 , θ 2 ] T , π = [π 1 , . . ., π g ] T , θ 1 = [α 11 , . . ., α gp 1 , β 11 , . . ., β gp 1 ] T , and θ 2 = [q 11 , . . ., q gp 2 ] T , where p 1 and p 2 each represents the dimension of the observations in gamma and bernoulli mixture models, respectively.X and Y denote the observations of gamma and bernoulli distributed data, respectively, function f of x, y as the density function of gamma and bernoulli distribution, and o = [x T , y T ] T .
ΓBMM is a joint mixture model of gamma and bernoulli distributions that assumes that for each component i the data of the two distributions are independent.Each component is assumed to be the product of p 1 independent gamma distributions, whose probability density function is defined as where Each component is assumed to follow a bernoulli distribution, with the probability density function for each gene defined as (θ 2i = [q i1 , . . ., q ip 2 ]) (2) The expectation maximization (EM) algorithm is applied to estimate the parameters θ iteratively.Eq. 3 is the data log-likelihood (natural logarithm is referred throughout the paper), given O = {o j : j = 1, ..., n} (its direct maximization, however, is difficult): To make the maximization of Eq. 3 tractable, the problem is casted in the framework of incomplete data.Since we assume that data from different distributions are independent, Lc can be factored as If we define c j ∈ {1, . . ., g} as the clustering membership of o j , then the complete data log-likelihood can be written as where χ(c j = i) is the indicator function of whether o j is from the i th component or not.In the EM algorithm, E step computes the expectation of the complete data log-likelihood where θ (m) represents the parameters estimated in the m th iteration.By computing the expectation, Eq. 4 becomes, according to Bayes' rule, where τ , We can assign each o j to its component based on {i 0 |τ ji 0 = max i τ ji }.
Eqs. 4 and 5 show that our assumption of the gamma and bernoulli distributed data being independent carries over to the expected log-likelihood as well.To derive the closed form or numerical optimization formula for updating parameters in ΓBMM, we used lagrange multipliers to solve this constrained optimization problem, with the lagrangian shown in Eq. 6.

Downloaded from
Determining noisy attractors with ΓBMM Parameters of ΓMM part, α and β, can be optimized by Newton-Raphson method.Let θ 1i = (α i , β i ), then the parameters are updated by and L(θ H −1 and ∇ θ 1 L(θ) can be obtained by plugging (1) in ( 6) taking the derivatives: and then forming .
Ψ and Ψ represent the digamma and trigamma functions, respectively, which are the first and second logarithmic derivatives of the gamma function.
Parameter updating of q is done by plugging Eq. 2 in Eq. 6, and taking the derivative of q Optimization of the prior probability of each gene's clustering membership, π, can be derived by and since

Model selection criteria
Four approximation-based model selection criteria, BIC (Pan, 2006;Schwarz, 1978), ICL (Ji et al., 2005), AIC (Biernacki et al., 1999;Akaike, 1974), and AIC3 (Biernacki et al., 1999;Bozdogan, 1987), are compared in ΓBMM, from which the best-performing one is chosen.Calculations for the above criteria are defined as where d is the number of free parameters, and M is the total amount of data (M = W w=1 M w , M w is the size of the data set w and W is the number of input data sets).(−2 n j=1 g i=1 τ ji log(τ ji )) is the estimated entropy of the fuzzy classification matrix C ji = (τ ji ) (Ji et al., 2005).The number of free parameters d in ΓBMM is d RB = 2p 1 g + p 2 g + g − 1 since there are p 1 g free α iu 's, p 1 g free β iu 's, p 2 g free q iv 's, and (g − 1) free π i 's.

Data generation
We use SGN Sim (Ribeiro and Lloyd-Price, 2007) to simulate the models of GRNs.GRNs consist of genes coupled via protein-protein and proteinoperator sites interactions.Transcription and translation are multiple time delayed events (Ribeiro et al., 2006).The dynamics is driven by the 'delayed SSA' (Roussel et al., 2006), a modified version of the Stochastic Simulation Algorithm (SSA) (Gillespie, 1977).The GRN models here studied (and how to implement them in SGN Sim) are given in supplementary material.
The time series of proteins, RNAs and promoter binding states are used as the inputs of ΓBMM.The outputs are the number of noisy attractors and the noisy attractor membership and log-likelihood of each data point (Fig. 1).

RESULTS AND DISCUSSION
We first compared the performance of ΓBMM with its component models, ΓMM and BMM, using controllable gamma and bernoulli distributed data (data and results in supplementary material).After observing an advantage in using multiple data sources compared with using single data sources alone, we used ΓBMM to detect the noisy attractors of two GRNs, a TS (Gardner et al., 2000) and the MeKS module, which is an excitable genetic circuit of Bacillus subtilis (Suel et al., 2006).

The Toggle Switch
The TS consists of reactions 7 to 12 (Ribeiro, 2008), where i ∈ {1, 2} (when only index i is present), or i, j ∈ {1, 2} with i = j (when both indices are present).Transcription (reaction 7) and translation (reaction 8) are modeled by time-delayed reactions (Ribeiro et al., 2006) where P ro i is the promoter of gene i, and R p is an RNA polymerase.The delays (τi) account for the time duration of the processes involved in transcription and translation and model how long it takes for a product to be produced after the reaction has occurred.Each promoter P ro i controls the expression of an RNA, Ri, i.e., once its ribosome binding site is formed, can start being translated by a ribosome Rib into a protein, P i .Reaction 11 models the binding and unbinding of the repressor protein to the other gene's promoter, defining the TS.Proteins decay via reactions 9, 10 and 12.
Values for delays (fixed in all simulations) are extracted from measurements (Yu et al., 2006).τ1 accounts for the transition from the closed promoter complex to the open promoter complex which is included in the transcription initiation process (McClure, 1980).We set τ 1 = 40 s, within the range from 10 s to several minutes (McClure, 1980).The length of the gene tsr-venus driven by a Lac promoter studied in E. coli is 2500 nucleotides.Since the average rate of transcriptional elongation in E. coli is ∼ 50 nt/s, it follows that τ2 is, on average, 90 s (= τ1+2500/50).The time of the ribosome binding site clearance in translation initiation, τ 3 , is set to 2 s, according to experimental estimations (Draper, 1996).The average translation rate is ∼15 amino acids/s, thus τ 4 = τ 3 + 2500 nt/(45 nt/s) = 58 s.The post-translational protein assembly process takes 420 ± 140 s in (Yu et al., 2006), thus the delay of Pi is drawn from a Gaussian distribution with a mean of τ 5 and standard deviation of τ 5std each time the reaction occurs.All delays, except for τ5, are assumed constant (Zhu et al., 2007) since measurements suggest that they are not highly variable between transcription events of a single gene (McClure, 1980).For example, τ1, although varies from one transcription event to another, measurements on an unrepressed Lac promoter (Lutz et al., 2001) suggest that it follows a Gaussian distribution with a mean of 40 s and a standard deviation of 4 s.We tested the effects of variable delays following a normal distribution with the same mean and comparably small standard deviations (with respect to the mean), and the TS dynamics was not significantly affected, which means, qualitatively, our conclusions hold when using a model with variable delays.
Reaction rates (in s −1 ) are set following the TS model in (Zhu et al., 2007), built to reproduce the dynamics of the genetic TS synthesized in E. coli (Gardner et al., 2000).Specifically, k tr = 0.00042, rbsd = 0.01, krep = 1, and kunrep = 0.1.Given these values, the transcription rate k t and protein decay rate k d are varied as discussed ahead.
The noise level of the TS can be varied by multiplying the transcription rate kt and protein decay rate k d with a pair of factors, while imposing equal mean values of (P 1 + P 2 ) to ensure that the noise varies from varying the toggling rate alone (Ribeiro, 2008).Both these rates are evolvable and vary from gene to gene.The initiation rate is sequence dependent while degradation rates are tunable and vary from one protein to another (Bernstein et al., 2002).The noise level affects the TS bistability.Excess of noise makes the TS monostable, with both genes expressing simultaneously (if transcription is more 'propense' than the reaction by which the promoter is repressed) (Ribeiro, 2008).Thus, the TS model allows testing the clustering algorithm for various levels of noise.
To test the ability of ΓBMM in quantizing genes' expression levels we simulate nine 'cell populations'.Each cell consists of a TS model.Cells of different populations have TSs with identical mean level of proteins but different noise level which is defined as the standard deviation over the mean of (P1 + P2) time series.Each cell population has a unique pair of multiplication factors {ν1, ω1, . . ., ν9, ω9} of the values of kt and k d , respectively.These multiplication factor pairs cannot be varied by constant amounts each time to attain equal mean proteins' levels due to the delay τ 1 at each transcription event, which accounts for the promoter occupancy time by an Rp (McClure, 1980).We first fix ν's, where ν = 1 is the 'baseline' and four inverse pairs are arranged on each side in ascending order; then ω's are searched by an heuristic method such that the mean value of (P1 + P2) doesn't differ more than 4% among cells of different populations.The nine pairs of multiplicative values of k t and k d are (0.01,0.2768), (0.1, 1.0935), (0.25, 1.25), (0.5, 1.2), (1, 1), (2, 0.9226), (4, 1.08), (10, 1.73), and (100, 3.4599), respectively.Originally, k t = 0.01 and k d = 0.001.
We generate 1000 cells for each population with SGN Sim.Given 10 6 s as the cell lifetime we sampled the state with a sampling frequency of 1000 s.This unrealistic lifetime allows better statistics without affecting the results qualitatively (alternatively one could follow the dynamics of cells for many generations).For each population, we quantize protein, RNA and promoter binary binding data of a cell by ΓBMM, and use the log-likelihood (Eq. 3) of each cluster within each noisy attractor and their summation to visualize the number of noisy attractors and the ergodic set's stability.
Note that for this study we simulate nine cell populations, each contains 1000 identically initialized cells, and cells of different population differ.The variation of parameter values across cell populations allows testing the robustness of the clustering algorithm under various degrees of GRN's internal noise.The simulation of many initially identical cells, each with the same internal parameter values and amount of each substance, is necessary to estimate the GRN's differentiation pattern for a given set of parameters' values, given the stochastic dynamics of the genetic networks driven by the delayed SSA.These patterns are then confronted with the quantization results for validation purposes.
Log-likelihood distribution of the ergodic set for nine cell populations are shown in Fig. 2 (a), where each log-likelihood is the sum of each state of each noisy attractor.The higher the loglikelihood is the more stable the quantized states are.Also shown is the noise in cells of the nine populations (Fig. 2 (b)).The , respectively.'NA' is short for 'noisy attractor'.'Ensemble' is the ensemble of the time series of proteins, RNAs and promoter binding state.
x-axis is time (unit is ×10 3 ), and y-axis is log-likelihood of each quantized data point.
TS noise increases as k t and k d increase due to the increase of toggling frequency, and then decreases due to the loss of bistability (Ribeiro, 2008).The turning point of log-likelihood (population 4, Fig. 2 (a)) corresponds to the loss of bistability and maximum noise (Fig. 2 (b)).
The quantization results of each noisy attractor are shown in Supplementary Fig. S2, with the most representative ones (populations 1, 3, 4 and 9) shown in Fig. 3.Note that there are two noisy attractors (of the same ergodic set) in population 1 (Fig. 3 (a)), which starts toggling as noise increases (Fig. 3 (b)), begins loosing its bistability in cells of population 4 (Fig. 3 (c)), and becomes monostable beyond this point (Fig. 3 (d)).Importantly, Fig. 3 (c) shows that cells of population 4 can be in three noisy attractors.The magenta and blue dots each represent one noisy attractor, and the red and green circles together represent a third noisy attractor.To observe changes in the number of attractors of the TS as noise increases we draw a transition graph (Fig. 4) by defining that two noisy attractors are distinct (and not a single noisy attractor) if the difference between their medians is larger than the mean of their standard deviations.The noisy attractors, defined as the set of states where the GRN spends most time when on an ergodic set, correspond to the states where sufficient data points share similar log-likelihood in Fig. 4. The result, which agrees with Fig. 3, shows that as noise increases, the TS changes from bistable to tristable and to monostable, verifying the claim in (Ribeiro, 2008) that the TS dynamical regime changes with noise level, and showing the suitability of ΓBMM to characterize the dynamics of stochastic GRNs.
To verify the existence of these regimes we consider the TS as a differentiation pathway switch (Monod and Jacob, 1961) by adding in model cells reactions 13 to 16.These are activated at the end of a cell's lifetime and the amount of P1 and P2 in the cell at that moment will stochastically determine the pathway ('Cell 0 ' represents an undifferentiated cell).For example, if only P 1 is present when these reactions become active, the cell will most likely differentiate into cell type 1 via reaction 13, if P 1 and P 2 are both present with significant amounts, Cell0 will most likely differentiate into cell type 3 (reaction 15), and if both proteins are absent, the cell chooses pathway 4 (reaction 16).
Fig. 5 shows the distribution of differentiation pathway choices of 1000 cells of four of the nine populations (complete results in Supplementary Fig. S3).As seen in Figs. 5 (b) and (c), cells of populations 4 and 5 each have three almost equally likely differentiation pathways, which accord well with the three noisy attractors found in cells of these populations (Fig. 4).Also in agreement, cells differentiate mostly into one or two type(s) when the TS only has one or two noisy attractor(s), respectively (Figs. 5(d) and (a)).
While the bistable and monostable regimes are visually detectable from the proteins' time series, the tristable regime is not, as seen in Supplementary Fig. S4 (c) and (d), where there is no visible difference between the TS dynamics of cells of populations 3 and 4 although the transition from bistable to tristable occurs at that point (Fig. 3 (b) and (c)).This is because the 'state', where both proteins express simultaneously, is unstable in the deterministic setting.Interestingly, for a given noise level, this unstable state is the most commonly occupied region of the state space by the TS (Fig. 5 (c)).

An excitable circuit of Bacillus subtilis
Under nutrient limitations, a minority of the cells of a population of B. subtilis becomes competent for DNA uptake while the majority commits to sporulation (Suel et al., 2006).Cells can switch probabilistically and transiently between competence and vegetative growth via a 'MeKS' module (Fig. 6) (Suel et al., 2006).In this module, protein ComK activates genes needed for competence.Protein ComS induces competence by allowing accumulation of ComK protein that indirectly enables the self-activation of gene comK by inhibiting ComK degradation by protein complex MecA.We model the 'MeKS' module in the delayed stochastic modeling framework following the ODE model proposed in (Suel et al., 2006) (both models are described in detail in the supplementary material).The rate constants and the initial amounts of all substances are heuristically tuned to match both the ODE model dynamics and the measurements of the time series of proteins P S and P G reported in (Suel et al., 2006).Time delays are set as previously described for the transcription and translation reactions in the TS since such information is not available for this genetic circuit.
Each cell is initialized with 100 MecA, 1 P ro i and no P i or R i , where 'i' stands for 'S', 'K' and 'G' for proteins, and represents 'S' and 'K' for RNAs and promoters.We simulated 10 3 cells for 100 hours and sampled each time series at each 10 3 seconds.The time series of proteins ComS and ComG of one cell is shown in Fig. 7 (a), where a negative correlation is apparent between them, in agreement with the experimental and modeling results in (Suel et al., 2006).
To validate the model we add reactions 18 and 17 to each cell to drive differentiation.Sporulation (incompetent cells) depends on the amounts of proteins ComS and ComG, and all cells differentiate during their life time.Cell S and Cell G represent the incompetent and competent cells, respectively.
Fig . 7 (b) shows the distribution of competent and incompetent cells (1000 cells in total).4% of the cells are competent, as in the experiments (3.6 ± 0.7%) (Suel et al., 2006).We applied ΓBMM to protein, RNA, and promoter binding time series of comS and comK since these are the genes regulating the transient differentiation.The quantization of one of the time series is shown in Fig. 8, where two noisy attractors are observed (in this case the competent state takes 4.04% of the sampled data).From a sample of 20 cells, 3.06% − 4.52% of the time is spent by the cells on the competent state, in the range of the experimental results, i.e., 2.9% − 4.3% (Suel et al., 2006), and in accord with the fraction of competent cells in the differentiation model.

CONCLUSION
In pluricellular organisms different cell types have identical GRNs.These, although being stochastic (Arkin et al., 1998)   to external perturbations, can stably maintain their type once fully differentiated.Reverse differentiation and metaplasias, while possible, are rare, indicating that cell types are ergodic sets.If the hypothesis of ergodic sets as cell types holds true then, from the dynamics of GRN, one ought to be able to predict patterns of differentiation and distributions of phenotypic expression of cell populations.
Here, we used a clustering method, ΓBMM, to detect noisy attractors of GRN from their dynamics.ΓBMM uses multiple time series of the GRN, and its applications to the TS and the MeKS module matched the patterns of cell differentiation.Data fusion improved the matching, suggesting that proteins levels alone might be insufficient to characterize the dynamics of stochastic GRNs.
The method is applicable to the analysis of real cells since the data required is experimentally attainable and the modeling strategy of GRNs (Ribeiro et al., 2006) used to test it matches gene expression measurements at the single RNA and protein levels (Zhu et al., 2007;Yu et al., 2006) and the model of differentiation of the MeKS module matched the measurements as well.While 1000 time series were used to test the accuracy of the algorithm, once tested, this amount of data is not needed in subsequent applications.The algorithm can then be applied to real time series directly, requiring only moderate amounts of experimental data.

Fig. 2 .Fig. 3 .
Fig. 2. Log-likelihood and noise distributions for nine toggle switches.xaxis represents the nine cell populations, y-axis shows the (a) summation of the log-likelihood of all the data points and (b) noise of each population.

Fig. 4 .
Fig. 4. Noisy attractors detected by ΓBMM in cells of each population.The number of bars in each of the nine populations equals the number of different noisy attractors detected (x-axis).The summation of the log-likelihood of the data points in each noisy attractor is shown in the y-axis.

Fig. 7 .
Fig. 7. Validation of the MeKS model.(a) Time series of ComS and ComG in one cell.(b) Distribution of differentiated cell types (of a total of 1000 cells).x-axis is (a) time (unit is ×10 5 ) and (b) cell types after differentiation.y-axis are (a) protein expression levels, and (b) number of cells.

Fig. 8 .
Fig.8.Quantization results from ΓBMM for the MeKS module.time in x-axis, unit is ×10 3 , and log-likelihood of each data point in y-axis.