Multi-modality in gene regulatory networks with slow promoter kinetics

TL;DR

This paper develops a rigorous mathematical framework to analyze phenotypic variability in gene regulatory networks with slow promoter kinetics, providing explicit formulas for metastable states and insights into how network parameters influence cell states.

Contribution

It introduces a novel analytical approach using singular perturbation theory to characterize multiple modes in stochastic gene regulatory networks with slow promoter kinetics, surpassing previous numerical methods.

Findings

Explicit formulas for metastable states in GRNs

Insights into the impact of network parameters on phenotypic states

Application to toggle switches, oscillators, and trans-differentiation models

Abstract

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.

Tables

Table 1. Table 1: Comparing the number of stable equilibria/modes for a self-regulating gene between stochastic with slow promoter kinetics and deterministic modelling frameworks. Details for reconstructing the table are given in SI-§5.3.

	Non-Cooperative		Cooperative
	Leaky	Non-Leaky	Leaky	Non-Leaky
Stochastic	2	1	2	1
(Slow promoter kinetics)		(at 0)		(at 0)
Deterministic	1	1	1-2	1-2

Equations

$${\mathrm{D}}_0^i\underset{{\alpha }_{-i}}{\overset{{\alpha }_i}{\rightleftharpoons }}{\mathrm{D}}_1^i.$$

$$\mathrm{TF}+{\mathrm{D}}_0^i\underset{{\alpha }_{-i}}{\overset{{\alpha }_i}{\rightleftharpoons }}{\mathrm{D}}_1^i,$$

$${\mathrm{TF}}_1+{\mathrm{D}}_{00}^i$$

$${\mathrm{TF}}_1+{\mathrm{D}}_{01}^i$$

$${\mathrm{TF}}_2+{\mathrm{D}}_{00}^i$$

$${\mathrm{TF}}_2+{\mathrm{D}}_{10}^i$$

$${\mathrm{D}}_j^i\stackrel{k_{ij}}{\to }{\mathrm{D}}_j^i+{\mathrm{X}}_i,$$

$${\mathrm{X}}_i\stackrel{k_{-i}}{\to }\varnothing .$$

$$n{\mathrm{X}}_i\underset{{\beta }_{-i}}{\overset{{\beta }_i}{\rightleftharpoons }}{\mathrm{X}}_{ic}.$$

$$R_j(z_i)=k_j\left(\genfrac{}{}{0ex}{}{z_i}{{\alpha }_{ij}}\right)=k_j\frac{z_i(z_i-1)...(z_i-{\alpha }_{ij})}{{\alpha }_{ij}!},$$

$$R_j(z_{i_1},z_{i_2})=k_j{z}_{i_1}{z}_{i_2}.$$

$$p_{d,x}(t)=\mathrm{Pr}[X(t)=x,D(t)=d],$$

$$p_d(t):=[p_{d{x}_0}(t),p_{d{x}_1}(t),....{]}^T,$$

$$p(t):=[p_0(t{)}^T,...,p_{L-1}^T(t){]}^T$$

$$\bar{\pi }=\underset{t\to \infty }{\lim }p(t).$$

$$\dot{p}(t)$$

$$\tilde{\Lambda }=\left[\begin{array}{ccc}{\Lambda }_0& & \\ & \ddots & \\ & & {\Lambda }_{L-1}\end{array}\right],\text{and}p(t)=\left[\begin{array}{c}{p}_0(t)\\ ⋮\\ p_{L-1}(t)\end{array}\right],$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\dot{p}_{X|d}=\Lambda_{d}p_{X|d},}$$

$$\varnothing \underset{k_{-i}}{\overset{k_{i{d}_i}}{\rightleftharpoons }}{\mathrm{X}}_i,$$

$$\varnothing \underset{k_{-i}}{\overset{k_{i{d}_i}}{\rightleftharpoons }}{\mathrm{X}}_i,n_i{\mathrm{X}}_i\underset{{\beta }_{-i}}{\overset{{\beta }_i}{\rightleftharpoons }}{\mathrm{X}}_{ic}.$$

$${\pi }_{X|d}^{(J)}(x)=\prod _{i=1}^N{\pi }_{X|di}(x_i),$$

$$\pi^{(J)}_{X|d_{i}}(x_{i})=\left\{\begin{array}[]{ll}\displaystyle\mathbf{P}\left(x_{i1},x_{i2};\frac{k_{id_{i}}}{k_{-i}},\frac{k_{id_{i}}^{n_{i}}\beta_{i}}{n_{i}!k_{-i}^{n_{i}}\beta_{-i}}\right)&\mbox{if}\ \mathrm{X}_{i}\ \mbox{is multimerized}\\ \mathbf{P}\left(x_{i};\frac{k_{id_{i}}}{k_{-i}}\right),&\mbox{otherwise}\end{array}\right.,$$

$${\pi }_{X|d}(x):=\sum _{i=1}^{M-N}\sum _{x_{i2}=0}^{\infty }{\pi }_{X|d}^{(J)}(x)=\prod _{i=1}^N\mathbf{P}\left(x_i;\frac{k_{i{d}_i}}{k_{-i}}\right),$$

$${\bar{\pi }}^{\epsilon }={\bar{\pi }}^{(0)}+{\bar{\pi }}^{(1)}\epsilon +o(\epsilon ).$$

$$\bar{\pi }(x,d)=\sum _{d=0}^{L-1}{\lambda }_d{\bar{\pi }}_{X|d}(x,d),$$

$${\Lambda }_r:=\left[\begin{array}{cccc}{1}^T& 0^T& ...& 0^T\\ 0^T& 1^T& ...& 0^T\\ & & \ddots & \\ 0^T& 0^T& ...& 1^T\end{array}\right]\hat{\Lambda }[{\bar{\pi }}_{X|0}{\bar{\pi }}_{X|1}...{\bar{\pi }}_{X|L-1}],$$

$${\pi }^{(J)}(x):=\underset{\epsilon \to 0^{+}}{\lim }\underset{t\to \infty }{\lim }\mathrm{Pr}[X(t)=x]=\sum _{d=0}^{L-1}{\lambda }_d{\pi }_{X|d}^{(J)}(x),$$

$$\pi (x):=\sum _{d=0}^{L-1}{\lambda }_d{\pi }_{X|d}(x)=\sum _{d=0}^{L-1}{\lambda }_d\prod _{i=1}^N\mathbf{P}\left(x_i;\frac{k_{i{d}_i}}{k_{-i}}\right).$$

$$\mathrm{TF}+{\mathrm{D}}_{d_i}^i\stackrel{\alpha }{\to }{\mathrm{D}}_{d_{i^{\prime }}}^i,$$

$${\mathrm{D}}_{d_i}^i\stackrel{\alpha \mathbb{E}[\text{TF}|D=d]}{\to }{\mathrm{D}}_{d_{i^{\prime }}}^i.$$

$$\mathbb{E}[\text{TF}|D=d]=\frac{\alpha }{n_i!}\frac{{\beta }_i}{{\beta }_{-i}}{\left(\frac{k_{i{d}_i}}{k_{-i}}\right)}^{n_i}.$$

$${\mathrm{D}}_0$$

Full text

Keywords: Gene Regulatory Networks $|$ Multi-modality $|$ Singular Perturbations $|$ Slow gene binding $|$ Slow Promoter Kinetics $|$ Markov chains $|$ Master Equation $|$ Cooperativity

A gene regulatory network (GRN) consists of a collection of genes that transcriptionally regulate each other through their expressed proteins. Through these interactions, including positive and negative feedback loops, GRNs play a central role in the overall control of cellular life [1, 2, 3]. The behavior of such networks is stochastic due to the random nature of transcription, translation, and post-translational protein modification processes, as well as the varying availability of cellular components that are required for gene expression. Stochasticity in GRNs is a source of phenotypic variation among genetically identical (clonal) populations of cells or even organisms [4], and is considered to be one of the mechanisms facilitating cell differentiation and organism development [5]. This phenotypic variation may also confer a population an advantage when facing fluctuating environments [6, 7]. Stochasticity due to randomness in cellular components and transcriptional and translational processes have been thoroughly researched [8, 9].

The fast equilibration of random processes sometimes allows stochastic behavior to be “averaged out” through the statistics of large numbers at an observational time-scale, especially when genes and proteins are found in large copy numbers. In those cases, an entire GRN, or portions of it, might be adequately described by a deterministic model. Stochastic effects that occur at a slower time scale, however, may render a deterministic analysis inappropriate and might alter the steady-state behavior of the system. This paper addresses a central question about GRNs: how many different “stable steady states” can such a system potentially settle upon, and how does stochasticity, or lack thereof, affect the answer? To answer this question, it is necessary to understand the possibly different predictions that follow from stochastic versus deterministic models of gene expression. Indeed, qualitative conclusions regarding the steady-state behavior of gene expression levels in a GRN are critically dependent on whether a deterministic or stochastic model is used (see [10] for a recent review). It follows that the mathematical characterization of phenomena such as non-genetic phenotype heterogeneity, switching behavior in response to environmental conditions, and lineage conversion in cells, will depend on the choice of the model.

In order to make the discussion precise, we must clarify the meaning of the term “stable steady state” in both the deterministic and stochastic frameworks. Deterministic models are employed when molecular concentrations are large, or if stochastic effects can be averaged out. They consist of systems of ordinary differential equations describing averaged-out approximations of the interactions between the various molecular species in the GRN under study. For these systems, steady states are the zeroes of the vector field defining the dynamics, and “stable” states are those that are locally asymptotically stable. The number of such stable states quantifies the degree of “multi-stability” of the system. Stochastic models of GRNs, in contrast, are based upon continuous-time Markov chains which describe the random evolution of discrete molecular count numbers. Their long-term behavior is characterized by a stationary probability distribution that describes the gene activity configurations and the protein numbers recurrently visited. Under weak ergodicity assumptions, this stationary distribution is unique [11], so multi-stability in the sense of multiple steady states of the Markov chain is not an interesting notion. A biologically meaningful notion of “multi-stability” in this context, and the one that we employ in our study, is “multi-modality,” meaning the existence of multiple modes (local maxima) of stationary distributions.

Intuitively, given a multi-stable deterministic system, adding noise may help to “shake” states, dislodging them from one basin of attraction of one stable state, and sending them into the basin of attraction of another stable state.

Therefore, in the long run, we are bound to see the various deterministic stable steady states with higher probability, that is to say, we expect that they will appear as modes in the stationary distribution of the Markov chain of the associated stochastic model. This is indeed a typical way in which modes can be interpreted as corresponding to stable states, with stochasticity responsible for the transitions between multiple stable states [12]. However, new modes could arise in the stationary distribution of a stochastic system besides those associated with stable states of the deterministic model, and this can occur even if the deterministic model had just a single stable state. This phenomenon of “stochastic multi-stability” has attracted considerable attention lately, both in theoretical and experimental work [8, 9, 13, 14, 15]. Stochastic multi-stability has been linked to behaviors such as transcriptional bursting/pulsing [16, 17] and GRN’s binary response [18]. Furthermore, multi-state gene transcription [4] has been used to propose explanations for phenotypic heterogeneity in isogenic populations.

A common assumption in gene regulation models is that transcription factor (TF) to gene binding/unbinding is significantly faster than the rate of protein production and decay [1]. However, it has been proposed [9, 19] that the emergence of new modes in stochastic systems in addition to those that arise from the deterministic model might be due to low gene copy numbers and slow promoter kinetics, which means that the process of binding and unbinding of TFs to promoters is slow. Thus, the emergence of multi-modality may be due to the slow TF-gene binding and unbinding. Already in prokaryotic cells, where DNA is more accessible to TF binding than in eukaryotic cells, some transcription factors can take several minutes to find their targets, comparable or even higher than the time required for gene expression [20],[21]. This is more relevant in eukaryotic cells, in which transcriptional regulation is often mediated by an additional regulation layer dictated by DNA methylation and histone modifications, commonly referred to as chromatin dynamics. For example, the presence of nucleosomes makes binding sites less accessible to TFs and therefore TF-gene binding/unbinding is modulated by the process of chromatin opening [22],[23, 9, 24, 25]. DNA methylation, in particular, has also been reported to slow down TF-gene binding/unbinding [26]. Several experiments have consolidated the role of the aforementioned complex transcription processes in slow promoter kinetics [17, 27, 28, 26].

In summary, new modes may appear in the stationary distribution that do not correspond to stable states in the deterministic model. Conversely, multiple steady states in the deterministic model may collapse, being “averaged out” by noise, with a single mode representing their mean. It is a well-established fact that, in general, multi-stability of the deterministic description of a biochemical network and multi-modality of the associated stochastic model do not follow from each other [29]. This is especially true in low copy number regimes with slow promoter kinetics. Figure 1 gives two examples for the emergence of new modes due to slow promoter kinetics, and it shows that equilibria derived from the corresponding deterministic model do not provide relevant information on the number and locations of the modes.

Here, we pursue a mathematical analysis of the role of slow promoter kinetics in producing multi-modality in GRNs and show analytically how the shape of the stationary distribution is dictated by key biochemical parameters. Previous studies of the chemical master equation (CME) for single genes has already observed the emergence of bimodality with slow TF-gene binding/unbinding [31, 19, 32, 33]. This phenomenon was also studied by taking the limit of slow promoter kinetics using the linear noise approximation [34] or hybrid stochastic models of gene expression [35, 36]. The basic constitutive gene expression model (refer to (Basic Example: Gene Bursting Model)) has been validated for transcriptional bursting [17]. However, and despite its application relevance, mathematical analysis of the CME for multi-gene networks with slow promoter kinetics has been missing, and only numerical solutions of the CME have been reported [37, 38].

In this work, an underlying theoretical contribution is the partitioning of the state space into weakly-coupled ergodic classes [11] which, in the limit of slow binding/unbinding, results in the reduction of the infinite-dimensional Markov chain into a finite-dimensional chain whose states correspond to “promoter states”. In this limit, the stationary distribution of the network can be expressed as a mixture of Poisson distributions, each corresponding to conditioning the chain on a certain promoter configuration. The framework proposed here enables us to analytically determine how the number of modes, their locations, and weights depend on the biophysical parameters. Hence, the proposed framework can be applied to GRNs to predict the different phenotypes that the network can exhibit with low gene copy numbers and slow promoter kinetics.

The results are derived by introducing a new formalism to model GRNs with arbitrary numbers of genes, based on continuous-time Markov chains. Then, we analyze the stationary solution of the associated CME through a systematic application of the method of singular perturbations [39]. Specifically, we study the slow promoter kinetics limit by letting the ratio of kinetic rate constants of the TF-gene binding/unbinding reactions with respect to protein reactions approach zero. The stationary solution is computed by applying the method of singular perturbations to the CME.

In order to illustrate the practical significance of our results, we work out several examples, some of which have not been studied before in the literature. As a first application, we discover that, with slow promoter kinetics, a self-regulating gene can exhibit bimodality even with non-cooperative binding to the promoter site. We then investigate the role of cooperativity. In contrast to deterministic systems, we find that cooperativity does not change the number of modes. Nevertheless, cooperativity adds extra degrees of freedom by allowing the network to tune the relative weight of each mode without changing its location.

As a second application, we revisit the classical toggle switch, under slow TF-gene binding/unbinding. It has been reported before that, with fast TF-gene binding/unbinding, the toggle switch with single-gene copies can be “bistable” without cooperative binding [40]. We show that this can also happen with slow promoter kinetics, and, moreover, that a new mode having both proteins at high copy numbers can emerge. We provide a method to calculate the weight of each mode and show that the third mode is suppressed for sufficiently high kinetic rates for the dimerization reactions.

A third application that we consider is a simplified model of synchronization of communicating toggle switches. In bacterial populations, quorum sensing has been proposed [41] as a way for bacterial cells to broadcast their internal states to other cells in order to facilitate synchronization. Quorum sensing communication has been adopted also as a tool in synthetic biology [42, 43]. Mathematical analysis of coupled toggle switches designs usually employs deterministic models [44]. We study a simplified stochastic model of coupled toggle switches with slow promoter kinetics and compare the resulting number of modes with deterministic equilibria.

Our final, and potentially most significant, application is motivated by cellular differentiation. A well-known metaphor for cell lineage specification arose from the 1957 work of Waddington [45], who imagined an “epigenetic landscape” with a series of branching valleys and ridges depicting stable cellular states. In that context, the emergence of new modes in cell fate circuits is often interpreted as the creation of new valleys in the epigenetic landscape, and (deterministic) multi-stability is employed to explain cellular differentiation [5]. However, an increasing number of studies have suggested stochastic heterogeneous gene expression as a mechanism for differentiation [13, 46, 47]. Numerical analysis of the CME for the canonical cell-fate circuit have shown the emergence of new modes due to slow promoter kinetics in such models [37, 48]. This general category of cell-fate circuits includes pairs such as PU.1:GATA1, Pax5:C/EBP$\rm\alpha$ and GATA3:T-bet [49]. Cell fate circuits are characterized by TF cross-antagonism. However, their behavior is affected by the promoter configurations available for binding, the cooperativity index of the TFs, and the relative ratio of production rates. Hence, we study two models that differ in the aforementioned aspects and we highlight the differences between our findings and the behavior predicted by the corresponding deterministic model. The first model employs independent cooperative binding. We show that such a network can exhibit more than four modes. In contrast, the deterministic model predicts up to four modes only with cooperativity [50] . The second network is a PU.1/GATA.1 network which employs non-cooperative binding and a restricted set of promoter configurations. The deterministic model is monostable, while the parameters of the stochastic model can be chosen to have additional modes including the cases of bistability and tristability.

Although we formulate our study in terms of steady state probability distributions, one may equally well view our results as describing the typical dynamic behavior of realizations of the stochastic process. These recapitulate the form of the steady state distributions: modes are reflected in metastable states along sample paths, states in which the system will stay for prolonged periods until switching to other states corresponding to alternative modes. In the SI, we provide Monte-Carlo simulations showing such metastable behavior along sample paths. We do so for the toggle switch as well as for a version of a well-studied genetic circuit [51] which exhibits oscillatory behavior along sample paths even though the corresponding deterministic model cannot admit oscillations.

The Reaction Network Structure

In this paper, a GRN will be formally defined as a set of nodes (genes) that are connected with each other through regulatory interactions via the proteins that the genes express. The regulatory proteins are called transcription factors (TFs). A TF regulates the expression of a gene by reversibly binding to the gene’s promoter and by either enhancing expression or repressing it.

The formalism we employ in order to describe GRNs at the elementary level is that of Chemical Reaction Networks (CRNs) [52]. A CRN consists of species and reactions, which we describe below.

Species:

The species in our context consist of promoter configurations for the various genes participating in the network, together with the respective TFs expressed from these genes and some of their multimers. A configuration of a promoter is characterized by the possible locations and number of TFs bound to the promoter at a given time. If a promoter is expressed constitutively, then there are two configurations specifying the expression activity state, active or inactive. A multimer is a compound consisting of a protein binding to itself several times. For instance, dimers and trimers are 2-mers and 3-mers, respectively. If a protein forms an $n^{\rm th}$-order multimer then we say that it has a cooperativity index of $n$. If species is denoted by $\rm X$, then its copy number is denoted by $X$.

For simplicity we assume the following:

1. A1

Each promoter can have up to two TFs binding to it. 2. A2

Each TF is a single protein that has a fixed cooperativity index, i.e, it cannot act as a TF with two different cooperativity indices. 3. A3

Each gene is present with only a single copy.

All the above assumptions can be relaxed. We make these assumptions only in order to simplify the notations and mathematical derivations. §SI-4,5 contain generalizations of the results to heterogenous TFs, and arbitrary numbers of gene copy numbers.

Consider the $i^{\rm th}$ promoter. The expression rate of a gene is dependent on the current configuration of its promoter. We call the set of all possible such configurations the binding-site set $B_{i}$. Each member of $B_{i}$ corresponds to a configuration that translates into a specific species $\mathrm{D}_{j}^{i},j\in B_{i}$. If a promoter has just one or no regulatory binding sites, then we let $B_{i}=\{0,1\}$. Hence, the promoter configuration can be represented by two species: the unbound species $\mathrm{D}_{0}^{i}$ and the bound species $\mathrm{D}_{1}^{i}$. If the promoter has no binding sites then the promotor configuration species are interpreted as the inactive and active configurations, respectively. On the other hand, if the promoter has two binding sites then $B_{i}=\{00,01,10,11\}$111We interpret the elements of the binding set as integers in binary representation.. The first digit in a member of $B_{i}$ specifies whether the first binding site is occupied, and the second digit specifies the occupancy of the second binding site. Hence, the promoter configuration can be represented by four species $\mathrm{D}_{00}^{i},\mathrm{D}_{10}^{i},\mathrm{D}_{01}^{i},\mathrm{D}_{11}^{i}$. Note that in general we need to define $2^{\kappa}$ species for a promoter with $\kappa$ binding sites.

The species that denotes the protein produced by the $i^{\rm th}$ gene is $\mathrm{X}_{i}$. A protein’s multimer is denoted by $\mathrm{X}_{ic}$. If protein $\mathrm{X}_{i}$ does not form a multimer then $\mathrm{X}_{ic}:=\mathrm{X}_{i}$.

Therefore, the set of species in the network is $\mathscr{S}=\bigcup_{i}\left(\{\mathrm{D}_{j}^{i},j\in B_{i}\}\cup\{\mathrm{X}_{i},\mathrm{X}_{ic}\}\right).$

Reactions:

In our context, the reactions consist of TFs binding and unbinding with promoters and the respective protein expression (with transcription and translation combined in one step), decay, and $n$-merization.

For each gene, we define a gene expression block. Each block consists of a set of gene reactions and a set of protein reactions as shown in Figure 2.

If promoter is constitutive, i.e. it switches between two configurations autonomously without an explicitly modeled TF-promoter binding, then $B_{i}=\{0,1\}$ and the gene reactions block consists of:

[TABLE]

We refer to $\mathrm{D}_{0}^{i}$ and $D_{1}^{i}$ as the inactive and active configurations, respectively. If the promoter has one binding site, then also $B_{i}=\{0,1\}$ and the gene reactions block consists of just two reactions:

[TABLE]

where $\mathrm{D}_{0}^{i}$ and $D_{1}^{i}$ are are the promoter configurations when unbound and bound to the TF, respectively. Note that we did not designate a specific species as the active one since it depends on whether the TF is an activator or a repressor. Specifically, when TF is an activator, $D_{1}^{i}$ will be the active configuration and $D_{0}^{i}$ will be the inactive configuration, and vice versa when TF is a repressor.

Finally, if the promoter has two TFs binding to it, then they can bind independently, competitively, or cooperatively. Cooperative binding is discussed in SI-§4.3.1. If they bind independently, then the promoter has two binding sites. Hence, $B_{i}=\{00,01,10,11\}$ and the gene block contains the following reactions:

[TABLE]

The activity of each configuration species is dependent on whether the TFs are activators or repressors, and on how they behave jointly. This can be characterized fully by assigning a production rate for each configuration as will be explained below.

In the case of competitive binding, two different TFs compete to bind to the same location. This can be modeled similarly to the previous case except that the transitions to $\mathrm{D}_{11}^{i}$, i.e. the configuration where both TFs are bound, are not allowed. Hence, the gene reactions block will have only two reactions: (3), (5), and the binding set reduces to $B_{i}=\{00,01,10\}$.

We assume that RNA polymerase and ribosomes are available in high copy numbers, and that we can lump transcription and translation into one simplified “production” reaction. The rate of production is dependent on the promoter’s configuration. So for each configuration $\mathrm{D}_{j}^{i},j\in B_{i}$ the production reaction is:

[TABLE]

where the kinetic constant $k_{ij}$ is a non-negative number. The case $k_{ij}=0$ means that when the promoter configuration is $\mathrm{D}_{j}^{i}$ there is no protein production, and hence $D_{j}^{i}$ is an inactive configuration. The promoter configuration can be ranked from the most active to the least active by ranking the corresponding production kinetic rate constants.

Consequently, the character of a TF is manifested as follows: if the maximal protein production occurs at a configuration with the TF being bound we say that the TF is activating, and if the reverse holds it is repressing. And, if the production is maximal with multiple configurations such that the TF is bound in some of them and unbound in others then the TF is neither repressing nor activating.

We model decay and/or dilution as a single reaction:

[TABLE]

The expressed proteins can act as TFs. They may combine to form dimers or higher order multimers before acting as TFs. The numbers of copies of the TF needed to form a multi-mer is called the the cooperativity index and we denote it by $n$. Hence, we model the cooperativity reactions as given in Figure 2 as follows:

[TABLE]

If the cooperativity index of $\mathrm{X}_{i}$ is 1, then the species $\mathrm{X}_{ic}:=\mathrm{X}_{i}$, and the multimerization reaction becomes empty.

Higher order multi-merization processes can be modelled as multi-step or sequential reactions [53]. We discuss how our theory includes this case in §SI-4.3.3, by showing how an equivalent one-step model with can be formulated.

Kinetics:

In order to keep track of molecule counts, each species $\mathrm{Z}_{i}\in\mathscr{S}$ is associated with a copy number $z_{i}\in\mathbb{Z}_{\geq 0}$.

To each reaction $\mathrm{R}_{j}$ one associates a propensity function $R_{j}$. We use Mass-Action Kinetics. If $\mathrm{R}_{j}$ has a single reactant species $Z_{i}$ with stoichiometry coefficient $\alpha_{ij}$, then [54]:

[TABLE]

where $k_{j}$ is a kinetic rate constant. Note if $\alpha_{ij}=1$, then $R_{j}(z_{i})=k_{j}z_{i}$.

The only bimolecular reaction we need is binding of a TF to a promoter, which has unity stoichiometry coefficients for each reactant species, i.e. the left side of the reaction is of the form $\mathrm{Z}_{i_{1}}+\mathrm{Z}_{i_{2}}$. In this case, the propensity function is:

[TABLE]

A gene regulatory network:

Consider a set of $N$ genes, binding sets $\{B_{i}\}_{i=1}^{N}$, and kinetic constants $k_{j}$’s. A gene expression block, as shown in Figure 2, is a set of gene reactions and protein reactions as defined above. Each gene block has an output that is either the protein or its $n$-mer, and it is designated by $\mathrm{X}_{ic}$. The input to each gene expression block is a subset of the set of the outputs of all blocks. Then, a GRN is an arbitrary interconnection of a gene expression blocks (Figure 2). SI-§4 defines a more general class of network that we can study.

A directed graph can be associated with a GRN as follows. Each vertex corresponds to a gene expression block. There is a directed edge from vertex A to vertex B if the output of A is an input to B. In order to simplify the presentation, we assume the following:

1. A4

The graph of gene expression blocks is connected.

Note that if A4 is violated, our analysis can be applied to each connected component.

Time-Scale Separation:

As mentioned in the introduction, we assume that the gene reactions (1)-(6) are considerably slower than the protein reactions (7)-(9). In order to model this assumption, we write the kinetic rates of gene reactions in the form $\varepsilon k_{j}$, where $0<\varepsilon\ll 1$ and assume that all other kinetic rates (for protein production, decay and multi-merization) are $\varepsilon^{-1}$-times faster.

Events in biological cells usually take place at different time-scales [1], and hence singular perturbation techniques are widely used in deterministic settings in order to reduce models for analysis. On the other hand, model-order reduction by time-scale separation in stochastic processes has been mainly used in the literature for computational purposes, for example to accelerate the stochastic simulation algorithm [55, 56], or to compute finite-space-projection solutions to the CME [30]. In this work, we use a singular perturbation approach for the analytical purpose of characterizing the form of the stationary distribution in the regimes of slow gene-TF binding/unbinding.

In the case of a finite Markov chain, the CME is a finite-dimensional linear ODE, and reduction methods for linear systems can be used [39] and applied to Markov chains [57, 58]. For continuous-time Markov chains on a countable space, as needed when analyzing gene networks, there are difficult and open technical issues. Exponential stochastic stability [59] needs to be established for the stationary solution in order to guarantee the existence of the asymptotic expansion in $\varepsilon$ [60]. Although it has been shown for a class of networks [61], the general problem needs further research. In this paper, we will not delve into technical issues of stochastic stability ; we assume that these expansions exist and that the solutions converge to a unique equilibrium solution.

Dynamics and the Master Equation

The dynamics of the network refers to the manner in which the state evolves in time, where the state $Z(t)\in\mathbf{Z}\subset\mathbb{Z}_{\geq 0}^{|\mathscr{S}|}$ is the vector of copy numbers of the species of the network at time $t$. The standard stochastic model for a CRN is that of a continuous Markov chain. Let $\mathbf{Z}$ denote the state space. Consider a time $t$ and let the state be $Z(t)=z\in\mathbf{Z}$. The relevant background is reviewed in SI-§1.1.

Let $p_{z}(t)=\Pr[Z(t)=z|Z(0)=z_{0}]$ be the stationary distribution for any given initial condition $z_{0}$. Its time evolution is given by the Chemical Master Equation (CME).

Since our species are either gene species or protein species, we split the stochastic process $Z(t)$ into two subprocesses: the gene process $D(t)$ and the protein process $X(t)$, as explained below.

For each gene we define one process $D_{i}$ such that $D_{i}(t)\in B_{i}$. $D_{i}(t)=j$ if and only the promoter configuration is encoded by $j\in B_{i}$. Collecting these into a vector, define the gene process $D(t):=[D_{1}(t),...,D_{N}(t)]^{T}$ where $D(t)\in\prod_{i=1}^{N}B_{i}$. The $i^{\rm th}$ gene can be represented by $|B_{i}|$ states, so $L:={\prod_{i=1}^{N}|B_{i}|}$ is the total number of promoter configurations in the GRN. With abuse of notation, we write also $D(t)\in\{0,..,L-1\}$ in the sense of the bijection between $\{0,..,L-1\}$ and $\prod_{i=1}^{N}B_{i}$ defined by interpreting $D_{1}...D_{N}$ as a binary representation of an integer. Hence, $d\in\{0,..,L-1\}$ corresponds to $(d_{1},...,d_{N})\in B_{1}\times..\times B_{N}$ and we write $d=(d_{1},..,d_{N})$.

Since each gene expresses a corresponding protein, we define $X_{i1}(t)\in\mathbb{Z}_{\geq 0},i=1,..,N$ protein processes. If the multimerized version of the $i^{\rm th}$ protein participates in the network as an activator or repressor then we define $X_{ic}(t)$ as the corresponding multimerized protein process, and we denote $X_{i}(t):=[X_{i1}(t),X_{ic}(t)]^{T}$. If there is no multimerization reaction then we define $X_{i}(t):=X_{i1}(t)$. Since not all proteins are necessarily multimerized, the total number of protein processes is $N\leq M\leq 2N$. Hence, the protein process is $X(t)=[X_{1}^{T}(t),..,X_{N}^{T}(t)]^{T}\in{\mathbb{Z}}_{\geq 0}^{M}$ and the state space can be written as $\mathbf{Z}={\mathbb{Z}}_{\geq 0}^{M}\times\prod_{i=1}^{N}B_{i}$.

Results

Decomposition of the Master Equation

It is crucial to our analysis to represent the linear system of differential equations given by the CME as an interconnection of weakly coupled linear systems. To this end, we present the appropriate notation in this subsection.

Consider the joint probability distribution:

[TABLE]

which represents the probability at time $t$ that the protein process $X$ takes the value $x\in{\mathbb{Z}}_{+}^{M}$ and the gene process $D$ takes the value $d\in\{0,..,L-1\}$. Recall that $x$ is a vector of copy numbers for the protein processes while $d$ encodes the configuration of each promoter in the network. Then, we can define for each fixed $d$:

[TABLE]

representing the vector enumerating the probabilities (10) for all values of $x$ and for a fixed $d$, where $x_{0},x_{1},..$ is an indexing of ${\mathbb{Z}}_{\geq 0}^{M}$. Note that $p_{d}(t)$ can be thought of as an infinite vector with respect to the aforementioned indexing. Finally, let

[TABLE]

representing a concatenation of the vectors (11) for $d=0,..,L-1$. Note that $p(t)$ is a finite concatenation of infinite vectors.

The joint stationary distribution $\bar{\pi}$ is defined as the following limit, which we assume to exist and is independent of the initial distribution:

[TABLE]

Note that $\bar{\pi}$ depends on $\varepsilon$.

Consider a given GRN. The CME is defined over a countable state space $\mathbf{Z}$. Hence, the CME can be interpreted as an infinite system of differential equations with an infinite infinitesimal generator matrix $\Lambda$ which contains the reaction rates (see SI-§1.1).

Consider partitioning the probability distribution vector as in (12). Recall that reactions have been divided into two sets: slow gene reactions (1)-(6) and fast protein reactions (7)-(9). This allows us to write $\Lambda$ as a sum of a slow matrix $\varepsilon\hat{\Lambda}$ and a fast matrix $\tilde{\Lambda}$, which we call a fast-slow decomposition. Furthermore, $\tilde{\Lambda}$ can be written as a block diagonal matrix with $L$ diagonal blocks which correspond to conditioning the Markov chain on a specific gene state $d$. This is stated in the following basic proposition (see SI-§2.1 for the proof):

Proposition 1**.**

Given a GRN. Its CME can be written as

[TABLE]

where

[TABLE]

where $\tilde{\Lambda}$ is the fast matrix, $\hat{\Lambda}$ is the slow matrix, and $\Lambda_{0},..,\Lambda_{L-1}$ are stochastic matrices.

Conditional Markov Chains

For each $d$, consider modifying the Markov chain $Z(t)$ defined in the previous section by replacing the stochastic process $D(t)$ by a deterministic constant process $D(t)=d$. This means that the resulting chain does not describe the gene process dynamics, it only describes the protein process dynamics conditioned on $d$. Henceforth, we refer to the resulting Markov chain as the Markov chain conditioned on $d$. The infinitesimal generator of a chain conditioned on $d$ is denoted by $\Lambda_{d}$, and is identical to the corresponding block on the diagonal of $\tilde{\Lambda}$ as given in (15). In other words, fixing $D(t)=d\in\{0,..,L-1\}$, the dynamics of the network can be described by a CME:

[TABLE]

where $p_{X|d}$ is a vector that enumerates the conditional probabilities $p_{x|d}=\Pr[X(t)=x|D(t)=d]$ for a given $d$. The conditional stationary distribution is denoted by: $\pi_{X|d}^{(J)}=\lim_{t\to\infty}p_{X|d}(t)$, where $(J)$ refers to the fact that it is joint in the protein and multimerized protein processes. Note that $\pi_{X|d}^{(J)}$ is independent of $\varepsilon$. This notion of a conditional Markov chain is useful since, at the slow promoter kinetics limit, $D(t)$ stays constant. It can be noted from (15) that when $\varepsilon=0$ the dynamics of $p_{d}$ decouples and becomes independent of $p_{\tilde{d}},\tilde{d}=0,..,L-1,\tilde{d}\neq d$.

We show below that each conditional Markov chain has a simple structure. Fixing the promoter configuration $D(t)=d=(d_{1},..,d_{N})$, the network consists of uncoupled birth-death processes. So for each $d_{i}$, the protein reactions (7)-(9) corresponding to the $i^{\rm th}$ promoter can be written as follows without multimerization:

[TABLE]

where the subscript $id_{i}$ refers to the production kinetic constant corresponding to the configuration species $\mathrm{D}_{d_{i}}^{i}$, or, if there is a multimerization reaction, it takes the form:

[TABLE]

Note that the stochastic processes $\mathrm{X}_{i}(t),i=1,..,N$ conditioned on $D(t)=d$ are independent of each other. Hence, the conditional stationary distribution $\pi_{X|d}^{(J)}$ can be written as a product of stationary distributions and the individual stationary distributions have Poisson expressions. The following proposition gives the analytic expression of the conditional stationary distributions: (see SI-§2.2 for proof)

Proposition 2**.**

Fix $d\in\{0,..,L-1\}$. Consider (16), then there exists a conditional stationary distribution $\pi_{X|d}^{(J)}$ and it is given by

[TABLE]

where

[TABLE]

where $(J)$ refers to the joint distribution in multimerized and non-multimerized processes, $x_{i1}$ refers to the copy number of $\mathrm{X}_{i}$, while $x_{i2}$ refers to the copy number of $\mathrm{X}_{ic}$, $\mathbf{P}(x;a):=\frac{a^{x}}{x!}e^{-a},\mathbf{P}(x_{1},x_{2};a_{1},a_{2}):=\frac{a_{1}^{x_{1}}}{x_{1}!}\frac{a_{2}^{x_{2}}}{x_{2}!}e^{-a_{1}-a_{2}}$.

Remark 1**.**

The conditional distribution in (19) is a joint distribution in the protein and multimerized protein processes. If we want to compute a marginal stationary distribution for the protein process only, then we average over the multimerized protein processes $X_{ic},i=1,..,N$ to get a joint Poisson in $N$ variables. Hence, the formulae (19) -(20) can be replaced by:

[TABLE]

where $M-N$ is the number of $n$-merized protein processes, and $\pi_{X|d}$ is the marginal stationary distribution for the protein process.

Decomposition of The Stationary Distribution

Recall the slow-fast decomposition of the CME in (14) and the joint stationary distribution (13). In order to emphasize the dependence on $\varepsilon$ we denote $\bar{\pi}^{\varepsilon}:=\bar{\pi}(\varepsilon)$. Hence, $\bar{\pi}^{\varepsilon}$ is the unique stationary distribution that satisfies $\Lambda_{\varepsilon}\bar{\pi}^{\varepsilon}=0$, $\pi^{\varepsilon}>0$, and $\sum_{z}\pi_{z}^{\varepsilon}=1$, where the subscript denotes the value of the stationary distribution at $z$.

Our objective is to characterize the stationary distribution as $\varepsilon\to 0$. Writing $\bar{\pi}_{\varepsilon}$ as an asymptotic expansion to first order in terms of $\varepsilon$, we have

[TABLE]

Our aim is to find $\bar{\pi}^{(0)}$. We use singular perturbations techniques to derive the following theorem (see SI-§2.3):

Theorem 3**.**

Consider a given GRN with $L$ promoter states with the CME (14). Writing (22), then the joint stationary distribution $\bar{\pi}:=\lim_{\varepsilon\to 0^{+}}\bar{\pi}^{\varepsilon}$ can be written as:

[TABLE]

where $\lambda=[\lambda_{0},..,\lambda_{L-1}]^{T}$ is the principal normalized eigenvector of:

[TABLE]

where $\bar{\pi}_{X|0},...,\bar{\pi}_{X|L-1}$ are the extended conditional stationary distributions defined as: $\bar{\pi}_{X|d}(x,d)=\pi_{X|d}(x)$, and $\bar{\pi}_{X|d}(x,d^{\prime})=0$ when $d^{\prime}\neq d$.

The result characterizes the stationary solution of (14) which is a joint distribution in $X$ and $D$. However, we are particularly interested in the marginal stationary distribution of the protein process $X$ and the marginal stationary distribution of the non-multimerized protein process, since these distributions are typically experimentally observable. Therefore, we can use Remark 1 to write the stationary distribution as mixture of $L$ Poisson distributions with weights $\{\lambda_{d}\}_{d=0}^{L-1}$:

Corollary 4**.**

Consider a given GRN with $L$ genes with the CME (14). Writing (22), let $\pi_{X|0},...,\pi_{X|L-1}$ be the conditional stationary distributions of $\Lambda_{0},...,\Lambda_{L-1}$, where explicit expressions are given in (19). Then, we can write the following:

[TABLE]

*where $\lambda=[\lambda_{0},..,\lambda_{L-1}]^{T}$ is as given Theorem 3.

Furthermore, the marginal stationary distribution of the non-multimerized protein process can be written as:*

[TABLE]

Remark 2**.**

In the remainder of the Results section, when we refer to the “stationary distribution” we mean the marginal stationary distribution of the non-multimerized protein process given in (25).

Remark 3**.**

If a mode is defined as a local maximum of a stationary distribution, then this does not necessarily imply that the stationary distribution has $L$ modes since the peak values of two Poisson distributions can be very close to each other. In the remainder of the paper we will call each Poisson distribution in the mixture as a “mode“ in the sense that it represents a component in the mixture distribution. The number of local maxima of a distribution can be found easily given the expression (25).

The Reduced-Order Finite Markov Chain

The computation of the weighting vector $\lambda$ in Theorem 3 requires computing the $L\times L$ matrix $\Lambda_{r}$ in (23) which can be interpreted as the infinitesimal generator of an $L$-dimensional Markov chain. The expression in (23) involves evaluating the product of infinite dimensional matrices. Since the structure of the GRN and the form of the conditional distribution in (19) are known, an easier algorithm to compute $\Lambda_{r}$ for our GRNs is given in Proposition SI-2. The algorithm provides an intuitive way to interpret Theorem 3 and can be informally described as follows.

Assume $D(t)=d$, the algorithm implies that each binding reaction of the form:

[TABLE]

gives the rate $\alpha\mathbb{E}[{TF}|D=d]$, where $\mathbb{E}$ denotes mathematical expectation. Hence it corresponds to a reaction of the following form in the reduced-order Markov chain:

[TABLE]

Using Proposition 2, we can write:

[TABLE]

Refer to §SI-2.4 for a more precise statement.

Basic Example: Gene Bursting Model

The simplest network is the autonomous TF-gene binding/unbinding model, and it has been used for transcriptional bursting [17] and studied using time-scale separation in [32, 62]. Consider:

[TABLE]

Referring to Figure 2, we identify a single gene block with two states. Using (19), the conditional stationary distributions are two Poissions at [math] and $k/k_{-}$. The reduced Markov chain is a binary Bernoulli process with a rate of $\alpha/(\alpha+\alpha_{-})$. Then the stationary distribution of $X$ can be written using (25) as: (see SI-§3.1)

[TABLE]

which is a bimodal distribution with peaks at 0 and $k/k_{-}$. The fast promoter kinetics model is obtained, instead, by reversing the time-scale separation such that the protein reactions become slow and gene reactions become fast. In that case, the resulting stationary distribution is a Poisson with mean $\frac{\alpha}{\alpha+\alpha_{-}}\frac{k}{k_{-}}$ which is the same as the deterministic equilibrium (with the conservation law $D_{1}(t)+\mathrm{D}_{0}(t)=1$). Although two models share the mean, the stationary distributions differ drastically.

The Role of Cooperativity

A TF is said to be cooperative if it acts only after it forms a dimer or a higher-order $n$-mer that binds to the gene’s promoter [53]. In standard deterministic modelling, a cooperative activation changes the form of the quasi-steady state activation rate from a Michaelis-Menten function into a Hill function. Cooperativity is often necessary for a network to have multiple equilibria in some kinetic parameter ranges. For example, a non-cooperative self- activating gene can only be mono-stable, while its cooperative counterpart can be multi-stable for some parameters.

Corollary 4 and (27) show that cooperativity plays in the context of slow promoter kinetics a role that is very different from the deterministic setting. This is since the stationary distribution is a mixture of $L$ Poisson processes, independent of whether the activations are multimerized or not, and (21) are also independent of the dimerization rates. Nevertheless, the multi-merization can tune the weighting coefficients in (27). In the non-cooperative case, a certain mode can be made more probable only by changing either the location of the mode or the dissociation ratio (the ratio of the binding to unbinding kinetic constants). On the other hand, a multimerized TF gives extra tuning parameters, namely the multimerization ratio and the cooperativity index. Hence, a certain mode can be made more or less probable by modifying the multimerization ratio and/or the cooperativity index without changing the location of the peaks or the dissociation ratio.

In order to illustrate the above idea, we analyze a self-regulating gene with slow promoter kinetics with and without cooperativity.

A Self-Regulating Gene

Consider a non-cooperative self-regulating gene:

[TABLE]

The network is activating if $k_{1}>k_{0}$, and repressing otherwise.

Referring to Figure 2, this is a single-gene block with two states. The reduced Markov chain is a binary Bernoulli process with the rate $\alpha k_{0}/(\alpha k_{-}+\alpha k_{0})$. Using (25) the stationary distribution is:

[TABLE]

where

[TABLE]

Next, consider the same reaction network, but now with cooperativity:

[TABLE]

In this case, the stationary distribution is still given by (30) but the weighting parameter changes to

[TABLE]

In both cases the distribution has modes at $\frac{k_{1}}{k_{-}}$ and $\frac{k_{0}}{k_{-}}$, where the height of the first mode is proportional to $\rho$.

Comparing (31) and (33), note that, in the non-cooperative case, if we want to increase the weight of the mode corresponding to the bound state keeping the dissociation ratio, then the mode location needs to be changed. On the other hand, the dimerization rates in (33) can be used in order to tune the weights freely while keeping the modes and the binding to unbinding kinetic constants ratio unchanged. For instance, the distribution can be made effectively unimodal with a sufficiently high dimerization ratio.

Comparison with the deterministic model:

Table 1 compares the number of stable equilibria in the deterministic model with the number of modes in the stochastic model in the case of a single gene copy. It can be noted that there is no apparent correlation between the numbers of deterministic equilibria and stochastic modes.

Figure 3 depicts the transition from a unique mode with fast promoter kinetics to multiple modes with slow kinetics with cooperativity and leakiness for a self-activating gene.

The Toggle Switch

A toggle switch is a basic GRN that exhibits deterministic multi-stability. It has two stable steady states and can switch between them with an external input or via noise. The basic design is a pair of two mutually repressing genes as in Figure 4-a. The ideal behavior is that only one gene is “on” at any moment in time. The network can be given by the following network with cooperativity indices $n,m$:

[TABLE]

$\mathrm{D}^{X}$, $\mathrm{D}^{Y}$ denote the states of the promoters of the two genes expressing $X,Y$, respectively.

For the case $n,m=1$, there is no multi-merization reaction. For consistency, we choose $\beta_{1}=\beta_{-1}=\beta_{2}=\beta_{-2}=1$ in that case.

Using the algorithm of Proposition SI-2 we get that the distribution has three modes only (see SI-§3.3). The stationary distribution for $X,Y$ is:

[TABLE]

where

[TABLE]

Since the stationary distribution has three modes, it deviates from the ideal behavior of a switch where at most two stable steady states, under appropriate parameter conditions, are possible. Nevertheless, a bimodal distribution can be achieved by minimizing the weight of the first mode at $(\frac{k_{10}}{k_{-1}},\frac{k_{20}}{k_{-2}})$. If we fix $\alpha_{1}/\alpha_{-1},\alpha_{2}/\alpha_{-2}$, then this can be satisfied by tuning $n,m,\beta_{\pm 1},\beta_{\pm 2}$ to maximize $\rho_{1},\rho_{2}$ in (35). Choosing higher cooperativity indices, subject to $n<k_{10}/k_{-1},m<k_{20}/k_{-2}$, achieves this. For instance, a standard design [63] uses $n=2,m=3$. Figure 4 depicts the effect of cooperativity on achieving the desired behavior with the same dissociation constant and production ratios, and dimerization ratios equal to one. Notice that cooperativity allow us to minimize or maximize the weight of the mode corresponding to both proteins at high concentrations.

The toggle switch has three modes regardless of the cooperativity index. This is unlike the deterministic model where only one positive stable state is realizable with non-cooperative binding, and two stable steady states are realizable with cooperative binding. §SI-3.3 contains further Monte-Carlo simulations that show that the predicted third mode appears with a 0.5 time scale separation. Experimentally, a recent pre-print has reported that the CRI-Cro toggle switch exhibits the third (high,high) mode and the authors proposed slow-promoter kinetics as a contributing mechanism [64].

Synchronization of interconnected toggle switches

We consider $N$ identical toggle switches:

[TABLE]

where $i=1,..,N$. We interconnect these systems through diffusion of the protein species $\mathrm{X}_{i},\mathrm{Y}_{i}$ among cells, modeled through reversible reactions with a diffusion coefficient $\Omega$:

[TABLE]

We study this model as a very simplified version of a more complex quorum sensing communication mechanism, in which orthogonal AHL molecules are produced and by cells and act as activators of TFs in receiving cells, as analyzed for example in [44].

Figure 5(a)a depicts a block diagram of such a network.

For a deterministic model, there exists a parameter range for which all toggle switches will synchronize into bistability for sufficiently high diffusion coefficient [44]. This implies each switch in the network behaves as a bistable switch, and it converges with all the other switches to the same steady-states.

Our aim is to analyze the stochastic model at the limit of slow promoter kinetics and compare it with the deterministic model.

This network is not in the form of the class of networks in Figure 1. Nevertheless, we show in SI-§4.1 that our results can be generalized to networks that admit weakly reversible deficiency zero conditional Markov chains.

There are $4^{N}$ conditional Markov chains, and using Theorem 3 and Proposition 5, the stationary distribution is a mixture of $4^{N}-1$ Poissons.

Consider now the case of a high diffusion coefficient. We show (see SI-§3.5) that as $\Omega\to\infty$, $X_{1},..,X_{N}$ will synchronize in the sense that the joint distribution of $X_{1},..,X_{N}$ is symmetric with respect to all permutations of the random variables. This implies that the marginal stationary distributions $p_{X_{i}},i=1,..,N$ are identical. Hence, for sufficiently large $\Omega$, the probability mass is concentrated around the region for which $X_{1},..,X_{N}$ are close to each other. Consequently, for large $\Omega$ we can replace the population of toggle switches with a single toggle switch with the synchronized protein processes $X(t),Y(t)$, which are defined, for the sake of convenience, as $X(t):=X_{1}(t),Y(t):=Y_{1}(t)$. Next, we describe the stationary distribution of $X(t),Y(t)$.

The state of synchronized toggle switches does not depend on individual promoter configurations, and it depends only on the total number of unbound promoter sites in the network. Hence, the number of modes will drop from $4^{N}-1$ to $(N+1)^{2}-1$. Note that similar to the single toggle switch, there are modes which have both $X,Y$ with non-zero copy number. On the other hand, there are many additional modes. Recall that in the case of a single toggle switch, we have tuned the cooperativity ratios such that the modes in which both genes are ON are suppressed. Similarly, the undesired modes can be suppressed by tuning the cooperativity ratio which can be achieved by choosing $\rho_{d_{i}}^{X},\rho_{d_{i}}^{Y},d=0,..,4^{N}-1$ sufficiently large. In particular, letting the multi-merization ratio $\beta_{x}/\beta_{-x},\beta_{y}/\beta_{-y}\to\infty$, the weights of modes in the interior of the positive orthant $\mathbb{R}_{+}^{2}$ approach zero.

In conclusion, for sufficiently high $\Omega$ and sufficiently high multimerization ratio the population behaves as a *multimodal switch *, which means that the whole network can have either the gene $\mathrm{X}$ ON, or the gene $\mathrm{Y}$ ON. And every gene can take $2N$ modes which are:

[TABLE]

Comparing to the low diffusion case, the network will have up to $2^{N}-1$ modes with sufficiently high multimerization ratio.

In order to illustrate the previous results, consider a population of three toggle switches ($N=3$) and cooperativity $n=2$. For $\Omega$ greater then a certain threshold, the deterministic system bifurcates into bistabiliy. This means that all toggle switches converge to the same exact equilibria if $\Omega$ is greater than the threshold. In contrast, the modes in the stochastic model of the toggle switches converge asymptotically to each other. Hence, we need to choose a threshold for $\Omega$ that constitutes “sufficient” synchronization. We choose to define this as the protein processes synchronizing within one copy number. In other words, we require the maximum distance between the modes to be less than 1. It can be shown (see SI-§3.5) that the diffusion coefficient needs to satisfy:

[TABLE]

The minimal $\Omega$ that satisfies the inequality is $\Omega=75$ in this example. The stationary distribution is depicted in Figure 5(a)d. The network has 15 modes, nine of which are in the interior are suppressed due to cooperativity. Comparing with the deterministic model, it bifurcates into synchronization for $\Omega>0.5$. The stable equilibria of synchronized switch are $(149.98,0.022),(0.022,149.98)$.

The stochastic model with slow promoter kinetics adds four additional modes at $(0,100)$, $(100,0)$, $(50,0)$, $(0,50)$. This can be interpreted in the following manner. In the stochastic model, the protein processes synchronize while the promoter configurations do not. The high states $(150,0),(0,150)$ correspond to the case when all the binding sites are empty. In the case when one binding site is empty, the first gene is producing while the second and the third are not. Due to diffusion, the first gene “shares” its expressed protein with the other two genes, which implies that each gene will receive a third of the total protein copy numbers produced in the network. A similar situation arises when two binding sites are empty.

Trans-Differentiation Network

We consider two networks for TF cross-antagonism in cell fate decision in this section. Both networks consist of two self-activating genes repressing each other as depicted in Figure 6-a [5]. The first network has independent cooperative binding of the TFs to the promoters. So it can be written as follows [37]:

[TABLE]

In order for the genes to be cross-inhibiting and self-activating we let: $k_{x1}=k_{y2}=0$. Also, $k_{x2}>k_{x0},k_{x3}$ and $k_{y1}>k_{y0},k_{y3}$.

The network can be analyzed with the proposed framework, as it consists of two genes each with two binding sites. Hence it can theoretically admit up to 16 modes according to (25). The stationary distribution is depicted in Figure 6-b for an example parameter set. Note that despite the fact that we have 16 modes, only eight of them contribute to most of the stationary distribution. This is to be contrasted with a deterministic model, which cannot produce more than 4 stable equilibria [50].

The second network that we study is a model of the PU.1/GATA.1 network, which is a lineage determinant in hematopoietic stem cells [65]. Diagrammatically, it can also be presented by Figure 6-a. However, it differs from the first network presented above in several ways. First, PU.1 needs GATA.1 to bind to the promoter of GATA.1 [66], and vice versa [67]. In our modelling framework this means that the promoter configurations $\mathrm{D}_{01}^{X},D_{10}^{Y}$ do not exist, where $\mathrm{X}$ stands for PU.1 and $Y$ stands for GATA.1. Hence, the network has nine gene states. Second, there is no evidence that PU.1 and GATA.1 form dimers to activate their own promoters cooperatively. In fact, it has been shown that self-activation for GATA-1 occurs primarily through monomeric binding [68]. Therefore, the PU.1/GATA.1 network can be written as follows:

[TABLE]

More detailed discussion of the model is included in [69].

With lack of cooperativity, a deterministic model is only monostable and cannot explain the emergence of bistability for the above network [69]. However, using our framework, up to nine modes can be realized. In order to simplify the landscape, we group the nine into four modes. This is possible since the states $\mathrm{D}_{11}^{X},D_{00}^{X}$, $D_{11}^{Y},D_{00}^{Y}$ have very low production rates. This gives a total of four modes which are (low,low),(high,low),(low,high),(high,high). Using our model, we choose the parameters to realize bistability and tristability. Figure 6-c depicts the stationary distribution for a set of parameters that satisfies the assumptions and give rise to a tristable distribution.

The Repressilator

A very different example is provided by a well-studied synthetic oscillator, the repressilator [51]. The slow gene activation model predicts oscillations even in the noncooperative binding regime in which the deterministic model as well as the adiabatic stochastic model do not oscillate. See §SI-3.4.

Methods

Numerical Simulation Software

All calculations were performed using MATLAB R2016a, except for the computation of deterministic solutions of the “quorum sensing” numerical example where we used Bertini 1.5, which is software for solving polynomials numerically via homotopy methods.

Discussion

Phenotypical variability in the absence of genetic variation is a phenomenon of great interest in current biological and translational research, as it plays an important role in processes as diverse as embryonic development [70], hematopoietic cell differentiation [71], and cancer heterogeneity [72]. A conceptual, and often proposed, unifying framework to explain non-genetic variability is to think of distinct phenotypes as multiple “metastable states” or “modes” in the complex energetic landscape associated to an underlying GRN. Following this point of view, we studied in this paper a general but simplified mathematical model of gene regulation. Our focus was on stochastic slow promoter kinetics, the time scale relevant when transcription factor binding and unbinding are affected by epigenetic processes such as DNA methylation and chromatin remodeling. In that regime, adiabatic approximations of promoter kinetics are not appropriate. In contrast to the existing literature, which largely confines itself to numerical simulations, in this work we provided a rigorous analytic characterization of multiple modes.

The general formal approach that we developed provides insight into the relative influence of model parameters on system behavior. It also allows making theoretical predictions of how changes in wiring of a gene regulatory network, be it through natural mutations or through artificial interventions, impact the possible number, location, and likelihood, of alternative states. We were able to tease out the role of cooperative binding in stochastic models in comparison to deterministic models, which is a question of great interest in both the analysis of natural systems and in synthetic biology engineering. Specifically, we found that, unlike deterministic systems, the number of modes is independent of whether the TF-promoter binding is cooperative or not; on the other hand, cooperative binding gives extra degrees of freedom for assigning weights to the different modes. More generally, we characterized the stationary distributions of CMEs for our GRNs as mixtures of Poisson distributions, which enabled us to obtain explicit formulas for the locations and probabilities of metastable states as a function of the parameters describing the system. One application of our mathematical results was to models of single or communicating “toggle switches” in bacteria, where we showed that, for suitable parameters, there are a very large number of metastable attractors. Indeed, stochastic effects have been shown before to lead to multi-stability in a population of enzymatic reactions [73].

This work was in fact motivated by our interest in hematopoietic cell differentiation, and in this paper we discussed two possible models of trans-differentiation networks in mammalian cells. In a first model, based on previous publications, we uncovered more modes than had been predicted with different analyses of the same model. This implies that in practice there could be unknown “intermediate” phenotypes that result from the network’s dynamics, which may be acquired by cells during the natural differentiation process or which one might be able to induce through artificial stimulation. The second model included only binding reactions that have been experimentally documented, and as such might be more biologically realistic than the first model. For this second model, a deterministic analysis predicts monostability, which is inconsistent with the fact that the network should control a switch between two stable phenotypes (erythroid and myeloid). This suggests that stochasticity, likely due to low copy numbers and/or slow promoter kinetics, might be responsible for the multiple attractors (phenotypes) that are possible in cell differentiation gene regulatory networks.

Our mathematical results, being quite generic, should also be useful in the analysis of networks that have been proposed for understanding aspects of cancer biology. For example, non-genetic heterogeneity has been recently recognized as an important factor in cancer development and resistance to therapy, with stochastic multistability in gene expression dynamics acting as a generator of phenotype heterogeneity, setting a balance between mesenchymal, epithelial, and cancer stem-cell-like states [74] [75] [76] [77], and nongenetic variability due to multistability arising from mutually repressing gene networks has been proposed to explain metastatic progression [78].

Acknowledgements

We thank Nithin S. Kumar for discussions regarding the PU.1/GATA.1 network and Cameron McBride for proofreading the manuscript. This work was supported by an AFOSR grant FA9550-14-1-0060.

Supporting Information

1 Review of the Master Equation and Markov Chains

1.1 The Chemical Master Equation

A generic reaction $\mathrm{R}_{j}\in\mathscr{R}$ takes the form:

[TABLE]

where $Z_{i}\in\mathscr{S}$, $\alpha_{ij},$ and $\beta_{ij}$ are positive integers. The reactions that we consider are limited to at most two reactants. The reverse reaction of $\mathrm{R}_{j}$ is the reaction in which the products and reactants are interchanged. If the network contains both the reaction and its reverse then we use the short-hand notation to denote both of them as

[TABLE]

The stoichiometry of a CRN can be summarized by a stoichiometry matrix $\Gamma$ which is defined element-wise as follows:

[TABLE]

The columns of the stoichiometry matrix $\gamma_{1},..,\gamma_{|\mathscr{R}|}$ are known as the stoichiometry vectors. We say that a nonzero nonnegative vector $d$ gives a conservation law for the stoichiometry if $d^{T}\Gamma=0$.

The dynamics of the network refers to the manner in which the state evolves in time, where the state $Z(t)\in\mathbf{Z}\subset\mathbb{Z}_{\geq 0}^{|\mathscr{S}|}$ is the vector of copy numbers of the species of the network at time $t$. Since the collision of molecules is random in nature, the time-evolution of states is described mathematically by a stochastic process. The standard stochastic model for a CRN is that of a continuous Markov chain. Let $\mathbf{Z}$ denote the state space. Consider a time $t$ and let the state be $Z(t)=z\in\mathbf{Z}$. Then, the probability that the $j^{\rm th}$ reaction fires in an interval $[t,t+\delta]$ is $R_{j}(z)\delta+o(\delta)$. If $\mathrm{R}_{j}$ fires, then the states changes from $z$ to $z+\gamma_{j}$, where $\gamma_{j}$ is the corresponding stoichiometric vector.

As $Z$ is a stochastic process we are interested in characterizing its qualitative behavior given by the joint probability distribution $p_{z}(t)=\Pr[Z(t)=z|Z(0)=z_{0}]$ for any given initial condition $z_{0}$. The time-evolution of the probability distribution can be shown [54] to be given by a system of linear ordinary differential equations known as the forward Kolmogorov equation or the Chemical Master Equation, given by:

[TABLE]

where $\gamma_{1},..,\gamma_{|\mathscr{R}|}$ are the columns of the stoichiometry matrix.

Since our species are either gene species or protein species, we split the stochastic process $Z(t)$ into two subprocesses: the gene process $D(t)$ and the protein process $X(t)$, as explained below.

Consider the $i^{\rm th}$ gene. For each configuration species $\mathrm{D}_{j}^{i}\in B_{i}$, let $D_{j}^{i}(t)\in\{0,1\}$ denote its occupancy, i.e. if $D_{j}^{i}(t)=1$, then at time $t$ the $i^{\rm th}$ gene is in a configuration $j\in B_{i}$. It can be seen from gene reactions (1)-(6) that the network always has a conservation law supported on $\{\mathrm{D}_{j}^{i},j\in B_{i}\}$, so that:

[TABLE]

which reflects the physical constraint that the promoter can be in only one configuration at any given time.

This conservation law enables us to introduce an equivalent reduced representation. For each gene we define one process $D_{i}$ such that $D_{i}(t)\in B_{i}$. $D_{i}(t)=j$ if and only if $D_{j}^{i}(t)=1$. Collecting these into a vector, define the gene process $D(t):=[D_{1}(t),...,D_{N}(t)]^{T}$ where $D(t)\in\prod_{i=1}^{N}B_{i}$. The $i^{\rm th}$ gene can be represented by $|B_{i}|$ states, so $L:={\prod_{i=1}^{N}|B_{i}|}$ is the total number of promoter configurations in the GRN. With abuse of notation, we write also $D(t)\in\{0,..,L-1\}$ in the sense of the bijection between $\{0,..,L-1\}$ and $\prod_{i=1}^{N}B_{i}$ defined by interpreting $D_{1}...D_{N}$ as a binary representation of an integer. Hence, $d\in\{0,..,L-1\}$ corresponds to $(d_{1},...,d_{N})\in B_{1}\times..\times B_{N}$ and we write $d=(d_{1},..,d_{N})$.

Since each gene expresses a corresponding protein, we define $X_{i1}(t)\in\mathbb{Z}_{\geq 0},i=1,..,N$ protein processes. If the multimerized version of the $i^{\rm th}$ protein participates in the network as an activator or repressor then we define $X_{ic}(t)$ as the corresponding multimerized protein process, and we denote $X_{i}(t):=[X_{i1}(t),X_{ic}(t)]^{T}$. If there is no multimerization reaction then we define $X_{i}(t):=X_{i1}(t)$. Since not all proteins are necessarily multimerized, the total number of protein processes is $N\leq M\leq 2N$. Hence, the protein process is $X(t)=[X_{1}^{T}(t),..,X_{N}^{T}(t)]^{T}\in{\mathbb{Z}}_{\geq 0}^{M}$ and the state space can be written as $\mathbf{Z}={\mathbb{Z}}_{\geq 0}^{M}\times\prod_{i=1}^{N}B_{i}$.

Consider the joint probability distribution:

[TABLE]

which represents the probability at time $t$ that the protein process $X$ takes the value $x\in{\mathbb{Z}}_{+}^{M}$ and the gene process $D$ takes the value $d\in\{0,..,L-1\}$. Recall that $x$ is a vector of copy numbers for the protein processes while $d$ encodes the configuration of each promoter in the network. Then, we can define for each fixed $d$:

[TABLE]

representing the vector enumerating the probabilities (43) for all values of $x$ and for a fixed $d$, where $x_{0},x_{1},..$ is an indexing of ${\mathbb{Z}}_{\geq 0}^{M}$. Note that $p_{d}(t)$ can be thought of as an infinite vector with respect to the aforementioned indexing. Finally, let

[TABLE]

representing a concatenation of the vectors (44) for $d=0,..,L-1$. Note that $p(t)$ is a finite concatenation of infinite vectors.

The joint stationary distribution $\bar{\pi}$ is defined as the following limit, which we assume to exist and be independent of the initial distribution:

[TABLE]

The stationary distribution $\bar{\pi}$ is a function of $\varepsilon$ also.

Consider a given GRN. The master equation (42) is defined over a countable state space $\mathbf{Z}$ which can be enumerated with an arbitrarily chosen order. Hence, the master equation can be interpreted as an infinite system of differential equations. Its infinite infinitesimal generator matrix $\Lambda$ can be written succinctly entry-wise as:

[TABLE]

where $\lambda_{z\tilde{z}}$ refers to the rate of transition from $z$ to $\tilde{z}$. The matrix $\Lambda$ is stochastic, which means that it is Metzler and $1^{T}\Lambda=0$. A Metzler matrix is a matrix whose off-diagonal elements are non-negative.

1.2 Irreducibility

An important property in the context of Markov chain analysis is that of irreducibility [11], and its significance stems from the fact that it is a necessary condition for the existence of a unique positive stationary distribution. Consider the Markov chain $Z(t)$ defined on $\mathbf{Z}$ with an associated infinitesimal generator $\Lambda$ as given in (47). Let $z,w\in\mathbf{Z}$. Then, it is said that $z$ leads to $w$ if there exist states $z_{0},...,z_{n}\in\mathbf{Z}$ such that $\lambda_{zz_{0}}\lambda_{z_{0}z_{1}}...\lambda_{z_{n}w}>0$. A set $U\subset\mathbf{Z}$ is said to be a communicating class if for every $z_{1},z_{2}\in U$, $z_{1}$ leads to $z_{2}$ and $z_{2}$ leads to $z_{1}$. The state space $\mathbf{Z}$ can always be partitioned into a disjoint union of communicating classes [11]. The Markov chain is said to be irreducible if the state space is a communicating class. A communicating class $U$ is said to be closed if $z\in U$, and $z$ leads to $w$ implies $w\in U$. A Markov chain is said to be weakly irreducible if it has a unique closed communicating class $U$, and for all $z\in\mathbf{Z}$, $z$ leads to some element $U$.

We state the following result, under assumption A4:

Proposition 5**.**

Consider a gene regulatory network that consists of $N$ gene expression blocks. Then the associated Markov chain is weakly irreducible.

Proof.

Consider the state $0\in\mathbf{Z}$. We first show that for all $z\in\mathbf{Z}$, $z$ leads to [math]. Let $z=(x_{1},d_{1},..,x_{N},d_{N})$. We list the set of reactions, i.e transitions, that will lead to 0. Consider $d_{i}\neq 0$, if $d_{i}=1$ then we apply either the reaction (1) or (2). If $d_{i}=10$, then we apply reaction (3) and if $d_{i}=01$ we apply reaction (5). If $d_{i}=11$, then we apply reactions (3) and (4). Hence, $z$ leads to a state of the form $(x_{1},0,x_{2},0,...,x_{N},0)$. Similarly, we can apply the decay reactions (8) and the reverse dimerization until we reach the origin.

Now we show that there exists a closed communicating class. If [math] does not lead to any state then $\{0\}$ is a closed communicating class. Otherwise, let $U$ be the smallest communicating class containing [math]. Note that $U$ is closed, since if there exists $z\in U$ that leads to $w$, then $w$ leads to $0\in U$.

In order to show that $U$ is unique, assume that there exists another closed communicating class $U^{\prime}$. But this contradicts with the fact that all $z\in U^{\prime}$ lead to $0\in U$. We have shown that for all $z\in\mathbf{Z}$, $z$ leads to 0. Hence $z$ leads to $U$. ∎

Remark 4**.**

For finite Markov chains, weak irreducibility with appropriate stochastic stability assumptions are sufficient for the existence of a nonnegative unique stationary distribution [58], while irreducibility is usually needed for the existence a positive stationary distribution. Note that not all GRNs are irreducible. However, our subsequent results require weak irreducibility only, and investigation of irreducibility is out of the scope of this paper. Nevertheless, necessary and sufficient graphical conditions for irreducibility can be developed and are subject to future work.

2 Proofs of the Main Results

In this section we include mathematical proofs of the main results in the main text.

2.1 Decomposition of the Master Equation

We include a proof for Proposition 1.

By the time-scale separation assumption, the gene reactions are slow and the protein reactions are fast. Then (47) can be written as:

[TABLE]

where $(f),(s)$ denote fast and slow, respectively.

Hence, the summation in (42) can decomposed into two terms. This implies that the system matrix can be written as a sum of a fast matrix $\tilde{\Lambda}$ and a slow matrix $\varepsilon\hat{\Lambda}$ as in Eq. (14).

We now show that Eq. (15) holds, which amounts to showing that $\tilde{\Lambda}$ is block diagonal. Assume $\exists j$ such that $\tilde{z}=z-\gamma_{j}$. Let $z=(x,d),\tilde{z}=(\tilde{x},\tilde{d})$ with $\tilde{d}\neq d$. As can be seen in Figure 2, protein reactions do not change the promoter configuration state $d$. Hence, the transition rate $\lambda_{z\tilde{z}}$ has terms corresponding to gene reactions rate only, i.e., $\lambda_{z\tilde{z}}=\varepsilon R_{j}^{(s)}(z)$. Hence, $\tilde{\Lambda}$ is block diagonal. $\blacksquare$

2.2 Analytic Expression of the Conditional Probability Distributions

We include a proof of Proposition 2.

As mentioned before, the stationary distribution is the product of the marginal stationary distributions, since the underlying conditional stochastic processes are independent. If $\mathrm{X}_{i}$ does not form a multimer then it is known that the stationary distribution of the reaction network (85) is Poisson with mean $k_{id_{i}}/k_{-i}$ as in (82),

Assume, instead, that $\mathrm{X}_{i}$ forms a multimer. In order to simplify notations, we drop the index $i$ and write $\emptyset\xrightleftharpoons[k_{-}]{k}\mathrm{X},\ n\mathrm{X}\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{X}_{n}.$ Let $x_{1},x_{2}$ denote the molecular counts of $\mathrm{X},\mathrm{X}_{n}$. Then, the master equation is

[TABLE]

We solve the recurrence equation assuming detailed balance, and then verify that the obtained solution, which is given in (82), solves (49) $\blacksquare$

2.3 The Stationary Distribution as a Mixture of Poisson Distributions

We include here the proof of Theorem 3.

Recall the slow-fast decomposition of the master equation in Eq. (14). Recall the joint stationary distribution (46). In order to emphasize the dependence on $\varepsilon$ we denote $\bar{\pi}^{\varepsilon}:=\bar{\pi}(\varepsilon)$. Hence, $\bar{\pi}^{\varepsilon}$ is the unique stationary distribution that satisfies $\Lambda_{\varepsilon}\bar{\pi}^{\varepsilon}=0$, $\pi^{\varepsilon}>0$, and $\sum_{z}\pi_{z}^{\varepsilon}=1$, where the subscript denotes the value of the stationary distribution at $z$.

Our objective is to characterize the stationary distribution as $\varepsilon\to 0$. Writing $\bar{\pi}_{\varepsilon}$ as an asymptotic expansion to first order in terms of $\varepsilon$, we have

[TABLE]

Our aim is to find $\bar{\pi}^{(0)}$. Substituting $\bar{\pi}^{\varepsilon}$ in Eq. (14), and equating the coefficients of the powers of $\varepsilon$ to zero we obtain the following two equations:

[TABLE]

where $\tilde{\Lambda}$ is given in Eq. (15). (53) implies that $\bar{\pi}^{(0)}\in\ker\tilde{\Lambda}$, where $\ker$ denotes the kernel of $\tilde{\Lambda}$. We next show how to compute $\ker\tilde{\Lambda}$.

Recall the conditional Markov chains with the associated infinitesimal generators as in Eq. (16). By the assumptions, for each $d\in\{0,..,L-1\}$ there exists a unique $\pi_{X|d}$ such that: $\Lambda_{d}\pi_{X|d}=0,$ $\pi_{X|d}>0$, and $\sum_{x}\pi_{X|d}(x)=1$. Recall that $\pi_{X|d}$ is the stationary distribution of the Markov chain conditioned on $D(t)=d$.

Defining the extended conditional distributions for $d=0,..,L-1$ as:

[TABLE]

The stationary distribution above can be interpreted as a function as follows: $\bar{\pi}_{X|d}(x,d)=\pi_{X|d}(x)$, and $\bar{\pi}_{X|d}(x,d^{\prime})=0$ when $d^{\prime}\neq d$.

Then $\ker\tilde{\Lambda}=\mbox{span}\{\bar{\pi}_{X|0},..,\bar{\pi}_{X|L-1}\}$. Hence, we can write:

[TABLE]

for some $\lambda_{0},...,\lambda_{L-1}\geq 0$. We normalize them to satisfy $\sum_{d=0}^{L-1}\lambda_{d}=1$.

In order to satisfy (54), we utilize the fact that each $\Lambda_{d}$ is an infinitesimal generator which satisfies $\mathbf{1}^{T}\Lambda_{d}=0$. Hence, we pre-multiply (54) by the vectors: $[\mathbf{1}^{T}\ \mathbf{0}^{T}\ ...\mathbf{0}^{T}]^{T}$, $[\mathbf{0}^{T}\ \mathbf{1}^{T}\ ...\mathbf{0}^{T}]^{T}$, $[\mathbf{0}^{T}\ \mathbf{0}^{T}\ ...\ \mathbf{1}^{T}]^{T}$ in order to get the following $L$-dimensional linear system:

[TABLE]

Furthermore, we need the following normalization equation to find $\lambda_{0},...,\lambda_{L-1}$ uniquely:

[TABLE]

This is equivalent to stating that $\lambda=[\lambda_{0},...,\lambda_{L-1}]$ is the principal eigenvector of $\Lambda_{r}$. $\blacksquare$

2.4 Computation of the Reduced-Order Markov Chain’s Generator

Recall that the generator of the reduced-order Markov chain can be written as follows:

[TABLE]

The $(d^{\prime},d)$ entry represents the probability of transition from the configuration $d^{\prime}$ to configuration $d$, and it can be interpreted as a weighted conditional expectation of $\pi_{d}$.

Consider the reduced chain, and fix a configuration $d$. Then, the maximum number of possible transitions out of $d$ is given by the number of reactions which is $\frac{1}{2}\sum_{i=1}^{N}|B_{i}|$. Hence, $\Lambda_{r}$ is a sparse matrix for large $N$. Computation of the infinite matrices and matrix product in (58) can be cumbersome for networks with multiple genes. Hence, we provide an algorithm for computing the nonzero entries in $\Lambda_{r}$. This can be achieved by considering all the possible transitions from a configuration $d=(d_{1},..,d_{N})$. Specifically, we consider a transition from $d$ to $d^{\prime}$ by a gene reaction modifying a single promoter configuration. For instance consider $D_{d_{i}}$. Then for a constitutive or single TF-gene binding/unbinding, there can be only one transition starting from $D_{d_{i}}$. This transition is either the forward or reverse reaction in (1) or (2), respectively. For the case of two TFs, there can be two reactions among (3)-(6).

The algorithm can be described as follows:

Proposition 6**.**

The matrix $\Lambda_{r}$ in (56) can be computed via the algorithm below.

• •

For each $d\in\{0,..,L-1\}$ write $d=(d_{1},..,d_{N})\in\prod_{i=1}^{N}B_{i}$. Using the previously discussed identification:

Let $\mathscr{R}_{d}=\{\rm R_{1},..,\rm R_{|\mathscr{R}_{d}|}\}$ the set of all gene reactions. Then, for each $j\in\{1,..,|\mathscr{R}_{d}|\}$:

Let $D_{d_{i}}^{i}$, and $D_{d_{i^{\prime}}}^{i}$ be the reactant and product configuration species of the $\rm R_{j}$. Hence, the reaction will cause a transition from $d$ to $d^{\prime}=(d_{1},..,d_{i^{\prime}},..,d_{N})$. Let $\alpha$ be the kinetic constant of $\rm R_{j}$. If $\rm R_{j}$ is a binding reaction, then let $\mathrm{X}_{\bar{i}}$, $\mathrm{X}_{\bar{i}c}$ denote the TF or the multi-merized TF, where $\bar{i}$ denotes the index of the gene that expresses the TF. 2. 2.

Then, the $(d^{\prime},d)$ entry of $\Lambda_{r}$ can be written as:

[TABLE]

• -

Set

[TABLE]

• •

*Set the rest of the entries of $\Lambda_{r}$ to zero. *

Proof. Recall that in Eq. (14), the matrix $\hat{\Lambda}$ represents the slow matrix, which corresponds to the gene binding reactions. Hence, $\hat{\Lambda}_{dd^{\prime}}$ represents the matrix corresponding to the transition between states of the form $(x,d)$ and $(x,d^{\prime})$. Assume that $D(t)=d=(d_{1},...,d_{N})$. Consider the $i$th block. Note that it has one or two reactions that can fire. Specifically, there are $\frac{1}{2}|B_{i}|$ gene reactions that can fire. Assume that such a reaction is in one of the forms:

[TABLE]

where $\mathrm{X}_{\bar{i}}$ is a TF for that block. Then,

[TABLE]

Now consider a reaction of the form:

[TABLE]

Then,

[TABLE]

The last equality follows from evaluating the mean value of the Poisson distribution in (83).

Finally, (60) holds since $\mathbf{1}^{T}\Lambda_{r}=0$, which follows from $\mathbf{1}^{T}\hat{\Lambda}=0$. $\blacksquare$

3 Detailed Discussion of Examples

3.1 The Gene Bursting Model

We start with the simplest form of network, which is the autonomous TF-gene binding/unbinding model. It has been verified as a model for transcriptional bursting [17]. This model has been studied analytically using time-scale separation [32],[62], Poisson-representations [79], and the exact steady solution is known [19], [80].

Consider:

[TABLE]

Referring to Figure 2, we identify a single gene block with two states. Using (82), the conditional stationary distributions are:

[TABLE]

In order to compute the stationary distribution $\pi$, we need to find the generator for the reduced chain. Since both reactions are monomolecular, we write the following using (59):

[TABLE]

Hence, the reduced Markov chain is a binary Bernoulli process with a rate of $\alpha/(\alpha+\alpha_{-})$. Then the stationary distribution of $X$ can be written using (88) as:

[TABLE]

which is a bimodal distribution with peaks at 0 and $k/k_{-}$. The fast promoter kinetics model is obtained, instead, by reversing the time-scale separation such that the protein reactions become slow and gene reactions become fast. In that case, the resulting stationary distribution can be shown to be a Poisson with mean $\frac{\alpha}{\alpha+\alpha_{-}}\frac{k}{k_{-}}$ which is the same as the deterministic equilibrium if we used the conservation law $D_{1}(t)+\mathrm{D}_{0}(t)=1$ for the above model. Finally, note that the mean of the slow promoter kinetics model is the same as in the fast kinetics model but the two stationary distributions differ drastically.

Figure 7 shows the transition from fast to slow promoter kinetics using the exact solution [80] and compares it to the expression (62).

3.2 A Self-Regulating Gene

Consider a non-cooperative self-regulating gene Eq. (91).

Referring to Figure 2, this is a single-gene block with two states. At the limit of slow promoter kinetics, Remark 4 implies that the gene binding/unbinding reaction can be written as in Eq. (29) as follows:

[TABLE]

Using (59), the reduced generator can be written as:

[TABLE]

Hence it defines a binary Bernoulli process with the rate $\alpha k_{0}/(\alpha k_{-}+\alpha k_{0})$. Using (88) the stationary distribution is a mixture of two Poisson distributions and can be written as:

[TABLE]

where

[TABLE]

Next, consider the same reaction network, but now with cooperativity:

[TABLE]

In this case, the gene process is still a Bernoulli process, but with a different rate. The stationary distribution for $X$ can be written as:

[TABLE]

where

[TABLE]

Both distributions (63), (65) have modes at $\frac{k_{1}}{k_{-}}$ and $\frac{k_{0}}{k_{-}}$. The height of the first mode is proportional to $\rho_{1}$ for (63), and is proportional to $\rho_{2}$ for (65). The network is activating if $k_{1}>k_{0}$, and repressing otherwise.

Comparing (63) and (65), note that, in the non-cooperative case, if we want to increase the weight of the mode corresponding to the bound state keeping the association ratio, then the mode location needs to be changed. On the other hand, the factor $\rho_{2}$ in the dimerization rates in (65) can be used in order to tune the weights freely while keeping the modes and the binding to unbinding kinetic constants ratio unchanged. For instance, we can make the distribution effectively unimodal with a sufficiently high dimerization ratio.

A non-cooperative self-regulating gene with slow promoter kinetics has been studied in literature by deriving close-form expression [31] and using time-scale separation [32] . However, the gene binding/unbinding reaction in both papers was approximated by an auto-catalytic reaction:

[TABLE]

This has the advantage of decoupling the slow and the fast processes. However, it is a simplification of the physical process. We did not use such simplifications.

Special cases:

The model above considers a network with possibly non-zero production rates for both the unbound and bound promoter configurations. We may also consider the special cases of pure self-activation or self-repression which are explained below:

Pure Self-Activation, i.e. $k_{0}=0$ in (91) and (64). Then the peak corresponding to the bound configuration disappears and we get only one peak at zero protein copy number as it forms an absorbing state. This is a manifestation of the Keizer’s paradox [81]. One way to circumvent this is to allow for a small transcriptional “leak”. This amounts to taking $k_{0}\ll k_{1}$, and it allows us to recover the second mode. 2. 2.

Pure Self-Repression, i.e. $k_{1}=0$ in (91) and (64). Then we get two modes: one at 0 and the other at $k_{0}/k_{-}$.

For both cases, in the non-cooperative case with a fixed dissociation ratio the choice of this kinetic rate determines completely the relative weight of the modes as in (63). Cooperativity allows us to tune the weight of the mode corresponding to the bound state without changing the location as mentioned before.

Comparison with Fast Promoter Kinetics:

In order to demonstrate that slow switching is responsible for the emergence of new modes compared to the deterministic model, consider the non-cooperative self-regulating gene network (91) with fast promoter kinetics modeled by letting $\varepsilon$ grow without bound in the first reversible TF-promoter binding/unbinding reactions. We state the following proposition which is proved in the Methods section:

Proposition 7**.**

As $\varepsilon\to\infty$, the stationary distribution of the network (91) is given by:

[TABLE]

where $w_{m}$ satisfies the following recurrence relation:

[TABLE]

where $w_{0}$ is chosen to satisfy $\sum_{m=0}^{\infty}w_{m}=1$.

Proof.

A decomposition dual to Eq. (14) can be written and it can be noted that the fast matrix is block-diagonal with respect to $\mathrm{X}+\mathrm{D}_{1}$ which is the slow variable, while the fast variable is $\mathrm{D}_{1}$.

Expanding asymptotically the stationary distribution in terms of $\varepsilon$, and taking the limit as $\varepsilon$ goes to zero, we can find the distribution for the slow variable as follows:

[TABLE]

where $w_{m}$ satisfies the following recurrence relation:

[TABLE]

The joint distribution can be given as:

[TABLE]

Hence we can compute the marginal density of $X$ as follows:

[TABLE]

∎

Since the ratio $w_{m+1}/w_{m}$ is a ratio of two polynomials and the denominator’s degree is higher than the numerator, then stationary distribution is unimodal, while slow TF-gene binding/unbinding was shown to give a bimodal distribution (see (63)).

3.3 The Toggle Switch

A toggle switch is a basic GRN that exhibits deterministic multi-stability. It has two stable steady states and can switch between them with an external input or via noise. The ideal behavior is that only one gene is “on” at any moment in time. We now study the network with the slow promoter kinetics. Consider the following network with cooperativity indices $n,m$:

[TABLE]

For the case $n,m=1$, there is no multi-merization reaction. For consistency, we choose $\beta_{1}=\beta_{-1},\beta_{2}=\beta_{-2}=1$ in that case.

Denote the promoter configuration species by $\mathrm{D}^{X},D^{Y}$. Then the network has four configurations $(\mathrm{D}^{X},\mathrm{D}^{Y})\in\{(0,0),(0,1)$ $,(1,0),(1,1)\}$. Using Corollary 4 we expect to have a stationary distribution with four modes $(k_{10}/k_{-1},k_{20}/k_{-2})$, $(k_{10}/k_{-1},0)$, $(0,k_{20}/k_{-2})$, $(0,0)$. Using the algorithm of Proposition 5, the reduced-order Markov chain infinitesimal generator is:

[TABLE]

where

[TABLE]

We notice immediately from the last row in the matrix (66) that the transition rates towards the configuration (1,1) are zero, which implies that the weight of the mode corresponding to $(1,1)$ is zero. Hence, we have three modes only. The weights corresponding to the modes can be found as the principal eigenvector of $\Lambda_{r}$ as given in Corollary 4. Hence, the stationary distribution for $X,Y$ is:

[TABLE]

Since the stationary distribution has three modes, it deviates from the ideal behavior of a switch where at most two stable steady states, under appropriate parameter conditions, are possible. Nevertheless, a bimodal distribution can be achieved by minimizing the weight of the first mode at $(\frac{k_{10}}{k_{-1}},\frac{k_{20}}{k_{-2}})$. If we fix $\alpha_{1}/\alpha_{-1},\alpha_{2}/\alpha_{-2}$, then this can be satisfied by tuning $n,m,\beta_{\pm 1},\beta_{\pm 2}$ to maximize $\rho_{1},\rho_{2}$ in (67). Choosing higher cooperativity indices, subject to $n<k_{10}/k_{-1},m<k_{20}/k_{-2}$, achieves this.

The toggle switch has three modes regardless of the cooperativity index. This is unlike the deterministic model where only one positive stable state is realizable with non-cooperative binding, and two stable steady states are realizable with cooperative binding. However, the toggle switch with fast switching can admit three modes in some parameter ranges. In contrast to the case of slow switching under consideration here, the third stable state is the (low,low) state [82].

Monte-Carlo simulations via Gillespie’s algorithm for the cooperative toggle switch have been performed to investigate the minimum time-scale separation to recover our predictions. We have chosen the parameters such that the (low,high), and (high,low) modes get 0.25 each, and the (high,high) modes gets 0.5 at the slow promoter limit. Figure 8-a shows that the three modes predicted by the analysis have been recovered with promoter kinetics being just two times slower than the protein decay rate. Figure 8-b shows a sample trajectory where the three modes are visible.

Alternative model of the toggle switch

Instead of modeling the toggle switch as two independent genes with their own promoters, an alternative model consists of a signal promoter regulating two operons [83], [84]. In our modelling framework this amounts to a single gene expression block with two transcription factors and two expressed proteins, i.e., two inputs and two outputs. Despite the fact that our formalism accounts for a single expressed protein, we find no difficulty in applying our methods as shown below.

The alternative model can be written as:

[TABLE]

Assuming slow-promoter kinetics, the stationary distribution can be decomposed into a mixture of Poisson distributions centered at $(\frac{k_{x}}{k_{-x}},\frac{k_{y}}{k_{-y}}),(\frac{k_{x}}{k_{-x}},0),(0,\frac{k_{y}}{k_{-y}})$. Despite the fact that model appears to be different, the resulting modes and the reduced-order Markov chain are similar to the model used before.

3.4 The Repressilator

The repressilator is a synthetic biological circuit that implements a ring oscillator [51], and it has been simulated with slow-promoter kinetics [21]. It is a canonical example of a GRN that exhibits a limit cycle, i.e. sustained oscillation. The block diagram is shown in Figure 9-a which shows three genes in a repression cycle. The list of reactions is given as follows:

[TABLE]

Deterministic analysis of the repressilator [3] reveals that it does not oscillate with non-cooperative binding. We proceed to apply our analysis methods to study the behaviour in the stochastic case with slow promoter kinetics.

Denote the promoter configuration species by $\mathrm{D}^{X},\mathrm{D}^{Y},\mathrm{D}^{Z}$. Then the network has eight configurations $(\mathrm{D}^{X},\mathrm{D}^{Y},\mathrm{D}^{Z})\in\{(0,0,0),(0,0,1)$ $,...,(1,1,1)\}$. Using Corollary 4 we expect to have a stationary distribution with eight modes. Using the algorithm of Proposition 5, the reduced-order Markov chain infinitesimal generator is:

[TABLE]

where $\rho=\left(\frac{k}{k_{-}}\right)^{n}\frac{\beta}{n!\beta_{-}}.$

The stationary distribution is:

[TABLE]

where $K=\frac{k}{k_{-}},w=2\rho\frac{\alpha}{\alpha_{-}},D=3w+3+\tfrac{1}{2w}$. Notice that the Poisson distribution centered at (0,0,0) has zero weight. Hence, the network can admit up to seven modes.

We will study the behaviour of this network in two cases: when $w\gg 1$, and $w\approx 1$.

The “Stochastic Oscillator”

Although the distribution is a mixture of seven modes, note that it reduces to three effective modes if

[TABLE]

where the three modes are located at $d_{x}=(K,0,0),d_{y}=(0,K,0),d_{z}=(0,0,K)$. The tri-modal stationary distribution is consistent with the classical oscillations of the repressilator, and this is independent of the cooperativity index. Note that this condition is analogous to the oscillation condition in the deterministic model [3] (but with cooperativity only) which also requires “large” production ratio.

In order to study whether the network oscillates, we need to define a notion of limit cycle for a stochastic system. Due to randomness, the time-series can not be periodic. Nevertheless, since the stationary distribution is tri-modal, we say that the network oscillates if the time trajectory continues to jump between the modes in the same order.

Assume $w\gg 1$. Let $d_{x},d_{y},d_{z}$ be the three dominant modes. We will show that if the reduced-order Markov chain is at mode $r_{x}$ then it is much more likely to transition to $d_{y}$ rather than to $d_{z}$. Similar arguments apply if we start from $d_{y},d_{z}$.

Let $Q(t)=e^{t\Lambda_{r}}$ be the probability transition matrix [11]. The $(i,j)$ entry in the transition matrix denotes the probability of being at $i$ at time $t$ if we start from $j$. Hence, $Q_{ij}(t)=\Pr[D(t)=i|D(0)=j]$.

We are interested in comparing the probabilities of transiting from $r_{x}$ to $r_{y},r_{z}$. Hence we study small $t\ll 1$. We utilize the power series expansion to evaluate $Q(t)$. Let $\Lambda_{d}=\mbox{diag}(\Lambda_{r})$ and $\Lambda_{+}=\Lambda_{r}-\Lambda_{d}$, where $\mbox{diag}(\Lambda_{r})$ is the diagonal matrix which contains the diagonal entries of $\Lambda_{r}$. Note that $\Lambda_{+}$ is positive. Since $\Lambda_{d}$ and $\Lambda_{+}$ commute, we can write

[TABLE]

We can approximate $Q(t)$ for small $t$ by writing the Taylor expansion of $e^{\Lambda_{+}}$. Hence, using (70) we write a third-order Taylor series for the transition probabilities that we need:

[TABLE]

where we assumed, w.l.o.g, that $\alpha_{-}=1$.

Using the expressions above, if the Markov chain is at $r_{x}$ then it is most likely to stay there. The transition is much more likely to happen to $d_{y}$ rather than $d_{z}$. Hence, we expect to see “long” periods of protein $X$ being expressed, and then it jumps to express protein $Y$, and then protein $Z$. Since the finite Markov chain is ergodic, the pattern repeats.

Note that the analysis above predicts that both the cooperative and the noncooperative repressilator are capable of oscillation with slow-promoter kinetics when $w\gg 1$. We performed Monte-Carlo simulations via the Gillespie algorithm for both fast and slow kinetics. The results are shown in Figure 10. We observe that the network always oscillates with cooperative binding. With non-cooperative binding, only the network with slow kinetics oscillates as predicted. The network with fast kinetics does not oscillate. Recall that the deterministic model with non-cooperative binding does not oscillate [3].

Multi-modality in the repressilator with slow-promoter kinetics

In the case that $w$ is close to 1, the seven modes in (71) share the probability almost equally. This is independent of the cooperativity index. Compare this the non-oscillatory deterministic model, which is mono-stable and it can’t admit multiple stable equilibria.

3.5 Synchronization of interconnected toggle switches

We consider $N$ identical toggle switches:

[TABLE]

where $i=1,..,N$. We interconnect these systems through diffusion of the protein species $\mathrm{X}_{i},\mathrm{Y}_{i}$ among cells, modeled through reversible reactions with a diffusion coefficient $\Omega$:

[TABLE]

For a deterministic model, there exists a parameter range for which all toggle switches will synchronize into bistability for sufficiently high diffusion coefficient [44]. This implies each switch in the network behaves as a bistable switch, and it converges with all the other switches to the same steady-states.

Our aim is to analyze the stochastic model at the limit of slow promoter kinetics and compare it with the deterministic model. The network has $4^{N}$ promoter configurations. Consider a promoter configuration

[TABLE]

where $d_{i}^{X},d_{j}^{Y}\in\{0,1\}$. As discussed before, we need to represent the Markov chain conditioned on $D(t)=d$. Each conditional Markov chain can be represented as follows:

[TABLE]

This conditional Markov is not the forms (85), (86) that arise from the class of GRNs defined previously as in Figure 1. Nevertheless, it can be observed that it is a reversible zero-deficiency network [85]. We show in the SI-§4.1 that our results can be generalized to networks that admit weakly reversible deficiency zero conditional Markov chains.

As shown in the SI-§4.1, there are $4^{N}$ conditional Markov chains, and their stationary distributions are Poisson distributions with the following set of modes:

[TABLE]

where

[TABLE]

in which $\bar{d}_{i}^{X}=1-d_{i}^{X}$. $\omega_{d_{i}}^{Y}$ is defined similarly.

Using Theorem 3, the stationary distribution is a mixture of $4^{N}$ Poisson distributions, and the weights can be found by finding the principal eigenvector of the infinitesimal generator of the reduced-order Markov chain as before. Proposition 5, gives the procedure to find the matrix $\Lambda_{r}$, where Eq. (27) is replaced by (78) with (82). As before, the Markov states are the configurations $d=(d^{X},d^{Y})\in\{0,..,4^{N}-1\}$. Assume that $D(t)=d$. Then, the state transitions are given by the following reactions:

[TABLE]

where

[TABLE]

and $\rho_{d_{i}}^{Y}$ is defined analogously.

Note that the mode corresponding to the state in which all the TFs are bound to the promoters has no incoming transitions in the reduced Markov chain and hence it has zero weight, hence the network has $4^{N}-1$ modes.

We consider now the case of a high diffusion coefficient. Note from (73), (74), (76) that as $\Omega\to\infty$, $X_{1},..,X_{N}$ will synchronize in the sense that the joint distribution of $X_{1},..,X_{N}$ is symmetric with respect to all permutations of the random variables. This implies that the marginal stationary distributions $p_{X_{i}},i=1,..,N$ are identical. Hence, for sufficiently large $\Omega$ the probability mass is concentrated around the region for which $X_{1},..,X_{N}$ are close to each other. Consequently, for large $\Omega$ we can replace the population of toggle switches with a single toggle switch with the synchronized protein processes $X(t),Y(t)$, which are defined, for the sake of convenience, as $X(t):=X_{1}(t),Y(t):=Y_{1}(t)$. Next, we describe the stationary distribution of $X(t),Y(t)$.

From (74) it can be seen that $\omega_{d_{i}}^{X}$ does not depend on $d_{i}$ for large $\Omega$. Instead it depends only on $\sum_{i=1}^{N}\bar{d}_{i}^{X}$, which is the total number of unbound promoter sites in the genes producing $\mathrm{X}_{1},..,\mathrm{X}_{N}$. The same holds for $\omega_{d_{i}}^{Y}$. Hence, the number of modes will drop from $4^{N}-1$ to $(N+1)^{2}-1$. Hence the joint distribution of $X,Y$ is a mixture of Poisson distributions with the following modes:

[TABLE]

Note that similar to the single toggle switch, there are modes which have both $X,Y$ with non-zero copy number. On the other hand, there are many additional modes. Recall that in the case of a single toggle switch, we have tuned the cooperativity ratios such that the modes in which both genes are ON are suppressed. Similarly, the undesired modes can be suppressed by tuning the cooperativity ratio which can be achieved by choosing $\rho_{d_{i}}^{X},\rho_{d_{i}}^{Y},d=0,..,4^{N}-1$ sufficiently large. In particular, letting the multi-merization ratio $\beta_{x}/\beta_{-x},\beta_{y}/\beta_{-y}\to\infty$, the weights of modes in the interior of the positive orthant $\mathbb{R}_{+}^{2}$ approach zero.

In conclusion, for sufficiently high $\Omega$ and sufficiently high multimerization ratio the population behaves as a *multimodal switch *, which means that the whole network can have either the gene $\mathrm{X}$ ON, or the gene $\mathrm{Y}$ ON. And every gene can take $2N$ modes which are:

[TABLE]

Comparing to the low diffusion case, the network will have up to $2^{N}-1$ modes with sufficiently high multimerization ratio.

Comparison with the deterministic model

For $\Omega$ greater then a certain threshold, the deterministic system bifurcates into bistabiliy. This means that all toggle switches converge to the same exact equilibria if $\Omega$ is greater than the threshold. As we have seen before, this is not the case for the stochastic system, since the toggle switches converges asymptotically to each other. Hence, we need to choose a a threshold for $\Omega$ that constitutes “sufficient” synchronization. We choose to define this as the protein processes synchronizing within one copy number. In other words, we require the maximum distance between the modes in (73) to be less than 1.

This amounts to requiring the diffusion coefficient needs to satisfy:

[TABLE]

We derive (77) below.

In order for the maximum distance between the modes in (73) to be less than 1, we need to satisfy this inequality:

[TABLE]

Let us consider first maximizing the term containing the variables related to the gene $\mathrm{X}$. Note from (74) that $\omega_{d_{i}}^{X}$ and $\omega_{d_{j}}^{X}$ can differ only in a $k_{-x}/\Omega$ term in the numerator. Hence, the maximization can be simplified to:

[TABLE]

Solving for $\Omega$ yields (77).

4 Extension to GRNs with Complex-Balanced Conditional Markov Chains

In the main text we have included assumptions to simplify the mathematical treatment and notations while retaining the same qualitative features of the problem. In this section, we elucidate the manner in which the results can be generalized.

The class of networks defined in the main text does not allow direct interactions between the proteins. However, our main results hold for a more general class of networks. Figure 2 depicts a generalized gene expression block where the proteins expressed by other genes participate in the gene reactions block. The gene expression block remains unchanged. The basic theory presented in the main text (Proposition 1, Theorem 3) remains unchanged. Proposition 5 needs can be modified by replacing Eq. (27) by the following:

[TABLE]

However, since there are no closed form formulae for computing the conditional expectation in general, the expression above is of limited utility. Nevertheless, we can define a more general class of network such that the conditional stationary distributions can be computed analytically as obtained in Proposition 2. In order to define our generalized class of networks, we need some further notation and background which are summarized in the following subsection.

4.1 Complex-Balancing and the Concept of Deficiency

Recall that a reaction network consists of the set of species $\mathscr{S}$ and the set of reactions $\mathscr{R}$. A reaction $\mathrm{R}_{j}\in\mathscr{R}$ can be written as Eq. (40).

The formal linear combination in the left-side of the reaction is called the reactant complex, while the linear combination on the right-side of the reaction is called the product complex. The set of all complexes in the network is denoted by $\mathscr{C}$. Each complex can be interpreted as a vector belonging to the Euclidean vector space whose unit vectors correspond to the species.

Define the following matrix:

[TABLE]

where $\alpha_{ij},\beta_{ij}$ are the stoichiometry coefficients. If we remove the redundant columns of $\tilde{\Gamma}$, this yields an $|\mathscr{S}|\times|\mathscr{C}|$ matrix $\Gamma_{c}$, where each column corresponds to a certain complex. Hence, it can be verified that there exists a matrix $Y\in\{0,\pm 1\}^{\mathscr{|}C|\times|\mathscr{R}|}$ such the that the stoichiometry matrix $\Gamma$ admits the following factorization:

[TABLE]

There is significant literature studying the reaction networks using the notion of complexes, initially for deterministic systems, [86, 85], and then for stochastic systems [87, 88]. The ordinary differential equation for a reaction network can be written as:

[TABLE]

where $x$ is the concentration vector, and $R$ is a reaction rate function which is commonly defined using mass-action kinetics as follows:

[TABLE]

where $k_{j}$ is reaction rate constant. A conservation law is a nonnegative vector $v$ such that $v^{T}\Gamma=0$. The stoichiometric compatibility class containing $x_{0}$ is defined as $(x_{\circ}+\Im\Gamma)\cap\mathbb{R}_{\geq 0}^{n}$. If there are no conservation laws then $\mathbb{R}_{\geq 0}^{n}$ is a stoichiometric compatibility class.

In order to state the subsequent results, we define a modified mass-action kinetics. For a given a reaction rate function $R$, the modified reaction rate function is given as

[TABLE]

As per (79), the complex formation rate is defined as $R_{c}(x):=YR(x)$. An equilibrium $\hat{x}$ of (80) is called complex balanced if it satisfies $R_{c}(x)=0$ also. The existence of a single complex-balanced equilibrium guarantees that all equilibria of the network are complex balanced [86].

Physically, reaction rates in $\ker Y$ do not change the complex-formation rate. Hence, if $\ker{\Gamma}=\ker{Y}$, then every equilibrium is complex-balanced. This ensures that a network is complex balanced regardless of kinetic rate constants. Hence, the deficiency of the network is defined as $\delta:=\mbox{dim}(\ker(\Gamma_{c})\cap\mbox{Image}(Y))$. If $\delta=0$ then every equilibrium is complex-balanced, and the network is said to be a zero-deficiency network.

A network is weakly reversible if: existence of a directed path from complex $C_{i}$ to complex $C_{j}$, implies the existence a directed path from $C_{j}$ to $C_{i}$. This is equivalent to the existence of a strictly positive vector in $\ker Y$.

For weakly-reversible networks, a simple graphical characterization of the deficiency is given by [85]:

[TABLE]

where $\ell$ is the number of strongly connected components in the graph of complexes.

A main result in the theory of deterministic complex-balanced networks is stated as:

Lemma 8** ([85]).**

If a reaction network with mass-action kinetics is weakly reversible and has deficiency zero, then there exists a unique positive equilibrium in every stoichiometric compatibility class.

Hence if the network is complex-balanced and without conservation laws, then it has a unique equilibrium for the deterministic system. A parallel result exists for the Markov chain model with associated chemical master equation. It is stated as follows:

Lemma 9** ([87, 88]).**

Assume a network is weakly reversible, has deficiency zero and has no conservation laws. Let $\hat{x}=[\hat{x}_{1},..,\hat{x}_{n}]^{T}$ be the unique complex balance equilibrium for the deterministic system with the modified mass-action kinetics. Then the stationary distribution for the corresponding master equation is:

[TABLE]

4.2 Extension of the Main Result

We are ready to define the generalized class of networks as gene regulatory networks whose conditional Markov chains are complex-balanced. In particular, if all the conditional Markov chains are weakly reversible and zero-deficiency, then the GRN can be treated with our methods.

Consider a set of $N$ genes, binding sets $\{B_{i}\}_{i=1}^{N}$, and kinetics constants $k_{j}$. A generalized gene expression block is as shown in Figure 11. Then, a gene regulatory network is an arbitrary interconnection of gene expression blocks subject to Assumption A2. A gene regulatory network admits conditional product-form distributions if, for each $d\in\{0,..,L-1\}$, the reaction network corresponding to the conditional Markov chain is weakly reversible and deficiency zero. Hence, we can restate Proposition 2:

Proposition 10**.**

Fix $d\in\{0,..,L-1\}$. Let $\hat{x}^{d},\hat{x}^{d2}$ be the complex-balanced equilibria for the proteins and multi-merized proteins, respectively. Consider eq. (16), then there exists a conditional joint stationary distribution $\pi_{X|d}^{(J)}$ and it is given by

[TABLE]

where

[TABLE]

where $(J)$ refers to the joint distribution in multimerized and non-multimerized processes, $x_{i1}$ refers to the copy number of $X_{i}$, while $x_{i2}$ refers to the copy number of $X_{ic}$, $\mathbf{P}(x;a):=\frac{a^{x}}{x!}e^{-a},\mathbf{P}(x_{1},x_{2};a_{1},a_{2}):=\frac{a_{1}^{x_{1}}}{x_{1}!}\frac{a_{2}^{x_{2}}}{x_{2}!}e^{-a_{1}-a_{2}}$.

4.3 Specific Examples

4.3.1 Cooperative Binding

Cooperative binding is the process in which a TF forms as a hetero-dimer consisting of two different proteins. Consider two proteins $\mathrm{X},\mathrm{Y}$, then the hetero-dimerization reaction:

[TABLE]

Consider a conditional Markov chain that consists of uncoupled birth-death processes, with multi-merization reactions. It is reversible and has deficiency zero. If a hetero-dimerization reaction is added, then this involves adding two complexes, adding a new strongly connected component to the complexes’ graph, and adding a new linearly-independent row to the stoichiometry matrix. According to (81) this will not change the deficiency. Hence, our results can be extended immediately to networks with cooperative binding.

4.3.2 Diffusion

We model diffusion between two gene expression blocks as a reversible reaction between the proteins as:

[TABLE]

Consider a conditional Markov chain that consists of uncoupled birth-death processes, with multi-merization reactions. It is reversible and has deficiency zero. If a diffusion reaction is added, then none of the number of complexes, the number of strongly connected components, or the rank of the stoichiometry matrix change. According to (81) this will not change the deficiency. Hence, our results can be extended immediately to networks with diffusion. We have already applied this to the communicating toggle switches in the main text.

4.3.3 Multi-Step Multi-merization

In the main text we have assumed that the multi-merization process occurs in one-step only as follows:

[TABLE]

As an alternative model, consider the following multi-step or sequential multi-merization process:

[TABLE]

Consider a conditional Markov chain, conditioned on $D=d$, that consists of the above network with a birth and death process:

[TABLE]

In order to evaluate the deficiency (81) of the network, note that $|\mathscr{C}|=2+2n,\ell=n+1,\mbox{rank}(\Gamma)=n+1$. Hence, $\delta=0$. Therefore, the network is complex balanced and the stationary distribution is a product of Poisson distributions.

Since $X_{n}$ acts as a transcription factor we are interested in finding $\mathbb{E}[X_{n}|D=d]$. Using Lemma SI-2, we can find the required quantity by finding the a complex-balanced equilibrium of the associated deterministic system with modified mass-action kinetics. Hence, we find:

[TABLE]

From the perspective of the slow promoter kinetics, the multi-step multi-merization reactions can be replaced with a single reaction (84) with $\beta/\beta_{-}=\frac{1}{2}n!\prod_{i=2}^{n}{\beta_{i}}/{\beta_{-i}}$.

5 Extension to Networks with Multiple Copies of the Genes

In the main text, we have assumed that there exists one copy only per gene. Nevertheless, our theory is not limited by this assumption. In this subsection we show how the theory can be generalized.

Consider the network in Figure 2. Assume that the $i$th gene has $M_{i}$ copies. Let $B_{i}$ be the corresponding binding-set, and denote $b_{i}=|B_{i}|-1$. We need to track the promoter configurations for all copies of the $i$th gene. Hence, we define the total binding set $P_{i}$. The new set $P_{i}$ describes all possible $|B_{i}|$-tuples of nonnegative integers such that they sum to $M_{i}$. It can be shown that $|P_{i}|=\binom{M_{i}+b_{i}}{b_{i}}$. To explain this, consider the following:

• •

If one TF binds to a promoter or the promoter changes its configuration autonomously, then $B_{i}=\{0,1\}$, while $P_{i}:=\{(0,M_{i}),(1,M_{i}),..,(M_{i},0)\}$. Hence $|P_{i}|=M_{i}+1$.

• •

If two promoters bind independently to a promoter then $B_{i}=\{00,01,10,11\}$, while $P_{i}:=\{(M_{i},0,0,0),(M_{i}-1,1,0,0),....,(0,0,0,M_{i})\}$, and $|P_{i}|=(M_{i}+3)(M_{i}+2)(M_{i}+1)/6$. Each tuple $(d_{i0},..,d_{i3})$ is interpreted as follows: there are $d_{i0}$ copies of the gene at configuration 00, $d_{i1}$ at configuration 01, etc. The tuple entries are ordered to correspond to 00,01,10,11, respectively.

• •

If two TFs bind competitively to a promoter then the binding set is $B_{i}:=\{00,10,01\}$, while $P_{i}:=\{(M_{i},0,0),(M_{i}-1,1,0),....,(0,0,M_{i})\}$, and $|P_{i}|=(M_{i}+2)(M_{i}+1)/2$. The tuple entries are ordered to correspond to 00,01,10, respectively.

Consider the master equation (42). Let us consider that the $i$th gene has $M_{i}$ copies. We again split the stochastic process $Z(t)$ into two subprocesses: the gene process $D(t)$ and the protein process $X(t)$. We will redefine $D(t)$ as explained below.

Consider the $i^{\rm th}$ gene. We define the process $D_{i}$ such that $D_{i}(t)\in P_{i}$, so it informs us exactly the copy numbers of genes at each promoter configuration. Collecting these into a vector, define the gene process $D(t):=[D_{1}(t),...,D_{N}(t)]^{T}$ where $D(t)\in\prod_{i=1}^{N}P_{i}$. The $i^{\rm th}$ gene can be represented by $|P_{i}|$ states, so $L:={\prod_{i=1}^{N}|P_{i}|}$ is the total number of configurations in the GRN. With abuse of notation, we write also $D(t)\in\{0,..,L-1\}$ in the sense of some fixed bijection between $\{0,..,L-1\}$ and $\prod_{i=1}^{N}P_{i}$. Hence, $d\in\{0,..,L-1\}$ corresponds to $(d_{1},...,d_{N})\in P_{1}\times..\times P_{N}$ and we write $d=(d_{1},..,d_{N})$. Note that $d_{i}=(d_{i0},..,d_{ib_{i}})$. With the redefined $D$, the protein processes and the master equation are defined in an identical manner to the main text.

In order to apply Proposition 2, we need to reexamine the conditional Markov chains. Fix $D(t)=d$, and consider the $i$th gene. Let $d_{i}=(d_{i0},..,d_{ib_{i}})$, which informs us the copy number of genes with a certain configuration. Hence, the production rate for the protein is the sum of the production rates of the genes weighted by their respective copy numbers. Then, the conditional Markov chain is given as:

[TABLE]

or, with multimerization, as:

[TABLE]

where

[TABLE]

Theorem 3 and Corollary 4 hold as in the main text. Hence we conclude that

[TABLE]

In order to find the weights $\lambda_{i}$’s, we need to restate Proposition SI-2 as follows:

Proposition 11**.**

The matrix $\Lambda_{r}$ in eq. (56) can be computed via the algorithm below.

• •

For each $d\in\{0,..,L-1\}$ write $d=(d_{1},..,d_{N})=((d_{10},...,d_{1b_{1}}),..,(d_{N0},...,d_{Nb_{N}}))\in\prod_{i=1}^{N}P_{i}$. Using the previously discussed identification:

Let $\mathscr{R}_{d}=\{\rm R_{1},..,\rm R_{|\mathscr{R}_{d}|}\}$ be the set of all gene reactions. Then, for each $\nu\in\{1,..,|\mathscr{R}_{d}|\}$:

Let $\mathrm{D}_{j}^{i}$, and $\mathrm{D}_{j^{\prime}}^{i}$ be the reactant and product configuration species of the $\rm R_{\nu}$, where $j\in B_{i}$. The reaction will cause a transition from $d$ to a state $d^{\prime}$ in which $d_{ij}^{{}^{\prime}}=d_{ij}-1$ $d_{ij^{\prime}}^{{}^{\prime}}=d_{ij^{\prime}}+1$, and $d_{ij}^{{}^{\prime}}=d_{ij}$ otherwise. Let $\alpha$ be the kinetic constant of $\rm R_{\nu}$. If $\rm R_{\nu}$ is a binding reaction, then let $\mathrm{X}_{\bar{i}}$, $\mathrm{X}_{\bar{i}c}$ denote the TF or the multi-merized TF, where $\bar{i}$ denotes the index of the gene that expresses the TF. 2. 2.

Then, the $(d^{\prime},d)$ entry of $\Lambda_{r}$ can be written as:

[TABLE]

where $k_{\bar{i}d_{\bar{i}}}$ is defined in (87).

• -

Set

[TABLE]

• •

*Set the rest of the entries of $\Lambda_{r}$ to zero. *

Example: The Self-Regulating Gene

Let us consider again the non-cooperative self-regulating gene,

[TABLE]

Assume that the gene has $M$ copies. Hence, the total binding set $P_{i}=\{(M,0)...,(0,M)\}$, and we have $M+1$ states. Using Proposition 11 the reduced order Markov chain can be represented as follows:

[TABLE]

The resulting stationary distribution is a mixture of Poissons with modes at $Mk_{0},(M-1)k_{0}+k_{1},...,Mk_{1}$. Figure 12 shows the stationary distribution for several values of $M$ for the network above. The number of modes increases with $M$, however, the stationary distribution converges for high $M$.

6 Additional Information for Tables and Figures in the Main Text

6.1 Figure 1-b

A self-repressing gene is given as:

[TABLE]

The ordinary differential equations for the deterministic model can be written as:

[TABLE]

with $d_{0}(t)+d_{1}(t)=1$ for all $t$.

Setting the derivatives to zero, and solving for $x$ gives the following charcterisitic equation:

[TABLE]

Since the third term has a negative sign, the network has a unique positive equilibrium for all parameter values. The parameters used in the simulation are:

[TABLE]

The deterministic equilibrium is $\hat{x}\approx 34.1641$. The curve for the fast promoter kinetic has been computed by applying the recurrence relation in Proposition 7. The curve for the slow promoter kinetic has been computed by Eq. (32). The remainder of the curves have been computed using a finite state projection truncated at $x=65$. This amounts to writing the master equation as a finite linear differential equation of the form:

[TABLE]

then evaluating the principal eigenvector of $\Lambda$ using the command null in MATLAB, and normalizing the resulting vector.

6.2 Figure 1-d

A hybrid repression-activation network is given as:

[TABLE]

The parameters used in the simulation are:

[TABLE]

Modelling the network deterministically, there exists a unique positive equilibrium where the equilibrium values of $X,Y$ are given as $\hat{x}\approx 25.4727,\hat{y}\approx 39.9216.$ The stationary distribution has been computed using Corollary 4 and Proposition 5.

6.3 Table 1

Let us consider the non-cooperative self-regulating gene. The reaction network is given as:

[TABLE]

The deterministic model gives the following characteristic equation for $x$:

[TABLE]

Since the third term has a negative sign, the network has a unique positive equilibrium for all parameter values. For the stochastic model: With no leakiness, i.e. $k_{0}=0$, the Markov chain has an absorbing state at $(X,D_{1})=(0,0)$. Hence, the network has a single mode at zero. With leakiness, we have shown in the main text that the slow promoter kinetics give rise to two modes (Corollary 4). While for fast promoter kinetics, we have mentioned that it gives rise to a uni-modal distribution after Proposition 7. Let us consider the cooperative self-regulating gene. The reaction network is given as:

[TABLE]

The deterministic model gives the following characteristic equation for $x$:

[TABLE]

Since the constant term is negative, this means that (96) has at least one positive real solution. The equation (96) has three solutions if

[TABLE]

and one solution otherwise. Using the Routh-Hurwitz cireterion, all solutions have positive real parts if and only if $k_{1}>k_{0}$. This means that the network can have either one or two positive stable equilibria only.

For the stochastic model: with no leakiness, the discussion is identical to the previous case. With leakiness, we have shown in the main text that the slow promoter kinetics give rise to two modes independent of the cooperativity index (Corollary 4). While for fast promoter kinetics, a single mode has been verified by solving the master equation numerically for many parameter sets including the ones that are multi-stable for the deterministic system.

6.4 Figure 3-b

The reaction network is given in Eq. (34). The parameters are:

[TABLE]

For these parameters, the deterministic system has a unique stable positive equilibrium at $\hat{x}\approx 9.9195$. The curve for the slow promoters limit has been computed using Eq. (36). The remaining curves have been computed by a numerical solution of the master equation as articulated in SI-5.1.

6.5 Figure 5b-e

We consider a population of three toggle switches ($N=3$) with $k_{x}=k_{y}=150,k_{-x}=k_{-y}=1,\beta_{x}=\beta_{y}=1,\beta_{-x}=\beta_{-y}=1,\alpha_{x}=\alpha_{y}=0.3\varepsilon,\alpha_{-x}=\alpha_{-y}=\varepsilon$ and $n=2$.

The surface plots have been computed using Eq. (25).

6.6 Figure 6b-c

The parameters for Figure 6b are as follows: The dissociation ratio is fixed at 1/2000 for all binding/unbinding reactions and we let the dimerization ratio be 1/90. The production ratios are $k_{x2}/k_{-x}=k_{y1}/k_{-y}=2700$. We assume that the inhibition and activation actions allow for leaks and hence we let $k_{x0}/k_{-x}=k_{y0}/k_{-y}=1080$, $k_{x3}/k_{-x}=k_{y3}/k_{-y}=675$ and $k_{x2}/k_{-x}=k_{y1}/k_{-y}=20$.

The parameters for Figure 6c are as follows: The dissociation ratios are $\alpha_{x0}/\alpha_{-x0}=10,\alpha_{x1}/\alpha_{-x1}=1/2700,\alpha_{y0}/\alpha_{-y0}=10,\alpha_{y1}/\alpha_{-x1}=1/2700$. The maximal production ratio is $k_{x1}/k_{-x}=k_{y1}/k_{-y}=2700$, while $k_{x0}/k_{-x}27,k_{x2}/k_{-y}=0$, $k_{y2}/k_{-y}=0,k_{y0}/k_{-y}=270$.

We ratio of the production to the leak were chosen as $k_{x1}:k_{x0}\approx 100:1$, and $k_{y1}:k_{y0}\approx 5-10:1$ [89],[90].