<title>Introduction</title> <p>The regulation of gene activity is essential for the proper functioning of cells,
which employ a variety of molecular mechanisms to control gene expression. Despite
this, there is considerable variation in the precise number and timing of protein
molecules that are produced for a given gene under any particular set of
circumstances. This is because gene expression is fundamentally a
“noisy” process, subject to a number of sources of randomness.
Some of these are <italic>intrinsic</italic> to the biochemical reactions that
comprise the transcription and translation of a particular gene ##REF##12183631##[1]##,##REF##16179466##[2]##. Several
of the reactions involve very small numbers of molecules. There are only one or two
copies of the DNA for the gene, and in its vicinity inside the cell there are likely
to be only a few copies of the relevant transcription factors and of RNA polymerase.
Similarly, for each mRNA molecule, the processes of ribosome binding and of mRNA
degradation are typically highly stochastic.</p> <p>Recent advances in experimental technology have shown that such single molecule
effects can lead to protein production occurring in bursts of varying size, each due
to a single transcription factor binding event ##REF##9023339##[3]##,##REF##11438714##[4]##. Other sources of
variability are <italic>extrinsic</italic> to the specific reactions, and include
fluctuations in relevant metabolites, polymerases, ribosomes, etc. ##REF##12183631##[1]##,##REF##16179466##[2]##. These
will not be considered further here.</p> <p>It is of considerable interest to determine the various contributions of such
different sources of variability. Within the last few years, experimental techniques
for addressing this question have increasingly become available. Elowitz et al.
##REF##12183631##[1]##
observed fluctuations in the expression level of genes tagged both with cyan and
yellow fluorescent proteins in monoclonal <italic>Escherichia coli</italic> cells
under identical environmental conditions. Similar work was carried out by Raser and
O'Shea ##REF##15166317##[5]## in the eukaryote <italic>Saccharomyces
cerevisiae</italic>. Such dual-reporter experiments are able to distinguish between
intrinsic and extrinsic sources of stochasticity. More recently, single molecule
data has become available ##REF##16543458##[6]##,##REF##16541077##[7]##, which monitors the expression of a gene a single
protein at a time and provides the distribution of the sizes of bursts. It had been
hoped that data of this kind would answer many of the remaining questions about the
origin of noise in gene expression and in particular distinguish between the
different contributions of transcription and translation to intrinsic noise.</p> <p>Intuitively, one might expect that randomness due to transcription would play the
more significant role than translation, since typically there will be more than one
mRNA molecule, and the fluctuations due to translation from each of these might to
some extent average out. To test this hypothesis and to put it on a quantitative
basis, it is necessary to employ mathematical models of gene expression. These also
provide a valuable tool for the analysis of experimental data, and in particular of
the burst size distributions reported in the literature, e.g., ##REF##16543458##[6]##,##REF##16541077##[7]##.</p> <p>A great deal of work has gone into modelling gene expression in both prokaryotic and
eukaryotic systems, with some of the earliest papers predicting fluctuations in mRNA
and protein levels published 30 years ago ##REF##607033##[8]##,##REF##96307##[9]##. McAdams and Arkin ##REF##9023339##[3]##
provided the first model of bursting at the translation level. They showed that the
number of protein molecules produced by a single mRNA transcript is described well
by a model which considers whether the next event is the production of a further
protein, or the degradation of the mRNA molecule. Such competitive binding between
ribosomes and RNase results in a geometric distribution for the protein number. Such
an analysis can also be applied to transcription following the binding of a
transcription factor to a gene and also results in a geometric distribution. The
joint analysis of these two stochastic processes forms the basis of the present
paper.</p> <p>The integration of simple stochastic (Markov) models of transcription factor, RNA
polymerase, ribosome and RNase binding leads to what is now widely regarded as the
standard model of gene expression for prokaryotes ##REF##11438714##[4]##. The analysis of this
model using a master equation allows the determination of the moments of the
distribution of the number of protein molecules when the system is in steady state.
Further analysis of this equilibrium distribution was carried out by Paulsson ##REF##14749823##[10]##–##REF##10852944##[12]## who used the master
equation and the fluctuation–dissipation theorem to obtain predictions
about the mean and variances of molecule numbers and lifetimes and the contribution
made by transcriptional and translational bursting. Other studies have been carried
out by Höfer ##REF##16362909##[13]## who used a rapid-equilibrium approximation to
compare mRNA levels for genes with one and two active alleles, and by Friedman et
al. ##REF##17155441##[14]##. The drawback of these approaches is that the master
equation that describes the temporal evolution of the probability distribution of
protein (and mRNA) numbers is too complex to be solved analytically. Furthermore,
the burst size distribution necessary for comparison with recent experimental data
##REF##16543458##[6]##,##REF##16541077##[7]## cannot be
obtained directly from the master equation. Such difficulties with master equation
based approaches are exacerbated in the case of more complex models of gene
expression such as multi-step models that include intermediate stages such as the
formation of DNA–RNA polymerase complexes, phosphorylation events, and
mRNA–ribosome binding. Both deterministic and stochastic simulation
studies of these models have been performed, e.g., ##REF##12087127##[15]## and ##REF##15721042##[16]##, but none of these
approaches have been useful for the analysis of burst size data.</p> <p>In the present work we avoid the problems associated with the master equation
approach, which are at least in part due to the explicit incorporation of time
evolution. Instead, we ignore time and directly derive an expression for the burst
size distribution by extending the analysis of ##REF##9023339##[3]##. In many ways this
approach is similar to that used for the analysis of multi-stage queues ##UREF##1##[17]##. The
distribution of the number of mRNA molecules produced in a single burst is geometric
and the distribution of the number of protein molecules produced by a single mRNA is
also geometric ##REF##9023339##[3]##. The overall burst size distribution is therefore
given by the compound distribution of two geometric distributions ##UREF##1##[17]##. This
can be readily computed using generating functions ##UREF##1##[17]## and is itself
<italic>not</italic> geometric. However, experimentally it is not possible to
detect bursts that produce no protein molecules at all, and therefore the published
data ##REF##16543458##[6]##,##REF##16541077##[7]## are in fact
the relevant conditional distributions, assuming at least one protein molecule is
produced in a burst. Surprisingly, it turns out that when we condition the compound
distribution in this way, we again obtain a geometric distribution. This is
determined by a single parameter, which we can derive in terms of physically
meaningful constants such as binding and unbinding rates. This shows that different
combinations of noise levels in the translation and transcription parts of the
process can give the same overall burst size distribution. Mathematically, this
means that the standard model of gene expression (described in detail below) is
<italic>nonidentifiable</italic>
##UREF##2##[18]##,##UREF##3##[19]## from
burst size data alone. This in turn implies that it is not possible to identify the
relative contributions of translation and transcription to the burst size
distribution of protein numbers only using this data.</p> <p>We also show that our approach is applicable to a variety of more detailed models
that incorporate additional steps to provide more realistic descriptions of
expression ##REF##15721042##[16]##. These still yield a single parameter geometric
conditional distribution. This shows that within the context of a very large class
of models, experimental burst size data on its own cannot identify the relative
contributions of different reactions to the overall noise level. However, by
simulating the equilibrium distribution of protein numbers for different parameter
combinations giving the same burst size distribution we demonstrate that a
combination of burst size distribution and equilibrium distribution can discern
different sources of noise. The difficulty with such an approach is that the
determination of the equilibrium distribution requires the knowledge of two
additional kinetic parameters: the transcription factor binding rate and the protein
degradation rate. Estimates of these are not easy to obtain independently, so that
we now have to estimate six unknown parameters from the combined burst size and
equilibrium distribution data. Initial simulations (not shown here) suggested that
it is difficult to do this reliably.</p> <p>It is possible however, by using independent estimates of one of the parameters to
reduce the parameter space from six to five dimensions. Using the relationship
between the remaining parameters determined from the burst size distribution allows
the elimination of a further parameter, leaving four kinetic parameters to be
estimated from the equilibrium distribution. We show below that by using the
Nelder-Mead algorithm to maximize the empirical likelihood, useful estimates of the
four remaining parameters can be obtained. We carry out this process twice, first
using independent measurements of the mRNA degradation rate and then of the protein
half-life. In the first case we obtain unrealistically short estimates of the
protein half-life, and in the second a considerably faster mRNA degradation. This
suggests that when in the repressed state, mRNA may be degraded at a faster rate
than when the gene is active.</p> <p>In principle, this method can be applied to any gene where burst size and equilibrium
distributions are available, providing a new approach to the estimation of
parameters estimates for the ever more sophisticated models increasingly being used
in computational biology.</p>
<title>Methods</title> <title>The Standard Model of Gene Expression</title> <p>In the so called “standard model” of gene expression, ##FIG##0##Figure 1##, an inactive gene can
be activated by a promoter or transcription factor. This allows molecules of RNA
polymerase to bind and produce mRNA. This in turn can bind to ribosomes leading
to the production of protein molecules. Eventually the transcription factor
unbinds, terminating the production of mRNA, and each mRNA molecule is degraded,
which stops protein production.</p> <p>Each of these processes is modelled as a transition in a continuous time Markov
chain with a particular rate. Such a rate is interpreted as the probability of
an event occurring in a unit time interval. Thus, if we denote the rate of
transcription factor binding by <italic>α</italic>
<sub>0</sub> then the
probability of this occurring in an interval of length
<italic>δt</italic>, assuming that the transcription factor is not bound
at the start of the interval, is
<italic>α</italic>
<sub>0</sub>
<italic>δt</italic>.
Integrating over time, this means that the probability of the event having
happened by time <italic>t</italic>, is , whilst the average time for the event to happen is
1/<italic>α</italic>
<sub>0</sub>. The same holds for the other
transitions in the model, with the rate of transcription factor unbinding
denoted by <italic>β</italic>
<sub>1</sub>. Whilst the transcription
factor is bound, RNA polymerase binds at a rate
<italic>α</italic>
<sub>1</sub>, and each such binding event is assumed
to produce one molecule of mRNA. More detailed models that allow the polymerase
to unbind before it has produced mRNA are considered later and will have no
effect on our overall conclusions.</p> <p>Each mRNA molecule binds to a ribosome at rate
<italic>α</italic>
<sub>2</sub> and is degraded at rate
<italic>β</italic>
<sub>2</sub>. When the last mRNA has decayed no
more protein will be produced. We define the number of proteins produced between
the transcription factor binding and the last mRNA decaying as a
“burst”. Note that since a burst begins once the
transcription factor has bound, we expect the distribution of burst sizes to be
independent of the transcription factor binding rate
<italic>α</italic>
<sub>0</sub>. This is confirmed by the rigorous
derivation below.</p> <p>Mathematically, the Standard Model of Gene Expression is a continuous time Markov
chain model. Each particular combination of number of mRNA molecules, number of
protein molecules and state of binding of the transcription factor constitutes a
single state of the model. It is possible to derive an (infinite) set of coupled
ordinary differential equations (called the Kolmogorov forward equations or
master equation) that govern the probability at any given time of the system
being in any given state. However, the analysis of a such a complex set of
equations is difficult. On the other hand, using the same approach as for
multi-stage queues, it is relatively easy to derive the distribution of protein
burst sizes.</p> <title>The Component mRNA and Protein Distributions</title> <p>We begin with the analysis of McAdams and Arkin ##REF##9023339##[3]## for the
distribution of the number of proteins produced by a single mRNA molecule. If a
certain number (possibly 0) of protein molecules has been produced, the
probability that the next event in which the mRNA molecule participates is the
production of another protein molecule is
<italic>p</italic> = <italic>α</italic>
<sub>2</sub>/<italic>α</italic>
<sub>2</sub>+<italic>β</italic>
<sub>2</sub>)
(see ##SUPPL##0##Text
S1## for derivation). Conversely, the probability that the next event is
the degradation of the mRNA molecule is
1−<italic>p</italic> = <italic>β</italic>
<sub>2</sub>/(<italic>α</italic>
<sub>2</sub>+<italic>β</italic>
<sub>2</sub>).
In order to produce precisely <italic>n</italic> molecules of protein, we need
<italic>n</italic> events of the first type to occur, followed by a final
degradation event. The probability of this happening is
<italic>p<sup>n</sup></italic>(1−<italic>p</italic>), giving the
distribution <italic>Q</italic>(<italic>n</italic>) of the number of protein
molecules produced by a single mRNA molecule
</p> <p>Here
<italic>A</italic>
<sub>2</sub> = <italic>α</italic>
<sub>2</sub>/<italic>β</italic>
<sub>2</sub>
is the expectation of <italic>Q</italic>. Contrasting this with ##REF##9023339##[3]##,
the parameter <italic>A</italic>
<sub>2</sub> defining the distribution is now
expressed in terms of physically measurable rate constants. Exactly the same
argument applies to the distribution of the number of RNA molecules produced
between the successive binding and unbinding of the transcription factor. In
particular, the probability of producing one more mRNA molecule before the
transcription factor unbinds is
<italic>α</italic>
<sub>1</sub>/(<italic>α</italic>
<sub>1</sub>+<italic>β</italic>
<sub>1</sub>)
and the probability of the transcription factor unbinding is
<italic>β</italic>
<sub>1</sub>/(<italic>α</italic>
<sub>1</sub>+<italic>β</italic>
<sub>1</sub>).
In order to produce precisely <italic>m</italic> mRNA molecules before the
transcription factor unbinds we need <italic>m</italic> independent production
events with probability
<italic>α</italic>
<sub>1</sub>/(<italic>α</italic>
<sub>1</sub>+<italic>β</italic>
<sub>1</sub>),
followed by the unbinding event with probability
<italic>β</italic>
<sub>1</sub>/(<italic>α</italic>
<sub>1</sub>+<italic>β</italic>
<sub>1</sub>).</p> <p>Thus the probability distribution, <italic>R</italic>(<italic>m</italic>), of the
number of mRNA molecules produced in one burst iswhere
<italic>A</italic>
<sub>1</sub> = <italic>α</italic>
<sub>1</sub>/<italic>β</italic>
<sub>1</sub>
is the expectation of <italic>R</italic>(<italic>m</italic>). In order to derive
the overall protein burst size distribution for the Standard Model in ##FIG##0##Figure 1## we need the
probability generating functions ##UREF##1##[17]## of the
distributions <italic>Q</italic>(<italic>n</italic>) and
<italic>R</italic>(<italic>m</italic>) which we denote as
<italic>Q</italic>*(<italic>z</italic>) and
<italic>R</italic>*(<italic>z</italic>), respectively. These are
simply obtained by summing the relevant geometric seriesand
</p>
<title>Results</title> <title>The Compound Protein Burst Size Distribution</title> <p>The distribution <italic>P</italic>(<italic>n</italic>) of the total number of
proteins produced in a single burst is simply the compound distribution of
<italic>R</italic> and <italic>Q</italic>
##UREF##1##[17]##.
This is easily computed using probability generating functions (see below), and
is not a geometric distribution. However, it is of relatively little interest
since it includes the possibility that the transcription factor unbinds before
any proteins have been produced (either because no mRNA is produced, or because
this mRNA is degraded before binding to a ribosome). Such events cannot be
observed in the experimental protocol used in ##REF##16543458##[6]##,##REF##16541077##[7]##, and hence
<italic>P</italic>(<italic>n</italic>) cannot be directly compared to the
data in these papers. However, we can re-scale
<italic>P</italic>(<italic>n</italic>) to give the probability distribution
<italic>Pˆ</italic>(<italic>n</italic>) = <italic>P</italic>(<italic>n</italic>)/(1−<italic>P</italic>(0))
of protein numbers conditional on at least one protein being produced. An
approximate calculation of this distribution was given in the supplementary
material of ##REF##16541077##[7]##. This replaced the discrete geometric distribution
<italic>Q</italic>(<italic>n</italic>) by a continuous exponential
distribution of the same mean and then used the Laplace transform to obtain the
(continuous approximation to the) compound distribution. Here we present an
exact derivation for the discrete distribution using generating functions (which
are closely related to the Laplace transform). Furthermore we relate the
parameter of the final burst size distribution to the original kinetic
parameters <italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub>.</p> <p>Thus, let <italic>X</italic>
<sup>(<italic>i</italic>)</sup> be the random
variable, with distribution <italic>Q</italic>(<italic>n</italic>), giving the
number of proteins produced by the <italic>i</italic>th mRNA transcript and let
<italic>Y</italic> be a random variable, with distribution
<italic>R</italic>(<italic>n</italic>) giving the number of mRNA molecules
produced. Then the random variablegives the total number of proteins in a burst. Denote the
distribution of <italic>X</italic> by <italic>P</italic>(<italic>n</italic>),
with generating function <italic>P</italic>*(<italic>z</italic>). Then a
standard result on generating functions of compound distributions ##UREF##1##[17]## givesTo obtain the distribution conditional on at least one protein
molecule being produced, we subtract <italic>P</italic>*(0) and
normalise (divide) by 1−<italic>P</italic>*(0) to giveThis is the generating function of a conditional geometric
distribution with (dimensionless) parameter
<italic>Â</italic>
<sub>2</sub> = <italic>A</italic>
<sub>2</sub>(1+<italic>A</italic>
<sub>1</sub>),
so that <italic>Pˆ</italic>(<italic>n</italic>) has the distributionwhere the parameter <italic>Â</italic>
<sub>2</sub> can
be expressed in terms of the mean number <italic>A</italic>
<sub>1</sub> of mRNA
molecules produced and the mean number <italic>A</italic>
<sub>2</sub> of protein
molecules produced from a single mRNA molecule as
We thus see that the burst size distribution is determined by a
single parameter, and that many different combinations of the parameters
<italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>, <italic>β</italic>
<sub>1</sub>,
and <italic>β</italic>
<sub>2</sub> will lead to the same burst size
distribution. In mathematical language this says that the Standard Model with
parameters <italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub> is <italic>nonidentifiable</italic> from
burst size data. In fact we can only estimate a single parameter (or a single
linear combination) and the three remaining parameters can be arbitrarily
chosen.</p> <title>Burst Distributions for Extensions of the Standard Model</title> <p>It might be hoped that such nonidentifiability is a particular pathology of the
Standard Model. We thus next consider a number of generalisations of this model,
which provide a more detailed description of the process of gene expression. We
find that for a wide range of generalisations we can still derive the burst size
distribution in a similar manner the above. It turns out to be geometric in each
case and hence all such models are also nonidentifiable.</p> <p>One common extension is to include an additional step in the model of the
transcription process ##REF##16362909##[13]##, as shown in ##FIG##1##Figure 2##. This accounts for the fact that
after the transcription factor has bound, one still requires the RNA polymerase
to bind to the transcription initiation complex, and this may not always happen
successfully. A similar modification could be made to the translation loop to
describe the binding of the mRNA transcript to the ribosome in more detail. Both
of these additions can be considered individually, or in combination.</p> <p>Doing this results in distributions <italic>R</italic> and <italic>Q</italic>
which are still geometric, but with the parameters
<italic>A</italic>
<sub>1</sub> and <italic>A</italic>
<sub>2</sub> given by more
complex combinations of the individual rates. We illustrate this for the
transcription loop, where we find that in order to produce exactly
<italic>m</italic> mRNA molecules, the system can pass through state
<italic>G</italic>* any number
<italic>i</italic>≥<italic>m</italic> times. On
<italic>i</italic>−<italic>m</italic> of these occasions the
polymerase unbinds before an mRNA molecule is produced, returning to
<italic>G</italic> with rate <italic>δ</italic>
<sub>1</sub>, and on the
remaining <italic>m</italic> occasions an mRNA molecule is produced, with rate
<italic>α</italic>
<sub>1</sub>. The <italic>m</italic> productive
steps can be interspersed in any order amongst the <italic>i</italic> visits,
giving  possible choices. The probability of producing
<italic>m</italic> mRNA molecules is thuswith <italic>A</italic>
<sub>1</sub> now given by
<italic>A</italic>
<sub>1</sub> = <italic>α</italic>
<sub>1</sub>
<italic>γ</italic>
<sub>1</sub>/<italic>β</italic>
<sub>1</sub>(<italic>α</italic>
<sub>1</sub>+<italic>δ</italic>
<sub>1</sub>).
A similar derivation holds for the translation loop. We see that carrying out
either or both of these modifications still results in a geometric distribution
in the form of Equation 4 for <italic>Pˆ</italic>(<italic>n</italic>),
with
<italic>Â</italic>
<sub>2</sub> = <italic>A</italic>
<sub>2</sub>(1+<italic>A</italic>
<sub>1</sub>),
but <italic>A</italic>
<sub>1</sub> and <italic>A</italic>
<sub>2</sub> now given
by
<italic>A</italic>
<sub>1</sub> = <italic>α</italic>
<sub>1</sub>
<italic>γ</italic>
<sub>1</sub>/<italic>β</italic>
<sub>1</sub>(<italic>α</italic>
<sub>1</sub>+<italic>δ</italic>
<sub>1</sub>)
and
<italic>A</italic>
<sub>2</sub> = <italic>α</italic>
<sub>2</sub>
<italic>γ</italic>
<sub>2</sub>/<italic>β</italic>
<sub>2</sub>(<italic>α</italic>
<sub>2</sub>+<italic>δ</italic>
<sub>2</sub>).
As a consequence the overall conditional protein size distribution,
<italic>Pˆ</italic>(<italic>n</italic>), will still be given by
Equation 4, with the parameter
<italic>Â</italic>
<sub>2</sub> = <italic>A</italic>
<sub>2</sub>
<italic>A</italic>
<sub>1</sub>+<italic>A</italic>
<sub>2</sub>
as before.</p> <p>An alternative generalisation is to add additional loops with the same structure
as the current transcription and translation loops. We prove in the Supporting
Information (##SUPPL##0##Text S1##) that if we have <italic>k</italic>−1 such loops,
the final conditional protein size distribution <italic>Pˆ</italic>
<italic><sub>k</sub></italic>(<italic>n</italic>) will still be geometric.</p> <p>We thus conclude that all of these models yield the same geometric protein burst
size conditional distribution, determined by a single parameter. In particular,
models which include additional steps to account for DNA–RNAP complex
formation and mRNA-ribosome complex formation give distributions that are
mathematically indistinguishable from those from the Standard Model. It is thus
impossible to differentiate between these models using experimentally observed
burst size distributions. Similarly we cannot use such data to differentiate
between the contributions to noisy gene expression from transcriptional versus
translational bursting.</p> <title>Comparison with Burst Size Data</title> <p>We can compare the probability distribution derived above directly with
experimental data. We consider recently published data of burst sizes for two
fluorescently tagged proteins in the bacterium <italic>Escherichia coli</italic>
##REF##16543458##[6]##,##REF##16541077##[7]##. In
##REF##16543458##[6]##, a
novel fluorescent imaging technique is used to determine the distribution of
protein molecules per transcription factor binding event in live <italic>E.
coli</italic> cells. The specific protein studied was a fusion of a yellow
fluorescent protein variant (Venus) with the membrane protein Tsr. The
<italic>tsr-venus</italic> gene is incorporated into the <italic>E.
coli</italic> chromosome, replacing the <italic>lacZ</italic> gene. This
modified gene is then under the control of the <italic>lac</italic> promoter. In
a second publication ##REF##16541077##[7]##, the same group used a different imaging
technique to determine the distribution of protein molecules per transcription
factor binding event of <italic>β</italic>-gal in live <italic>E.
coli</italic> cells.</p> <p>Such experimental data can be compared to the predicted distribution
<italic>Pˆ</italic>(<italic>n</italic>) in two ways. One
possibility is to use maximum likelihood estimation to find the value of
<italic>Â</italic>
<sub>2</sub> for which
<italic>Pˆ</italic>(<italic>n</italic>) best fits the data. This
is illustrated in ##FIG##2##Figure 3##,
which shows that it is possible to obtain excellent agreement between the
theoretical and experimental distributions. The estimated value of
<italic>Â</italic>
<sub>2</sub> for Tsr-Venus is
<italic>Â</italic>
<sub>2</sub> = 3.57,
whilst for <italic>β</italic>-gal,
<italic>Â</italic>
<sub>2</sub> = 20.96.
The difference in magnitude between these two estimates may be partially due to
the fact that <italic>β</italic>-gal is only active as a tetramer. Thus,
each burst of activation measured experimentally (and thus available for
fitting) corresponds to the production of 4 monomers. The disadvantage of
fitting the model in this way is it can only provide an estimate of the single
parameter <italic>Â</italic>
<sub>2</sub>, but not of the underlying
kinetic parameters <italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub>.</p> <p>An alternative approach to verifying the model would be to obtain independent
estimates of the model parameters from which we can calculate
<italic>Â</italic>
<sub>2</sub> using Equation 6. The resulting
geometric distribution can then be compared to the observed burst size data.
Unfortunately, as is common for most models in cell and molecular biology,
direct experimental measurements of many of these rates are not available. For
the <italic>β</italic>-gal data, <italic>β</italic>
<sub>2</sub>
can be obtained from the reported mRNA half life ##REF##16541077##[7]##,##UREF##4##[20]##, but the other
three parameters corresponding to the off-rate of the transcription factor and
to the binding rates of RNA polymerase to DNA and of mRNA to ribosome
respectively are not available.</p> <title>Application to Experimental Data</title> <title>Incorporating steady state distribution data</title> <p>We thus conclude that we can neither estimate all the kinetic parameters
<italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub> from the burst size data, nor measure
them by other means. However, experience suggests that by supplementing the
burst size distribution with other experimental data it may be possible to
overcome the nonidentifiability of these parameters. This is reinforced by
the observation that parameter combinations that lead to the same
<italic>Â</italic>
<sub>2</sub> and hence the same burst size
distribution can yield quite different steady-state distributions, as shown
for example in ##FIG##3##Figure 4##.
The two steady state distributions shown have different choices of
<italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub> that yield the same value for
<italic>Â</italic>
<sub>2</sub>, and hence the same burst size
distribution. However, the two steady state distributions are clearly
different. This shows that steady state distribution data should allow us to
distinguish between different combinations of parameters with the same
<italic>Â</italic>
<sub>2</sub>, and hence potentially identify
some or all of these parameters.</p> <title>Empirical likelihood estimation</title> <p>The main difficulty with such an approach is the lack of analytic expressions
for the steady state distribution, making it impossible to derive an
explicit formula for the likelihood. Instead one has to compute an estimate
of the equilibrium distribution using simulations of the reaction network
##UREF##4##[20]## and then use these to derive an empirical
likelihood by comparing to the experimental data. This can then be maximized
in the usual way.</p> <p>We applied this approach to the data from ##REF##16541077##[7]##, which presents
both burst size and steady state distributions for the same experimental
system. In order to fully specify the steady state distribution, we need two
additional parameters: the rate of transcription factor binding
<italic>α</italic>
<sub>0</sub> and the rate of protein decay
<italic>β</italic>
<sub>3</sub>. These do not enter into the
expressions for the burst size distribution, and were assumed to be known
(and fixed) for the simulations shown in ##FIG##3##Figure 4##. In the absence of independent
estimates of these parameters for the <italic>β</italic>-gal system,
we explored the possibility of estimating these from the data in ##REF##16541077##[7]##
directly by computing an empirical likelihood using simulation of the model
(see below). We attempted both to maximize this empirical likelihood
directly, and to obtain its distribution using Markov Chain Monte Carlo
sampling. Neither of these approaches were successful with the full six
parameter model (results not shown).</p> <p>We can however, make use of independent estimates of parameters in the model
to reduce the dimensionality of the parameter space. In effect this
constrains the orginal optimization to a lower dimensional sub-space. We
applied this approach with two different choices of parameter: the mRNA
degradation rate <italic>β</italic>
<sub>2</sub> and the protein
degradation rate <italic>β</italic>
<sub>3</sub>.</p> <title>Constraining on the mRNA degradation rate</title> <p>We chose first to make use of the wide availability of estimates of the value
of <italic>β</italic>
<sub>2</sub>, the rate of mRNA degradation.
Since we also have the burst size data, we first estimate
<italic>Â</italic>
<sub>2</sub> and then use Equation 6 to
obtain an expression for <italic>α</italic>
<sub>2</sub> in terms of
<italic>α</italic>
<sub>1</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub>. We are left with the four
dimensional parameter space <italic>α</italic>
<sub>0</sub>,
<italic>α</italic>
<sub>1</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>3</sub>. At each point in this space, we
simulate the model using the Gillespie algorithm to given an empirical
estimate of the probability
<italic>P<sub>n</sub></italic>(<italic>α</italic>
<sub>0</sub>,<italic>α</italic>
<sub>1</sub>,<italic>β</italic>
<sub>1</sub>,<italic>β</italic>
<sub>3</sub>)
of observing <italic>n</italic> proteins at equilibrium. This gives the
empirical log-likelihoodwhere  is the number of times that <italic>n</italic> proteins
are observed in the experimentally data .</p> <p>This empirical likelihood function can be maximized using any suitable
optimization algorithm. Because it is computed using stochastic simulations
any particular realization of the function is not smooth, making algoritms
that use gradient (or Hessian) information unsuitable. We therefore opted to
use the Nelder-Mead (simplex) method ##REF##16677373##[21]## with the
results shown in ##FIG##4##Figure
5##.</p> <p>We can see that the likelihood has a local maximum at
<italic>L</italic>≈−1800, and the simplex method frequently
gets stuck in this region. However the majority of runs (73.4%)
converge to the presumed global maximum. The means and standard deviations
of the estimated parameter values are shown in ##TAB##0##Table 1##. The value obtained for
<italic>β</italic>
<sub>3</sub> is 0.0297
s<sup>−1</sup>, which corresponds to a half-life of 23 seconds.
This appears to be unrealistically short, since
<italic>β</italic>-galactosidase is a stable protein with a reported
lifetime of hours.</p> <p>One possible explanation of this discrepancy is that the reported mRNA
degradation rate
<italic>β</italic>
<sub>2</sub> = 7.2×10<sup>−3</sup>
is always measured in experimental conditions where the gene is active. On
the other hand, the burst size and equilibrium distributions in ##REF##16541077##[7]## are
obtained under conditions where the gene is suppressed. It is possible that
the mRNA degradation rates are significantly different in the two cases. To
explore this hypothesis, we approached the problem from an alternative
direction, fixing the protein degradation rate
<italic>β</italic>
<sub>3</sub> to correspond to a half-life of 1
hour, and estimating the remaining five parameters, including mRNA
degradation rate <italic>β</italic>
<sub>2</sub>.</p> <title>Constraining on the protein degradation rate</title> <p>We therefore fixed <italic>β</italic>
<sub>3</sub> to
1.92×10<sup>−4</sup> s<sup>−1</sup>,
corresponding to a protein half-life of one hour, and then used same method
as described above to estimate the other parameters
<italic>α</italic>
<sub>0</sub>,
<italic>α</italic>
<sub>1</sub>,
<italic>α</italic>
<sub>2</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub>. We ran 10,000 simulations, as a
relatively low number of runs converged (23.37%), with the others
becoming trapped in a region with physically unrealistic (negative) reaction
rates, and a log-likelihood of <italic>L</italic>≈−2100. Of
the runs which converged, 2057 (88%) converged to a local maximum
at <italic>L</italic>≈−9150, while 279 (12%)
converged to the presumed global maximum at
<italic>L</italic>≈−1100. The results for the runs which
converged can be seen in ##FIG##5##Figure
6##, whilst summary statistics for the runs which converged to the
presumed global maximum are presented in ##TAB##1##Table 2##.</p> <title>Comparison of the estimates</title> <p>The transcription factor binding rate <italic>α</italic>
<sub>0</sub>
is almost unchanged under both assumptions. When we fix the protein
degradation rate in the second set of estimates to a value much lower than
estimated in the first set we find that the transcription rate (i.e., rate
of RNA polymerase binding) <italic>α</italic>
<sub>1</sub>, is
approximately one third of the previous value, decreasing from 0.0017
s<sup>−1</sup> to 0.0006 s<sup>−1</sup>, whilst the
translation rate, <italic>α</italic>
<sub>2</sub> shows an
approximate two-fold increase, from 0.1538 s<sup>−1</sup> to
0.3352 s<sup>−1</sup>. It is intuitively reasonable that such a
combination of decreasing transcription and increasing translation leads to
the same overall level of protein expression. The parameter
<italic>β</italic>
<sub>1</sub>, the rate at which the
transcription factor unbinds increases slightly, leading to shorter bursts.
The increase in the mRNA degradation rate,
<italic>β</italic>
<sub>2</sub> from the original assumption of 0.007
s<sup>−1</sup> to the estimate of 0.0161
s<sup>−1</sup> (corresponding to an mRNA half life of
approximately 43 seconds) suggests that when expression of the gene is being
strongly repressed as in this situation, there may well be active
degradation of the mRNA. It would be interesting to experimentally
investigate this biologically significant prediction.</p>
<title>Discussion</title> <p>We have shown that it is possible to use results from queuing theory to derive the
burst size distribution of protein molecules produced by a single transcription
factor binding event in terms of physically measurable kinetic rate constants for
both the simplest model of gene expression, the so-called Standard Model, and for a
number of natural extensions.</p> <p>Furthermore, we have shown that the mathematical form of these models is
nonidentifiable, and all such burst size distributions are actually determined by a
single parameter. This implies that it is impossible to use burst size data alone to
determine the relative contributions of transcription and translation to the
variability in gene expression.</p> <p>One possible way of overcoming this limitation is to use a combination of burst size
data and steady-state data. However, this requires estimates of a further two
parameters (which are not needed when using burst-size data alone). We were unable
to estimate all six parameters directly from the combined data. However, using
independent estimates of either the mRNA lifetime or the protein lifetime reduces
the number of parameters by one, and enables successfully estimation of the
remaining five parameters by maximizing an empirical likelihood using the
Nelder-Mead simplex algorithm. Although this suffers from the common problem of
occasional convergence to a local maximum, by using computing repeated estimates it
was possible to identify and exclude such cases and hence obtain good estimates of
the desired five kinetic parameters under the different constraints.</p>
<p>Conceived and designed the experiments: PJI MPHS JS. Wrote the paper: PJI JS.
Performed the analysis: PJI.</p> <p>Over the last few years, experimental data on the fluctuations in gene activity
between individual cells and within the same cell over time have confirmed that
gene expression is a “noisy” process. This variation is in
part due to the small number of molecules taking part in some of the key
reactions that are involved in gene expression. One of the consequences of this
is that protein production often occurs in bursts, each due to a single promoter
or transcription factor binding event. Recently, the distribution of the number
of proteins produced in such bursts has been experimentally measured, offering a
unique opportunity to study the relative importance of different sources of
noise in gene expression. Here, we provide a derivation of the theoretical
probability distribution of these bursts for a wide variety of different models
of gene expression. We show that there is a good fit between our theoretical
distribution and that obtained from two different published experimental
datasets. We then prove that, irrespective of the details of the model, the
burst size distribution is always geometric and hence determined by a single
parameter. Many different combinations of the biochemical rates for the
constituent reactions of both transcription and translation will therefore lead
to the same experimentally observed burst size distribution. It is thus
impossible to identify different sources of fluctuations purely from protein
burst size data or to use such data to estimate all of the model parameters. We
explore methods of inferring these values when additional types of experimental
data are available.</p> <title>Author Summary</title> <p>Recent experimental data showing fluctuations in gene activity between individual
cells and within the same cell over time confirm that gene expression is a
“noisy” process. This variation is partly due to the small
number of molecules involved in gene expression. One consequence is that protein
production often occurs in bursts, each due to the binding of a single
transcription factor. Recently, the distribution of the number of proteins
produced in such bursts has been experimentally measured, offering a unique
opportunity to study the relative importance of different sources of noise in
gene expression. We derive the theoretical probability distribution of these
bursts for a wide variety of gene expression models. We show a good fit between
our theoretical distribution and experimental data and prove that, irrespective
of the model details, the burst size distribution always has the same shape,
determined by a single parameter. As different combinations of the reaction
rates lead to the same observed distribution, it is impossible to estimate all
kinetic parameters from protein burst size data. When additional data, such as
protein equilibrium distributions, are available, these can be used to infer
additional parameters. We present one approach to this, demonstrating its
application to published data.</p>
<title>Supporting Information</title>
<fig id="pcbi-1000192-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g001</object-id><label>Figure 1</label><caption><title>The standard gene expression model.</title><p>An inactive sequence of DNA and a transcription factor bind to produce an
active gene <italic>G</italic>. This produces mRNA, denoted by
<italic>M</italic> at a rate <italic>α</italic>
<sub>1</sub>, and
in turn the mRNA produces protein at rate
<italic>α</italic>
<sub>2</sub>. Eventually, the
transcription factor will unbind (at rate
<italic>β</italic>
<sub>1</sub>), and the gene will become
inactive again. Each copy of mRNA produced will also be degraded (at
rate <italic>β</italic>
<sub>2</sub>).</p></caption></fig> <fig id="pcbi-1000192-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g002</object-id><label>Figure 2</label><caption><title>Diagram of the generalised situation in which intermediate,
reversible stages are introduced.</title><p>Here, G represents an active gene, G* an active gene with a bound
RNA polymerase, M an mRNA molecule, M* an mRNA molecule bound to
a ribosome, P a protein, and <italic>S</italic>
<sub>0</sub> states which
correspond to transcription factor unbinding and mRNA transcript
decay.</p></caption></fig> <fig id="pcbi-1000192-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g003</object-id><label>Figure 3</label><caption><title>Comparison of the distribution of experimentally measured burst sizes
for the proteins <italic>Tsr-Venus</italic> (A) ##REF##16543458##[6]## and for the
<italic>β-gal</italic> (B) ##REF##16541077##[7]## with the
standard model of gene expression.</title><p>In both cases the blue line shows the best fit of the model to the data,
obtained using the method of maximum likelihood giving
<italic>Â</italic>
<sub>2</sub> = 3.57
for <italic>Tsr-Venus</italic> and
<italic>Â</italic>
<sub>2</sub> = 20.96
for <italic>β-gal</italic>. The error bars show the upper and
lower bounds of the 95% confidence interval for the fitted
parameter.</p></caption></fig> <fig id="pcbi-1000192-g004" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g004</object-id><label>Figure 4</label><caption><title>Simulations of steady-state protein expression levels.</title><p>For A we have
<italic>α</italic>
<sub>1</sub> = 0.018
and
<italic>β</italic>
<sub>1</sub> = 0.086
and for B we have
<italic>α</italic>
<sub>1</sub> = 0.009
and
<italic>β</italic>
<sub>1</sub> = 0.043,
resulting in the same <italic>Â</italic>
<sub>2</sub> and
hence identical burst size distributions. Other parameters were
<italic>α</italic>
<sub>0</sub> = 0.012,
<italic>α</italic>
<sub>2</sub> = 0.013,
<italic>β</italic>
<sub>2</sub> = 0.0039,
and
<italic>β</italic>
<sub>3</sub> = 0.0007,
based on previous simulation studies ##REF##11062240##[22]##. The
distributions shown are for a run of 10,000 seconds using the Stocks
implementation of Gillespie's method ##REF##4604615##[23]##, after an initial transient of 10,000
seconds. Previous studies have indicated that the steady state is in
fact attained in under 1000 seconds.</p></caption></fig> <fig id="pcbi-1000192-g005" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g005</object-id><label>Figure 5</label><caption><title>Parameter estimation results with fixed mRNA degradation rate.</title><p>The results of 1000 runs of the Nelder-Mead maximisation of the
log-likelihood for the parameters
<italic>α</italic>
<sub>0</sub>,
<italic>α</italic>
<sub>1</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>3</sub>, with
<italic>α</italic>
<sub>2</sub> determined by the
relationship in Equation 6, and with the mRNA degradation rate
<italic>β</italic>
<sub>2</sub> set to
7.2×10<sup>−3</sup>, corresponding to a
half-life for <italic>β</italic>-gal mRNA of 1.6 mins ##UREF##5##[24]##. The panels in column A show the
estimates of the values of the parameters and the percentage of
times the Nelder-Mead algorithm converged to those values. The
panels in column B are scattergrams of the values of the parameter
estimates against the value of the log-likelihood. Each simulation
is run 10,000 times to simulate a population of 10,000 cells, and
each simulation is run for 5000 reaction steps. The starting values
for the optimisation routine are:
<italic>α</italic>
<sub>0</sub> = 0.01
s<sup>−1</sup>,
<italic>α</italic>
<sub>1</sub> = 0.02
s<sup>−1</sup>,
<italic>β</italic>
<sub>1</sub> = 0.1
s<sup>−1</sup>, and
<italic>β</italic>
<sub>3</sub> = 0.0007
s<sup>−1</sup>, and are based on previous simulation
studies ##REF##15721042##[16]##.</p></caption></fig> <fig id="pcbi-1000192-g006" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.g006</object-id><label>Figure 6</label><caption><title>Parameter estimation results with fixed protein degradation rate.</title><p>The results of 10,000 runs of the Nelder-Mead maximisation of the
log-likelihood for the parameters
<italic>α</italic>
<sub>0</sub>,
<italic>α</italic>
<sub>1</sub>,
<italic>β</italic>
<sub>1</sub>, and
<italic>β</italic>
<sub>2</sub>, with
<italic>α</italic>
<sub>2</sub> determined by the
relationship in Equation 6, and with the protein degradation rate
set, <italic>β</italic>
<sub>3</sub> set to
2.77×10<sup>−4</sup>, consistent with a
half-life for <italic>β</italic>-gal 60 mins. The panels in
column A show the estimates of the values of the parameters and the
percentage of times the Nelder-Mead algorithm converged to those
values. The panels in column B are scattergrams of the values of the
parameter estimates against the value of the log-likelihood. Each
simulation is run 10,000 times to simulate a population of 10,000
cells, and each simulation is run for 5000 reaction steps. The
starting values for the optimisation routine are:
<italic>α</italic>
<sub>0</sub> = 0.01
s<sup>−1</sup>,
<italic>α</italic>
<sub>1</sub> = 0.02
s<sup>−1</sup>,
<italic>β</italic>
<sub>1</sub> = 0.1
s<sup>−1</sup>, and
<italic>β</italic>
<sub>2</sub> = 0.007
s<sup>−1</sup>, and are based on previous simulation
studies ##REF##15721042##[16]##.</p></caption></fig>
<table-wrap id="pcbi-1000192-t001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.t001</object-id><label>Table 1</label><caption><title>Means and standard deviations of estimates of the other
parameters when the mRNA degradation rate
<italic>β</italic>
<sub>2</sub> set to
7.2×10<sup>−3</sup>.</title></caption><alternatives><table frame="hsides" rules="groups"><colgroup span="1"><col align="left" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/></colgroup><thead><tr><td align="left" rowspan="1" colspan="1"/><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>0</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>1</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>2</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>β</italic>
<sub>1</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>β</italic>
<sub>3</sub>
</td></tr></thead><tbody><tr><td align="left" rowspan="1" colspan="1">
<italic>μ</italic>
</td><td align="left" rowspan="1" colspan="1">0.0049</td><td align="left" rowspan="1" colspan="1">0.0017</td><td align="left" rowspan="1" colspan="1">0.1538</td><td align="left" rowspan="1" colspan="1">0.1210
(<italic>τ</italic>
<sub>1/2</sub> = 5.7
s)</td><td align="left" rowspan="1" colspan="1">0.0297
(<italic>τ</italic>
<sub>1/2</sub> = 23.3
s)</td></tr><tr><td align="left" rowspan="1" colspan="1">
<italic>σ</italic>
</td><td align="left" rowspan="1" colspan="1">0.0052</td><td align="left" rowspan="1" colspan="1">0.0003</td><td align="left" rowspan="1" colspan="1">0.0088</td><td align="left" rowspan="1" colspan="1">0.0098</td><td align="left" rowspan="1" colspan="1">0.0098</td></tr></tbody></table></alternatives></table-wrap> <table-wrap id="pcbi-1000192-t002" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000192.t002</object-id><label>Table 2</label><caption><title>Means and standard deviations of the parameter estimates, when
the protein degradation rate <italic>β</italic>
<sub>3</sub>
set to 1.92×10<sup>−4</sup>
(<italic>τ</italic>
<sub>1/2</sub> = 3600
s).</title></caption><alternatives><table frame="hsides" rules="groups"><colgroup span="1"><col align="left" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/></colgroup><thead><tr><td align="left" rowspan="1" colspan="1"/><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>0</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>1</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>α</italic>
<sub>2</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>β</italic>
<sub>1</sub>
</td><td align="left" rowspan="1" colspan="1">
<italic>β</italic>
<sub>2</sub>
</td></tr></thead><tbody><tr><td align="left" rowspan="1" colspan="1">
<italic>μ</italic>
</td><td align="left" rowspan="1" colspan="1">0.0048</td><td align="left" rowspan="1" colspan="1">0.0006</td><td align="left" rowspan="1" colspan="1">0.3352</td><td align="left" rowspan="1" colspan="1">0.1314
(<italic>τ</italic>
<sub>1/2</sub> = 5.3
s)</td><td align="left" rowspan="1" colspan="1">0.0161
(<italic>τ</italic>
<sub>1/2</sub> = 43.1
s)</td></tr><tr><td align="left" rowspan="1" colspan="1">
<italic>σ</italic>
</td><td align="left" rowspan="1" colspan="1">0.0033</td><td align="left" rowspan="1" colspan="1">0.0004</td><td align="left" rowspan="1" colspan="1">0.1476</td><td align="left" rowspan="1" colspan="1">0.0078</td><td align="left" rowspan="1" colspan="1">0.0071</td></tr></tbody></table></alternatives></table-wrap>
<inline-formula></inline-formula> <disp-formula><label>(1)</label></disp-formula> <disp-formula><label>(2)</label></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula><label>(3)</label></disp-formula> <disp-formula></disp-formula> <disp-formula><label>(4)</label></disp-formula> <disp-formula><label>(5)</label></disp-formula> <disp-formula><label>(6)</label></disp-formula> <inline-formula></inline-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <inline-formula></inline-formula> <inline-formula></inline-formula>
<supplementary-material content-type="local-data" id="pcbi.1000192.s001"><label>Text S1</label><caption><p>Derivation of probabilities.</p><p>(0.10 MB PDF)</p></caption></supplementary-material>
<table-wrap-foot><fn id="nt101"><p>These statistics are based on those runs that approached the
global maximum (73.4% of all runs). These were
selected by imposing a threshold at
<italic>L</italic> = −1350
and only considering those runs converging to a larger (more
positive) likelihood. All reaction rates have units of
s<sup>−1</sup>.</p></fn></table-wrap-foot> <table-wrap-foot><fn id="nt102"><p>These statistics are based on those runs which approached the
global maximum (12% of converging runs). These were
selected by imposing a threshold at
<italic>L</italic> = −2500
and only considering those runs converging to a larger (more
positive) likelihood. All reaction rates have units of
s<sup>−1</sup>.</p></fn></table-wrap-foot> <fn-group><fn fn-type="COI-statement"><p>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p>PJI and MPHS acknowledge financial support from the Wellcome Trust. MPHS is
grateful to the Carlsberg Foundation and the Royal Society, UK, for their
generous support. MPHS receives further support through an EMBO Young
Investigator Award. JS is supported by the UK Biotechnology and Biological
Sciences Research Council via the Centre for Integrative Systems Biology at
Imperial College, BB/C519670/1.</p></fn></fn-group>
<graphic xlink:href="pcbi.1000192.g001"/> <inline-graphic xlink:href="pcbi.1000192.e001.jpg" mimetype="image"/> <graphic xlink:href="pcbi.1000192.e002.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e003.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e004.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e005.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e006.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e007.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e008.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e009.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e010.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.e011.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.g002"/> <inline-graphic xlink:href="pcbi.1000192.e012.jpg" mimetype="image"/> <graphic xlink:href="pcbi.1000192.e013.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000192.g003"/> <graphic xlink:href="pcbi.1000192.g004"/> <graphic xlink:href="pcbi.1000192.e014.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pcbi.1000192.e015.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000192.e016.jpg" mimetype="image"/> <graphic xlink:href="pcbi.1000192.g005"/> <graphic id="pcbi-1000192-t001-1" xlink:href="pcbi.1000192.t001"/> <graphic xlink:href="pcbi.1000192.g006"/> <graphic id="pcbi-1000192-t002-2" xlink:href="pcbi.1000192.t002"/>
<media xlink:href="pcbi.1000192.s001.pdf"><caption><p>Click here for additional data file.</p></caption></media>
[{"label": ["11"], "element-citation": ["\n"], "surname": ["Paulsson"], "given-names": ["J"], "year": ["2005"], "article-title": ["Models of stochastic gene expression."], "source": ["Phys Life Rev"], "volume": ["2"], "fpage": ["157"], "lpage": ["175"]}, {"label": ["17"], "element-citation": ["\n"], "surname": ["Bailey"], "given-names": ["NTJ"], "year": ["1964"], "source": ["The Elements of Stochastic Processes"], "publisher-loc": ["New York"], "publisher-name": ["Wiley"]}, {"label": ["18"], "element-citation": ["\n"], "surname": ["Timmer", "Muller", "Swameye", "Sandra", "Klingmuller"], "given-names": ["J", "TG", "I", "O", "U"], "year": ["2004"], "article-title": ["Modelling the nonlinear dynamics of cellular signal transduction."], "source": ["Int J Bifurcat Chaos"], "volume": ["14"], "fpage": ["2069"], "lpage": ["2079"]}, {"label": ["19"], "element-citation": ["\n"], "surname": ["Sontag"], "given-names": ["ED"], "year": ["2002"], "article-title": ["For differential equations with r parameters, 2r+1\nexperiments are enough for identification."], "source": ["J Nonlinear Sci"], "volume": ["12"], "fpage": ["553"], "lpage": ["583"]}, {"label": ["20"], "element-citation": ["\n"], "surname": ["Gillespie"], "given-names": ["DT"], "year": ["1977"], "article-title": ["Exact stochastic simulation of coupled chemical reactions."], "source": ["J Phys Chem"], "volume": ["81"], "fpage": ["2340"], "lpage": ["2361"]}, {"label": ["24"], "element-citation": ["\n"], "surname": ["Nelder", "Mead"], "given-names": ["JA", "R"], "year": ["1965"], "article-title": ["A simplex method for function minimization."], "source": ["Comput J"], "volume": ["7"], "fpage": ["308"], "lpage": ["313"]}]
PLoS Comput Biol. 2008 Oct 10; 4(10):e1000192