# Multi-modality in gene regulatory networks with slow promoter kinetics

M. Ali Al-Radhawi  Department of Mechanical Engineering, Massachusetts
Institute of Technology, Cambridge, MA 02139-4307, USA.  Domitilla Del Vecchio
Department of Mechanical Engineering, Massachusetts Institute of Technology,
Cambridge, MA 02139-4307, USA.  Eduardo D. Sontag  Departments of
Bioengineering and of Electrical and Computer Engineering, Northeastern
University, Boston, MA 02115, USA. Email<EMAIL_ADDRESS>Harvard Program
in Therapeutic Science, Harvard Medical School, Boston, MA 02115, USA.

(This manuscript was compiled on )

###### 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.

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.

Figure 1: Emergence of multi-modality due to slow promoter kinetics. (a) A
diagram of a self-repressing gene. (b) The stationary probability distribution
for different $\varepsilon$ which shows the transition from fast promoter
kinetics, i.e., $\varepsilon\to\infty$, to slow promoter kinetics, i.e.,
$\varepsilon\to 0$, in a non-cooperative self-repressing gene. The stationary
distribution is bimodal for small $\varepsilon$ and unimodal for large
$\varepsilon$. The deterministic equilibrium coincides with the fast kinetics
mode, refer to SI-§5.1. The slow kinetic limit is calculated via (25), while
the fast kinetics limit is calculated by iterating the recurrence relation in
Proposition SI-2. The remaining curves are computed by a finite projection
solution [30] of the master equation. (c) A diagram of a repression-activation
two-node network. (d) Slow promoter kinetics gives rise to four modes while
the deterministic model admits a unique stable equilibrium, refer to SI-§5.2.
The surface is plotted using (25).

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.

Figure 2: A Gene Expression Block. A generic regulatory network that consists
of gene expression blocks. A block consists of a gene reactions block and a
protein reactions block. The gene reactions are described in the text. TF is a
vector of TFs which can be monomers, dimers, or higher order multimers.
$\mathrm{D}^{i}$ is a vector whose components consist of the
$\mathrm{D}_{j}^{i}$’s. The dimension of TF is equal to the number of binding
sites of the gene.

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:

$\mathrm{D}_{0}^{i}\xrightleftharpoons[\alpha_{-i}]{\alpha_{i}}\mathrm{D}_{1}^{i}.$
(1)

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:

${\rm
TF}+\mathrm{D}_{0}^{i}\xrightleftharpoons[\alpha_{-i}]{\alpha_{i}}\mathrm{D}_{1}^{i},$
(2)

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:

$\displaystyle{\rm TF}_{1}+\mathrm{D}_{00}^{i}$
$\displaystyle\xrightleftharpoons[\alpha_{-i1}]{\alpha_{i1}}\mathrm{D}_{10}^{i}$
(3) $\displaystyle{\rm TF}_{1}+\mathrm{D}_{01}^{i}$
$\displaystyle\xrightleftharpoons[\alpha_{-i2}]{\alpha_{i2}}\mathrm{D}_{11}^{i},$
(4) $\displaystyle{\rm TF}_{2}+\mathrm{D}_{00}^{i}$
$\displaystyle\xrightleftharpoons[\alpha_{-i3}]{\alpha_{i3}}\mathrm{D}_{01}^{i},$
(5) $\displaystyle{\rm TF}_{2}+\mathrm{D}_{10}^{i}$
$\displaystyle\xrightleftharpoons[\alpha_{-i4}]{\alpha_{i4}}\mathrm{D}_{11}^{i}.$
(6)

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:

$\mathrm{D}_{j}^{i}\xrightarrow{k_{ij}}\mathrm{D}_{j}^{i}+\mathrm{X}_{i},$ (7)

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:

$\mathrm{X}_{i}\xrightarrow{k_{-i}}\emptyset.$ (8)

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:

$n\mathrm{X}_{i}\xrightleftharpoons[\beta_{-i}]{\beta_{i}}\mathrm{X}_{ic}.$
(9)

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]:

$R_{j}(z_{i})=k_{j}\binom{z_{i}}{\alpha_{ij}}=k_{j}\frac{z_{i}(z_{i}-1)...(z_{i}-\alpha_{ij})}{\alpha_{ij}!},$

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:

$R_{j}(z_{i_{1}},z_{i_{2}})=k_{j}z_{i_{1}}z_{i_{2}}.$

##### 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:

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

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$:

$p_{d}(t):=[p_{dx_{0}}(t),p_{dx_{1}}(t),....]^{T},$ (11)

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

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

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:

$\bar{\pi}=\lim_{t\to\infty}p(t).$ (13)

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

$\displaystyle\dot{p}(t)$
$\displaystyle=\Lambda_{\varepsilon}p(t)=\left(\tilde{\Lambda}+\varepsilon\hat{\Lambda}\right)p(t),$
(14)

where

$\displaystyle\tilde{\Lambda}=\begin{bmatrix}\Lambda_{0}&&\\\ &\ddots&\\\
&&\Lambda_{L-1}\end{bmatrix},\ \text{and}\ p(t)=\begin{bmatrix}p_{0}(t)\\\
\vdots\\\ p_{L-1}(t)\end{bmatrix},$ (15)

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:

${\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},}$
(16)

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:

$\emptyset\xrightleftharpoons[k_{-i}]{k_{id_{i}}}\mathrm{X}_{i},$ (17)

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:

$\emptyset\xrightleftharpoons[k_{-i}]{k_{id_{i}}}\mathrm{X}_{i},\
n_{i}\mathrm{X}_{i}\xrightleftharpoons[\beta_{-i}]{\beta_{i}}\mathrm{X}_{ic}.$
(18)

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

$\pi^{(J)}_{X|d}(x)=\prod_{i=1}^{N}\pi_{X|di}(x_{i}),$ (19)

where

$\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.,$
(20)

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:

$\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_{id_{i}}}{k_{-i}}\right),$
(21)

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

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

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:

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

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

$\Lambda_{r}:=\begin{bmatrix}\mathbf{1}^{T}&\mathbf{0}^{T}&...&\mathbf{0}^{T}\\\
\mathbf{0}^{T}&\mathbf{1}^{T}&...&\mathbf{0}^{T}\\\ &&\ddots&\\\
\mathbf{0}^{T}&\mathbf{0}^{T}&...&\mathbf{1}^{T}\end{bmatrix}\hat{\Lambda}\,[\bar{\pi}_{X|0}\
\bar{\pi}_{X|1}\ ...\ \bar{\pi}_{X|L-1}],$ (23)

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:

$\pi^{(J)}(x):=\lim_{\varepsilon\to
0^{+}}\lim_{t\to\infty}\Pr[X(t)=x]=\sum_{d=0}^{L-1}\lambda_{d}\pi_{X|d}^{(J)}(x),$
(24)

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:

$\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_{id_{i}}}{k_{-i}}\right).$
(25)

######  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:

$\mathrm{TF}+\mathrm{D}_{d_{i}}^{i}\xrightarrow{\alpha}\mathrm{D}_{d_{i^{\prime}}}^{i},$

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:

$\mathrm{D}_{d_{i}}^{i}\xrightarrow{\alpha\mathbb{E}[\text{TF}|D=d]}\mathrm{D}_{d_{i^{\prime}}}^{i}.$
(26)

Using Proposition 2, we can write:

$\mathbb{E}[\text{TF}|D=d]={\frac{\alpha}{n_{i}!}\frac{\beta_{i}}{\beta_{-i}}\left(\frac{k_{id_{i}}}{k_{-i}}\right)^{n_{i}}}.$
(27)

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:

$\displaystyle\mathrm{D}_{0}$
$\displaystyle\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}$
$\displaystyle\mathrm{D}_{1}$
$\displaystyle\mathop{\rightarrow}\limits^{k}\mathrm{X}+\mathrm{D}_{1},$ (28)
$\displaystyle\mathrm{X}$ $\displaystyle\mathop{\rightarrow}\limits^{k_{-}}0.$

Referring to Figure 2, we identify a single gene block with two states. Using
(19), the conditional stationary distributions are two Poissions at $0$ 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)

$\pi(x)=\frac{\alpha_{-}}{\alpha+\alpha_{-}}\mathbf{P}(x;0)+\frac{\alpha}{\alpha+\alpha_{-}}\mathbf{P}(x;k/k_{-}),$

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:

$\displaystyle\mathrm{X}+\mathrm{D}_{0}$
$\displaystyle\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}$
$\displaystyle\mathrm{D}_{0}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{0}}\mathrm{D}_{0}+\mathrm{X},$
(29) $\displaystyle\mathrm{D}_{1}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{1}}\mathrm{D}_{1}+\mathrm{X},$
$\displaystyle\mathrm{X}$ $\displaystyle\mathop{\rightarrow}\limits^{k_{-}}0.$

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:

$\pi_{1}(x)=\frac{\alpha\rho_{1}}{\alpha_{-}+\alpha\rho}\mathbf{P}(x;k_{1}/k_{-})+\frac{\alpha_{-}}{\alpha_{-}+\alpha\rho}\mathbf{P}(x;k_{0}/k_{-}),$
(30)

where

$\rho=\mathbb{E}[X_{2}|D=0]=k_{0}/k_{-}.$ (31)

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

$\displaystyle\mathrm{X}_{2}+\mathrm{D}_{0}$
$\displaystyle\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}$
$\displaystyle\mathrm{D}_{0}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{0}}D_{0}+\mathrm{X},$ (32)
$\displaystyle\mathrm{D}_{1}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{1}}\mathrm{D}_{1}+\mathrm{X},$
$\displaystyle\mathrm{X}$ $\displaystyle\mathop{\rightarrow}\limits^{k_{-}}0$
$\displaystyle 2\mathrm{X}$
$\displaystyle\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{X}_{2}.$

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

$\rho_{2}=\mathbb{E}[X_{2}|D=0]=\frac{k_{0}^{2}\beta}{2k_{-}^{2}\beta_{-}}.$
(33)

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.

| 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
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.

Figure 3: More modes emerge due to slow promoter kinetics. (a) A self-
activating gene. (b) The stationary probability distribution for different
$\varepsilon$ which shows the transition from fast promoter kinetics to slow
promoter kinetics in a leaky cooperative self-activation of a gene given by (A
Self-Regulating Gene). The slow kinetic limit is calculated via (25), while
the remaining curves are computed by a finite projection solution [30] of the
CME.  The details and parameters are given in SI-§5.4 .

### 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$:

$\displaystyle\begin{array}[]{rl}\mathrm{Y}_{m}+\mathrm{D}_{0}^{X}&\xrightleftharpoons[\alpha_{-1}]{\alpha_{1}}\mathrm{D}_{1}^{X}\\\
\mathrm{D}_{0}&\mathop{\rightarrow}\limits^{k_{10}}\mathrm{D}_{0}^{X}+\mathrm{X},\\\
\mathrm{X}&\mathop{\rightarrow}\limits^{k_{-1}}0\\\
n\mathrm{X}&\xrightleftharpoons[\beta_{-1}]{\beta_{1}}\mathrm{X}_{n}\\\
\end{array}\begin{array}[]{rl}\mathrm{X}_{n}+\mathrm{D}_{0}^{Y}&\xrightleftharpoons[\alpha_{-2}]{\alpha_{2}}\mathrm{D}_{1}^{Y}\\\
\mathrm{D}_{0}^{Y}&\mathop{\rightarrow}\limits^{k_{20}}D_{0}^{Y}+\mathrm{Y},\\\
\mathrm{Y}&\mathop{\rightarrow}\limits^{k_{-2}}0\\\
mX&\xrightleftharpoons[\beta_{-2}]{\beta_{2}}\mathrm{Y}_{m}.\end{array}$

$\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:

$\pi(x,y)=\frac{1}{\frac{\alpha_{1}}{\alpha_{-1}}\rho_{2}+\frac{\alpha_{2}}{\alpha_{-2}}\rho_{1}+1}\left(\mathbf{P}(y;{\tfrac{k_{20}}{k_{-2}}})\mathbf{P}(x;{\tfrac{k_{10}}{k_{-1}}})+\frac{\alpha_{1}}{\alpha_{-1}}\rho_{2}\mathbf{P}(y;{\tfrac{k_{20}}{k_{-2}}})\delta(x)+\frac{\alpha_{2}}{\alpha_{-2}}\rho_{1}\mathbf{P}(x;{\tfrac{k_{10}}{k_{-1}}})\delta(y)\right).$
(34)

where

$\rho_{1}=\left(\frac{k_{10}}{k_{-1}}\right)^{n}\frac{\beta_{1}}{n!\beta_{-1}},\quad\rho_{2}=\left(\frac{k_{20}}{k_{-2}}\right)^{m}\frac{\beta_{2}}{m!\beta_{-2}}.$
(35)

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.

|
---|---
|

(a) Diagram of the toggle switch. (b) Non-cooperative,  $n,m=1$.

(c) Cooperative, $\beta_{1}/\beta_{-1}=\beta_{2}/\beta_{-2}=1$. (d)
Cooperative, $\beta_{1}/\beta_{-1}=\beta_{2}/\beta_{-2}=0.01$.

Figure 4: Cooperativity enables tuning of modes’ weights. Comparison of the
stationary distribution between non-cooperative and cooperative binding. For
all cases:
$\alpha_{1}/\alpha_{-1}=\alpha_{2}/\alpha_{-2}=1/200,k_{10}/k_{-1}=k_{20}/k_{-2}=40$.
(a) Diagram of the toggle switch. (b) The stationary distribution
corresponding to the non-cooperative case. (c) The stationary distribution
corresponding to the cooperative case with
$m=n=2,\beta_{1}/\beta_{-1}=\beta_{2}/\beta_{-2}=1$. (d) The stationary
distribution corresponding to the cooperative case with
$m=n=2,\beta_{1}/\beta_{-1}=\beta_{2}/\beta_{-2}=0.01$. All surfaces are
plotted using (25).

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:

$\begin{array}[]{rcl}\mathrm{Y}_{ic}+\mathrm{D}_{0}^{xi}&\xrightleftharpoons[\varepsilon\alpha_{-x}]{\varepsilon\alpha_{x}}&\mathrm{D}_{1}^{xi}\\\
\mathrm{D}_{0}^{xi}&\mathop{\rightarrow}\limits^{k_{x}}&\mathrm{X}_{i}+\mathrm{D}_{0}^{xi},\\\
\mathrm{X}_{i}&\mathop{\rightarrow}\limits^{k_{-x}}&0,\\\
n\mathrm{X}_{i}&\xrightleftharpoons[\beta_{-x}]{\beta_{x}}&\mathrm{X}_{ic},\end{array}\begin{array}[]{rcl}\mathrm{X}_{ic}+\mathrm{D}_{0}^{yi}&\xrightleftharpoons[\varepsilon\alpha_{-y}]{\varepsilon\alpha_{y}}&\mathrm{D}_{1}^{yi}\\\
\mathrm{D}_{0}^{yi}&\mathop{\rightarrow}\limits^{k_{y}}&\mathrm{Y}_{i}+\mathrm{D}_{0}^{yi},\\\
\mathrm{Y}_{i}&\mathop{\rightarrow}\limits^{k_{-y}}&0,\\\
n\mathrm{Y}_{i}&\xrightleftharpoons[\beta_{-y}]{\beta_{y}}&\mathrm{Y}_{ic},\end{array}$

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$:

$\mathrm{X}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{X}_{j},\quad\mathrm{Y}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{Y}_{j},\quad
i\neq j,\;i,j=1,..,N.$ (36)

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:

$\left\\{\left(\frac{ik_{x}}{Nk_{-x}},0\right),\left(0,\frac{ik_{y}}{Nk_{-y}}\right):i=1,..,N\right\\}.$

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:

$\Omega\geq\frac{1}{N}\max\\{k_{x}-k_{-x},k_{y}-k_{-y}\\}.$ (37)

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.

(a) Network of toggle switches Figure 5: Slow promoter kinetics lead to the
emergence of a multi- modal toggle switch (a) A diagram of population of
toggle switches. Arrows between blocks represent reversible diffusion
reactions (36). Each block contains a toggle switch. The remaining subfigures
show stationary probability distributions for a population of three identical
cooperative toggle switches. Due to the symmetries we plot joint distributions
of $X_{1},Y_{1}$ and $X_{1},X_{2}$ only. Subplots (b),(c) depict the uncoupled
toggle switches. It can be noted that the $X_{1}$ and $X_{2}$ are not
synchronized. Subplots (d), (e) depict a high diffusion case. The toggle
switches synchronize into a multi-modal toggle switch. Parameters and details
for calculating the plots are given in SI-§5.5

### 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]:

$\begin{array}[]{rl}\mathrm{X}_{2}+\mathrm{D}_{00}^{X}\xrightleftharpoons[]{}\mathrm{D}_{10}^{X},&\mathrm{X}_{2}+\mathrm{D}_{01}^{X}\xrightleftharpoons[]{}\mathrm{D}_{11}^{X},\\\
\mathrm{Y}_{2}+\mathrm{D}_{00}^{X}\xrightleftharpoons[]{}\mathrm{D}_{01}^{X},&\mathrm{Y}_{2}+\mathrm{D}_{10}^{X}\xrightleftharpoons[]{}\mathrm{D}_{11}^{X},\\\
\mathrm{X}_{2}+\mathrm{D}_{00}^{Y}\xrightleftharpoons[]{}\mathrm{\mathrm{missing}}D_{10}^{Y},&\mathrm{X}_{2}+\mathrm{D}_{01}^{Y}\xrightleftharpoons[]{}\mathrm{D}_{11}^{Y},\\\
\mathrm{Y}_{2}+\mathrm{D}_{00}^{Y}\xrightleftharpoons[]{}\mathrm{D}_{01}^{Y},&\mathrm{Y}_{2}+\mathrm{D}_{10}^{Y}\xrightleftharpoons[]{}\mathrm{D}_{11}^{Y},\\\
2\mathrm{X}\xrightleftharpoons[\beta_{-x}]{\beta_{x}}\mathrm{X}_{2},&2\mathrm{Y}\xrightleftharpoons[\beta_{-y}]{\beta_{y}}\mathrm{Y}_{2}\end{array}\begin{array}[]{rl}\mathrm{D}_{00}^{X}\mathop{\rightarrow}\limits^{k_{x0}}\mathrm{D}_{00}^{X}+\mathrm{X},&\mathrm{D}_{01}^{X}\mathop{\rightarrow}\limits^{k_{x1}}\mathrm{D}_{01}^{X}+\mathrm{X},\\\
\mathrm{D}_{10}^{X}\mathop{\rightarrow}\limits^{k_{x2}}\mathrm{D}_{10}^{X}+\mathrm{X},&\mathrm{D}_{11}^{X}\mathop{\rightarrow}\limits^{k_{x3}}\mathrm{D}_{11}^{X}\mathop{\rightarrow}\limits^{k_{x3}}\mathrm{D}_{11}^{X}+\mathrm{X},\\\
\mathrm{D}_{00}^{Y}\mathop{\rightarrow}\limits^{k_{y0}}\mathrm{D}_{00}^{Y}+\mathrm{Y},&\mathrm{D}_{01}^{Y}\mathop{\rightarrow}\limits^{k_{y1}}\mathrm{D}_{01}^{Y}+\mathrm{Y},\\\
\mathrm{D}_{10}^{Y}\mathop{\rightarrow}\limits^{k_{y2}}\mathrm{D}_{10}^{Y}+\mathrm{Y},&\mathrm{D}_{11}^{Y}\mathop{\rightarrow}\limits^{k_{y3}}\mathrm{D}_{11}^{Y}\mathop{\rightarrow}\limits^{k_{y3}}\mathrm{D}_{11}^{Y}+\mathrm{Y},\\\
\mathrm{X}\mathop{\rightarrow}\limits^{k_{-x}}\emptyset,&\mathrm{Y}\mathop{\rightarrow}\limits^{k_{-y}}\emptyset.\\\
\end{array}$ (38)

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}$.

(a) A cell fate network (b) The stationary distribution of cell-fate circuit
(38)  (c) The stationary distribution of a PU.1/GATA.1 network (39)

Figure 6: The cell-fate decision network with slow promoter kinetics has more
modes than what a deterministic model predicts. (a) A diagram of a generic
cell-fate circuit that can describe both networks (38), (39), (b) The joint
probability of an archetypical cell-fate circuit (38) computed using Theorem
3. (c) The stationary distribution of a PU.1/GATA.1 circuit (39), where
$\mathrm{X}$ denotes PU.1 and $\mathrm{Y}$ denotes GATA.1. Three modes can be
seen. The parameters are given in SI-§5.6.

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:

$\begin{array}[]{rl}\mathrm{X}_{2}+\mathrm{D}_{00}^{X}\xrightleftharpoons[\alpha_{-x0}]{\alpha_{x0}}\mathrm{D}_{10}^{X},&\mathrm{Y}_{2}+\mathrm{D}_{10}^{X}\xrightleftharpoons[\alpha_{-x1}]{\alpha_{x1}}\mathrm{D}_{11}^{X},\\\
\mathrm{Y}_{2}+\mathrm{D}_{00}^{Y}\xrightleftharpoons[\alpha_{-y0}]{\alpha_{y0}}\mathrm{D}_{01}^{Y},&\mathrm{X}_{2}+\mathrm{D}_{01}^{Y}\xrightleftharpoons[\alpha_{-y1}]{\alpha_{y1}}\mathrm{D}_{11}^{Y},\\\
\mathrm{X}\mathop{\rightarrow}\limits^{k_{-x}}\emptyset,&\mathrm{Y}\mathop{\rightarrow}\limits^{k_{-y}}\emptyset,\\\
\end{array}\begin{array}[]{rl}\mathrm{D}_{00}^{X}\mathop{\rightarrow}\limits^{k_{x0}}\mathrm{D}_{00}^{X}+\mathrm{X},&\mathrm{D}_{10}^{X}\mathop{\rightarrow}\limits^{k_{x1}}\mathrm{D}_{10}^{X}+\mathrm{X},\\\
\mathrm{D}_{00}^{Y}\mathop{\rightarrow}\limits^{k_{y0}}\mathrm{D}_{00}^{Y}+\mathrm{Y}&\mathrm{D}_{01}^{Y}\mathop{\rightarrow}\limits^{k_{y1}}\mathrm{D}_{01}^{Y}+\mathrm{Y},\\\
\end{array}$ (39)

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:

$\mathrm{R}_{j}:\
\sum_{i=1}^{|\mathscr{S}|}\alpha_{ij}\mathrm{Z}_{i}\rightarrow\sum_{i=1}^{|\mathscr{S}|}\beta_{ij}\mathrm{Z}_{i},$
(40)

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

$\mathrm{R}_{j}:\
\sum_{i=1}^{|\mathscr{S}|}\mathrm{\alpha}_{ij}\mathrm{Z}_{i}\rightleftharpoons\sum_{i=1}^{|\mathscr{S}|}\beta_{ij}\mathrm{Z}_{i}.$
(41)

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

$[\Gamma]_{ij}=\beta_{ij}-\alpha_{ij}.$

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:

$\dot{p}_{z}(t)=\sum_{j=1}^{|\mathscr{R}|}R_{j}(z-\gamma_{j})p_{z-\gamma_{j}}(t)-R_{j}(z)p_{z}(t),{z\in\mathbf{Z}},$
(42)

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:

$\sum_{j\in B_{i}}D_{j}^{i}(t)=1,$

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:

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

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$:

$p_{d}(t):=[p_{dx_{0}}(t),p_{dx_{1}}(t),....]^{T},$ (44)

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

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

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:

$\bar{\pi}=\lim_{t\to\infty}p(t).$ (46)

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:

$\lambda_{z\tilde{z}}:=\left\\{\begin{array}[]{ll}R_{j}(z)&\mbox{if}\,\exists
j\ \mbox{such that}\ \tilde{z}=z-\gamma_{j}\\\ -\sum_{\tilde{z}\neq
z}\lambda_{z\tilde{z}}=-\sum_{j=1}^{|\mathscr{R}|}R_{j}(z)&\mbox{if}\
\tilde{z}=z\\\ 0&\mbox{otherwise}\end{array}\right.,$ (47)

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 $0$. 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 $0$ does not
lead to any state then $\\{0\\}$ is a closed communicating class. Otherwise,
let $U$ be the smallest communicating class containing $0$. 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:

$\lambda_{z\tilde{z}}:=\left\\{\begin{array}[]{ll}R_{j}^{(f)}(z)+\varepsilon
R_{j}^{(s)}(z)&\mbox{if}\,\exists j\ \mbox{such that}\
\tilde{z}=z-\gamma_{j}\\\ -\sum_{\tilde{z}\neq
z}\lambda_{z\tilde{z}}=-\sum_{j=1}^{|\mathscr{R}|}R_{j}^{(f)}(z)+\varepsilon
R_{j}^{(s)}(z)&\mbox{if}\ \tilde{z}=z\\\
0&\mbox{otherwise}\end{array}\right.,$ (48)

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

$\displaystyle\dot{p}_{x_{1},x_{2}}$
$\displaystyle=\left(kp_{x_{1}-1,x_{2}}-k_{-}x_{1}p_{x_{1},x_{2}}\right)+\left(k_{-}(x_{1}+1)p_{x_{1}+1,x_{2}}-kp_{x_{1},x_{2}}\right)$
(49)
$\displaystyle+\left(\beta_{-2}(x_{2}+1)p_{x_{1}-2,x_{2}+1}-\tfrac{1}{n!}\beta\prod_{k=0}^{n-1}(x_{1}-k)p_{x_{1},x_{2}}\right)$
(50)
$\displaystyle+\left(\tfrac{1}{n!}\beta\prod_{k=1}^{n}(x_{1}+k)p_{x_{1}+n,x_{2}-1}-\beta_{-}x_{2}p_{x_{1},x_{2}}\right)..$
(51)

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

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

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:

$\displaystyle{\tilde{\Lambda}\bar{\pi}^{(0)}}$ $\displaystyle=0$ (53)
$\displaystyle{\tilde{\Lambda}\bar{\pi}^{(1)}+\hat{\Lambda}\bar{\pi}^{(0)}}$
$\displaystyle=0$ (54)

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:

$\bar{\pi}_{X|d}:=[\overbrace{\mathbf{0}^{T}\ ...\
\mathbf{0}^{T}}^{d-1}\pi_{X|d}^{T}\overbrace{\mathbf{0}^{T}\ ...\
\mathbf{0}^{T}}^{L-d}]^{T}.$ (55)

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:

$\bar{\pi}^{(0)}=\sum_{i=0}^{L-1}\lambda_{i}\bar{\pi}_{X|i},$

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:

$\Lambda_{r}\lambda:=\begin{bmatrix}\mathbf{1}^{T}&\mathbf{0}^{T}&...&\mathbf{0}^{T}\\\
\mathbf{0}^{T}&\mathbf{1}^{T}&...&\mathbf{0}^{T}\\\ &&\ddots&\\\
\mathbf{0}^{T}&\mathbf{0}^{T}&...&\mathbf{1}^{T}\end{bmatrix}\hat{\Lambda}\,[\bar{\pi}_{X|0}\
\bar{\pi}_{X|1}\ ...\ \bar{\pi}_{X|L-1}]\begin{bmatrix}\lambda_{0}\\\
\vdots\\\ \lambda_{L-1}\end{bmatrix}=0.$ (56)

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

$\lambda_{0}+...+\lambda_{L-1}=1.$ (57)

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:

$\Lambda_{r}=\begin{bmatrix}\mathbf{1}^{T}\hat{\Lambda}_{00}\pi_{0}&\dots&\mathbf{1}^{T}\hat{\Lambda}_{0{L-1}}\pi_{L-1}\\\
\vdots&\ddots&\\\
\mathbf{1}^{T}\hat{\Lambda}_{{L-1}0}\pi_{0}&\dots&\mathbf{1}^{T}\hat{\Lambda}_{(L-1)(L-1)}\pi_{L-1}\end{bmatrix}.$
(58)

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}|\\}$:

      1. 1.

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:

$[\Lambda_{r}]_{d^{\prime}d}=\left\\{\begin{array}[]{ll}\alpha,&\mbox{if the
reaction is monomolecular}\\\ \\\
{\frac{\alpha}{n_{\bar{i}}!}\frac{\beta_{\bar{i}}}{\beta_{-\bar{i}}}\left(\frac{k_{\bar{i}d_{\bar{i}}}}{k_{-\bar{i}}}\right)^{n_{\bar{i}}}},&\mbox{if
the reaction is bimolecular}\end{array}\right.$ (59)

    * -

Set

$\mathbf{[}\Lambda_{r}]_{dd}=-\sum_{i^{\prime}\neq
i}\mathbf{1}^{T}\hat{\Lambda}_{d_{i^{\prime}}d_{i}}\pi_{d}.$ (60)

  * •

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:

$\mathrm{D}_{d_{i}}^{i}\mathop{\to}\limits^{\alpha}\mathrm{D}_{d_{i^{\prime}}},\mathrm{D}_{d_{i}}^{i}\mathop{\to}\limits^{\alpha}\mathop{\to}\limits^{\alpha}\mathrm{D}_{d_{i^{\prime}}}+\mathrm{X}_{\bar{i}},$

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

$[\Lambda_{r}]_{d^{\prime}d}=\mathbf{1}^{T}\Lambda_{d^{\prime}d}\pi_{d}=\alpha\mathbf{1}\pi_{d}=\alpha.$

Now consider a reaction of the form:

$\mathrm{X}_{\bar{i}}+\mathrm{D}_{d_{i}}^{i}\mathop{\to}\limits^{\alpha}\mathrm{D}_{d_{i^{\prime}}}^{i}.$

Then,

$[\Lambda_{r}]_{d^{\prime}d}=\mathbf{1}^{T}\Lambda_{d^{\prime}d}\pi_{d}=\alpha\sum_{x_{\bar{i}}=0}^{\infty}x_{\bar{i}}\pi_{d\bar{i}}\sum_{x_{i\neq\bar{i}}}\pi_{di}=\alpha\mathbb{E}[X_{\bar{i}}(t)|D(t)=d]=\alpha\frac{k_{\bar{i}d_{\bar{i}}}^{n_{\bar{i}}}\beta_{\bar{i}}}{n_{\bar{i}}!k_{-\bar{i}d_{\bar{i}}}^{n_{\bar{i}}}\beta_{-\bar{i}}}.$

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:

$\displaystyle\mathrm{D}_{0}$
$\displaystyle\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}$
$\displaystyle\mathrm{D}_{1}$
$\displaystyle\mathop{\rightarrow}\limits^{k}\mathrm{X}+\mathrm{D}_{1},$ (61)
$\displaystyle\mathrm{X}$ $\displaystyle\mathop{\rightarrow}\limits^{k_{-}}0.$

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

$\displaystyle\pi_{0}(x)$ $\displaystyle=\delta(x)$ $\displaystyle\pi_{1}(x)$
$\displaystyle=\mathbf{P}(x;k/k_{-}).$

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):

$\Lambda_{r}=\begin{bmatrix}-\alpha&\alpha_{-}\\\
\alpha&-\alpha_{-}\end{bmatrix}.$

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:

$\pi(x)=\frac{\alpha_{-}}{\alpha+\alpha_{-}}\mathbf{P}(x;0)+\frac{\alpha}{\alpha+\alpha_{-}}\mathbf{P}(x;k/k_{-}),$
(62)

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).

Figure 7: The stationary probability distribution for different $\varepsilon$
which shows the transition from fast promoter kinetics, i.e.,
$\varepsilon\to\infty$, to slow promoter kinetics, i.e., $\varepsilon\to 0$,
in a single unregulated gene. The stationary distribution is bimodal for small
$\varepsilon$, i.e. $\varepsilon\leq 1$, and unimodal for large $\varepsilon$.
The deterministic equilibrium coincides with the fast kinetics mode at
$\frac{\alpha}{\alpha+\alpha_{-}}\frac{k}{k_{-}}$. The slow kinetic limit is
calculated via (62), the fast kinetics limit is a Poisson centered at the
deterministic equilibrium, while the remaining curves are computed by
evaluating the exact solution given in [80] . The parameters are
$\alpha_{-}=0.1,\alpha=1,k_{-}=2,k=20$.

###  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:

$\mathrm{D}_{0}\xrightleftharpoons[\alpha_{-}]{\alpha
k_{0}/k_{-}}\mathrm{D}_{1}.$

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

$\Lambda_{r}=\begin{bmatrix}-\alpha k_{0}/k_{-}&\alpha_{-}\\\ \alpha
k_{0}/k_{-}&-\alpha_{-}\end{bmatrix}.$

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:

$\pi_{1}(x)=\frac{\alpha\rho_{1}}{\alpha_{-}+\alpha\rho_{1}}\mathbf{P}(x;k_{1}/k_{-})+\frac{\alpha_{-}}{\alpha_{-}+\alpha\rho_{1}}\mathbf{P}(x;k_{0}/k_{-}),$
(63)

where

$\rho_{1}=\mathbb{E}[X_{2}|D=0]=k_{0}/k_{-}.$

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

$\displaystyle\mathrm{X}_{2}+\mathrm{D}_{0}$
$\displaystyle\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}$
$\displaystyle\mathrm{D}_{0}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{0}}D_{0}+\mathrm{X},$ (64)
$\displaystyle\mathrm{D}_{1}$
$\displaystyle\mathop{\rightarrow}\limits^{k_{1}}\mathrm{D}_{1}+\mathrm{X},$
$\displaystyle\mathrm{X}$ $\displaystyle\mathop{\rightarrow}\limits^{k_{-}}0$
$\displaystyle 2\mathrm{X}$
$\displaystyle\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{X}_{2}.$

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:

$\pi_{2}(x)=\frac{\alpha\rho_{2}}{\alpha_{-}+\alpha\rho_{2}}\mathbf{P}(x;k_{1}/k_{-})+\frac{\alpha_{-}}{\alpha_{-}+\alpha\rho_{2}}\mathbf{P}(x;k_{0}/k_{-}),$
(65)

where

$\rho_{2}=\mathbb{E}[X_{2}|D=0]=\frac{k_{0}^{2}\beta}{2k_{-}^{2}\beta_{-}}.$

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:

$\mathrm{X}+\mathrm{D}_{0}\xrightleftharpoons[\varepsilon\alpha_{-}]{\varepsilon\alpha}\mathrm{D}_{1}+\mathrm{X}.$

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:

  1. 1.

_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:

$\pi(m)=\lim_{t\to\infty}\Pr[X(t)=m]=\frac{\alpha_{-}w_{m}+\alpha
mw_{m+1}}{\alpha m+\alpha_{-}},$

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

$\displaystyle w_{m+1}$
$\displaystyle=\frac{((k_{1}m+(\alpha_{-}/\alpha)k_{0})(m+(\alpha_{-}/\alpha)+1)}{k_{-}(m+1)(m+(\alpha_{-}/\alpha))^{2}}w_{m},\
m\geq 0.$

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:

$\Pr[X+D_{1}=m]=w_{m},$

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

$\displaystyle w_{m+1}$
$\displaystyle=\frac{((k_{1}m+(\alpha_{-}/\alpha)k_{0})(m+(\alpha_{-}/\alpha)+1)}{k_{-}(m+1)(m+(\alpha_{-}/\alpha))^{2}}w_{m},\
m\geq 0.$

The joint distribution can be given as:

$\displaystyle\Pr[X+D_{1}=m,D_{1}=0]$
$\displaystyle=w_{m}\frac{\alpha_{-}}{\alpha m+\alpha_{-}},$
$\displaystyle\Pr[X+D_{1}=m,D_{1}=1]$ $\displaystyle=w_{m}\frac{\alpha
m}{\alpha m+\alpha_{-}}.$

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

$\displaystyle\Pr[X=m]$ $\displaystyle=\Pr[X+D=m,D=0]+\Pr[X+D=m+1,D=1]$
$\displaystyle=\frac{\alpha_{-}w_{m}+\alpha mw_{m+1}}{\alpha m+\alpha_{-}}.$

∎

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$:

$\displaystyle\begin{array}[]{rl}\mathrm{Y}_{m}+\mathrm{D}_{0}^{X}&\xrightleftharpoons[\alpha_{-1}]{\alpha_{1}}\mathrm{D}_{1}^{X}\\\
\mathrm{D}_{0}&\mathop{\rightarrow}\limits^{k_{10}}\mathrm{D}_{0}^{X}+\mathrm{X},\\\
\mathrm{X}&\mathop{\rightarrow}\limits^{k_{-1}}0\\\
n\mathrm{X}&\xrightleftharpoons[\beta_{-1}]{\beta_{1}}\mathrm{X}_{n}\\\
\end{array}\begin{array}[]{rl}\mathrm{X}_{n}+\mathrm{D}_{0}^{Y}&\xrightleftharpoons[\alpha_{-2}]{\alpha_{2}}\mathrm{D}_{1}^{Y}\\\
\mathrm{D}_{0}^{Y}&\mathop{\rightarrow}\limits^{k_{20}}D_{0}^{Y}+\mathrm{Y},\\\
\mathrm{Y}&\mathop{\rightarrow}\limits^{k_{-2}}0\\\
m\mathrm{Y}&\xrightleftharpoons[\beta_{-2}]{\beta_{2}}\mathrm{Y}_{m}.\end{array}$

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:

$\Lambda_{r}=\begin{bmatrix}-\alpha_{1}\rho_{2}-\alpha_{2}\rho_{1}&\alpha_{-2}&\alpha_{-1}&0\\\
\alpha_{2}\rho_{1}&-\alpha_{-2}&0&\alpha_{-1}\\\
\alpha_{1}\rho_{2}&0&-\alpha_{-1}&\alpha_{-2}\\\
0&0&0&-\alpha_{-2}-\alpha_{-1}\end{bmatrix},$ (66)

where

$\rho_{1}=\left(\frac{k_{10}}{k_{-1}}\right)^{n}\frac{\beta_{1}}{n!\beta_{-1}},\quad\rho_{2}=\left(\frac{k_{20}}{k_{-2}}\right)^{m}\frac{\beta_{2}}{m!\beta_{-2}}.$
(67)

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:

$\pi(x,y)=\frac{1}{\frac{\alpha_{1}}{\alpha_{-1}}\rho_{2}+\frac{\alpha_{2}}{\alpha_{-2}}\rho_{1}+1}\left(\mathbf{P}(y;{\tfrac{k_{20}}{k_{-2}}})\mathbf{P}(x;{\tfrac{k_{10}}{k_{-1}}})+\frac{\alpha_{1}}{\alpha_{-1}}\rho_{2}\mathbf{P}(y;{\tfrac{k_{20}}{k_{-2}}})\delta(x)+\frac{\alpha_{2}}{\alpha_{-2}}\rho_{1}\mathbf{P}(x;{\tfrac{k_{10}}{k_{-1}}})\delta(y)\right).$
(68)

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:

$\begin{array}[]{llll}\mathrm{X}_{c}+\mathrm{D}_{00}\xrightleftharpoons[]{}\mathrm{D}_{10},&\mathrm{X}_{c}+\mathrm{D}_{01}\xrightleftharpoons[]{}\mathrm{D}_{11},&\mathrm{Y}_{c}+\mathrm{D}_{00}\xrightleftharpoons[]{}\mathrm{D}_{01},&\mathrm{Y}_{c}+\mathrm{D}_{10}\xrightleftharpoons[]{}\mathrm{D}_{11},\\\
\mathrm{D}_{00}\mathop{\rightarrow}\limits^{k_{x}}\mathrm{D}_{00}+\mathrm{X},&\mathrm{D}_{00}\mathop{\rightarrow}\limits^{k_{y}}\mathrm{D}_{00}+\mathrm{Y},&\mathrm{D}_{01}\mathop{\rightarrow}\limits^{k_{x}}\mathrm{D}_{01}+\mathrm{X},&\quad\mathrm{D}_{10}\mathop{\rightarrow}\limits^{k_{y}}\mathrm{D}_{10}+\mathrm{Y},\\\
\mathrm{X}\mathop{\rightarrow}\limits^{k_{-x}}\emptyset,&\mathrm{Y}\mathop{\rightarrow}\limits^{k_{-y}}\emptyset,&n\mathrm{X}\xrightleftharpoons[\beta_{-x}]{\beta_{x}}\mathrm{X}_{c},&\quad
m\mathrm{Y}\xrightleftharpoons[\beta_{-y}]{\beta_{y}}\mathrm{Y}_{c},\par\end{array}$
(69)

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.

Figure 8:  Analysis predicts the stationary distribution with 0.5 time scale
separation. (a) The empirical stationary distribution has been computed based
on a realization of the stochastic process. The realization has been computed
via Gillespie’s algorithm and has been averaged over $70\times 10^{6}$
iterations. (b) A sample trajectory of the toggle switch. (c) A zoomed plot
showing the three modes Observe that between times 18 and 18.5 both modes are
high. The parameters are:
$\alpha_{-}=10,\alpha=0.1,k_{-}=20,k=2000,\beta=20,\beta_{-}=2000$. Note that
the unbinding rate is two times slower than the protein decay rate.

###  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:

$\displaystyle\begin{array}[]{rl}\mathrm{Z}_{n}+\mathrm{D}_{0}^{X}&\xrightleftharpoons[\alpha_{-}]{\alpha}\mathrm{D}_{1}^{X}\\\
\mathrm{D}_{0}&\mathop{\rightarrow}\limits^{k}\mathrm{D}_{0}^{X}+\mathrm{X},\\\
\mathrm{X}&\mathop{\rightarrow}\limits^{k_{-}}0\\\
n\mathrm{X}&\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{X}_{n}\\\
\end{array}\begin{array}[]{rl}\mathrm{X}_{n}+\mathrm{D}_{0}^{Y}&\xrightleftharpoons[\alpha_{-}]{\alpha}\mathrm{D}_{1}^{Y}\\\
\mathrm{D}_{0}^{Y}&\mathop{\rightarrow}\limits^{k}D_{0}^{Y}+\mathrm{Y},\\\
\mathrm{Y}&\mathop{\rightarrow}\limits^{k_{-}}0\\\
n\mathrm{Y}&\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{Y}_{n}.\end{array}\begin{array}[]{rl}\mathrm{Y}_{n}+\mathrm{D}_{0}^{Z}&\xrightleftharpoons[\alpha_{-}]{\alpha}\mathrm{D}_{1}^{Z}\\\
\mathrm{D}_{0}^{Z}&\mathop{\rightarrow}\limits^{k}D_{0}^{Z}+\mathrm{Z},\\\
\mathrm{Z}&\mathop{\rightarrow}\limits^{k_{-}}0\\\
n\mathrm{Z}&\xrightleftharpoons[\beta_{-}]{\beta}\mathrm{Z}_{n}.\end{array}$

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:

$\Lambda_{r}=\begin{bmatrix}-3\alpha\rho&\alpha_{-}&\alpha_{-}&0&\alpha_{-}&0&0&0\\\
\alpha\rho&-\alpha_{-}-\alpha\rho&0&\alpha_{-}&0&\alpha_{-}&0&0\\\
\alpha\rho&0&-\alpha_{-}-\alpha\rho&\alpha_{-}&0&0&\alpha_{-}&0\\\
0&\alpha\rho&0&-2\alpha_{-}&0&0&0&\alpha_{-}\\\
\alpha\rho&0&0&0&-\alpha_{-}-\alpha\rho&\alpha_{-}&\alpha_{-}&0\\\
0&0&0&0&\alpha\rho&-2\alpha_{-}&0&\alpha_{-}\\\
0&0&\alpha\rho&0&0&0&-2\alpha_{-}&\alpha_{-}\\\
0&0&0&0&0&0&0&-3\alpha_{-}.\end{bmatrix},$ (70)

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

The stationary distribution is:

$\displaystyle\pi(x,y,z)$
$\displaystyle=\frac{1}{D}\Big{(}w\mathbf{P}(x,y,z;K,0,0)+w\mathbf{P}(x,y,z;0,K,0)+w\mathbf{P}(x,y,z;0,0,K)$
(71)
$\displaystyle\left.+\frac{1}{2w}\mathbf{P}(x,y,z;K,K,K)+\mathbf{P}(x,y,z;K,K,0)+\mathbf{P}(x,y,z;0,K,K)+\mathbf{P}(x,y,z;K,0,K)\right),$

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

$w=2\frac{\alpha}{\alpha_{-}}\left(\frac{k}{k_{-}}\right)^{n}\frac{\beta}{n!\beta_{-}}\gg
1,$

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

$Q(t)=e^{\Lambda_{d}}e^{\Lambda_{+}}.$

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:

$\displaystyle Q_{d_{x}d_{x}}$
$\displaystyle=e^{-2t}+e^{-2t}o(t)=e^{-2t}(1+wt^{2})+e^{-2t}o(t^{3})$
$\displaystyle Q_{d_{y}d_{x}}$
$\displaystyle=e^{-2t}o(t)=e^{-2t}(\tfrac{1}{2}wt^{2})+e^{-2t}o(t^{3})$
$\displaystyle Q_{d_{z}d_{x}}$ $\displaystyle=e^{-2t}o(t^{3}),$

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.

Figure 9:  A diagram of the repressilator.

(a) Cooperative, Slow-Promoter Kinetics. (b) Noncooperative, Slow-Promoter
Kinetics. (c) Cooperative, Fast-Promoter Kinetics. (d) Noncooperative, Fast-
Promoter Kinetics.

Figure 10:  The noncooperative repressilator oscillates.  (a) A time-series
for the cooperative repressilator with cooperativity index 2, and slow
promoter kinetics. (b) A time-series for the cooperative repressilator with
cooperativity index 2, and fast promoter kinetics. (c) A time-series for the
noncooperative repressilator, and slow promoter kinetics. (d) A time-series
for the noncooperative repressilator, and fast promoter kinetics. The plots
were generated by stochastic simulation via the Gillespie algorithm. For all
the figures, the parameters are:
$\alpha=5\varepsilon,\alpha_{-}=1\varepsilon,k=2000,k_{-}=20,\beta_{\pm}=1$,
where $\varepsilon=0.1$ for slow kinetics, and $\varepsilon=1000$ for fast
kinetics.

###  3.5 Synchronization of interconnected toggle switches

We consider $N$ identical toggle switches:

$\begin{array}[]{rcl}\mathrm{Y}_{ic}+\mathrm{D}_{0}^{xi}&\xrightleftharpoons[\varepsilon\alpha_{-x}]{\varepsilon\alpha_{x}}&\mathrm{D}_{1}^{xi}\\\
\mathrm{D}_{0}^{xi}&\mathop{\rightarrow}\limits^{k_{x}}&\mathrm{X}_{i}+\mathrm{D}_{0}^{xi},\\\
\mathrm{X}_{i}&\mathop{\rightarrow}\limits^{k_{-x}}&0,\\\
n\mathrm{X}_{i}&\xrightleftharpoons[\beta_{-x}]{\beta_{x}}&\mathrm{X}_{ic},\end{array}\begin{array}[]{rcl}\mathrm{X}_{ic}+\mathrm{D}_{0}^{yi}&\xrightleftharpoons[\varepsilon\alpha_{-y}]{\varepsilon\alpha_{y}}&\mathrm{D}_{1}^{yi}\\\
\mathrm{D}_{0}^{yi}&\mathop{\rightarrow}\limits^{k_{y}}&\mathrm{Y}_{i}+\mathrm{D}_{0}^{yi},\\\
\mathrm{Y}_{i}&\mathop{\rightarrow}\limits^{k_{-y}}&0,\\\
n\mathrm{Y}_{i}&\xrightleftharpoons[\beta_{-y}]{\beta_{y}}&\mathrm{Y}_{ic},\end{array}$

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$:

$\mathrm{X}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{X}_{j},\quad\mathrm{Y}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{Y}_{j},\quad
i\neq j,\;i,j=1,..,N.$ (72)

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

$d=(d^{X},d^{Y}):=(d_{1}^{X},..,d_{N}^{X},d_{1}^{Y},...,d_{N}^{Y})\in\left\\{0,..,4^{N}-1\right\\},$

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:

$\displaystyle\begin{array}[]{rcl}\emptyset&\xrightleftharpoons[k_{-x}]{d_{i}^{X}k_{x}}&\mathrm{X}_{i},\\\
n\mathrm{X}_{i}&\xrightleftharpoons[\beta_{-x}]{\beta_{x}}&\mathrm{X}_{ic},\end{array}\begin{array}[]{rcl}\emptyset&\xrightleftharpoons[k_{-y}]{d_{i}^{Y}k_{y}}&\mathrm{Y}_{i},\\\
n\mathrm{Y}_{i}&\xrightleftharpoons[\beta_{-y}]{\beta_{y}}&\mathrm{Y}_{ic},\end{array}$
$\displaystyle\mathrm{X}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{X}_{j},\quad\mathrm{Y}_{i}\xrightleftharpoons[\Omega]{\Omega}\mathrm{Y}_{j},\quad
i,j=1,..,N.$

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:

$\left\\{\left(\frac{k_{x}}{k_{-x}}\omega_{d_{1}}^{X},\dots,\frac{k_{x}}{k_{-x}}\omega_{d_{N}}^{X},\frac{k_{y}}{k_{-y}}\omega_{d_{1}}^{Y},\dots,\frac{k_{y}}{k_{-y}}\omega_{d_{N}}^{Y}\right):d=(d^{X},d^{Y})\in\left\\{0,..,4^{N}-1\right\\}\right\\},$
(73)

where

$\omega_{d_{i}}^{X}=\frac{\sum_{i=1}^{N}\bar{d}_{i}^{X}+\bar{d}_{i}^{X}k_{-x}/\Omega}{N+k_{-x}/\Omega},i=1,..,N,$
(74)

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:

$\mathrm{D}_{0}^{xi}\xrightleftharpoons[\alpha_{-x}]{\alpha_{x}\rho_{d_{i}}^{X}}\mathrm{D}_{1}^{xi},\quad\mathrm{D}_{0}^{yi}\xrightleftharpoons[\alpha_{-y}]{\alpha_{y}\rho_{d_{i}}^{Y}}\mathrm{D}_{1}^{yi},$
(75)

where

$\rho_{d_{i}}^{X}=\mathbb{E}[\mathrm{X}_{ic}|D(t)=d]=\frac{1}{n!}\frac{\beta_{x}}{\beta_{-x}}\left(\frac{k_{x}}{k_{-x}}\omega_{d_{i}}^{X}\right)^{n},$