The origin of the Poisson distribution in stochastic dynamics of gene expression

Gene regulatory network(GRN) controls the life process by producing and degrading various kinds of proteins that perform important biological functions. It is a well-known fact that the time evolution of mRNA and/or protein molecules in such a network is stochastic [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. Even when various external conditions such as cellular environments are identical, there is always intrinsic noise due to remaining uncontrolled factors that influence the GRN of interest, making its dynamics stochastic. It has been suggested that biological organisms may have evolved to take advantage of such fluctuations [13].

In simple theoretical models of stochastic GRN dynamics, the Poisson distribution,

often appears as a probability distribution for the number of mRNA or protein molecules, both as time-dependent and stationary distributions [1, 3, 15, 16, 21, 23, 24]. In fact, the Poisson distribution arises from a Poisson process, where the probability of an event occurring during a short time interval $[t,t+dt]$ is independent of events in other time regions [25, 26, 27]. Now, consider a simple transcription process where a gene is always active, and an mRNA molecule $X$ is transcribed with a rate $\alpha$,

It is clear that the creation events are indeed independent events, forming a Poisson process. Therefore, if we start from zero molecule at $t=0$, then the number $n$ of mRNA molecules at any later time $t$ is exactly the same as the number of creation events up to that time-point, and consequently, the Poisson distribution is actually the distribution of the number of creation events during the time interval $[0,t]$. However, consider a model where a degradation

is also included. Now, the degradation is not a Poisson process, although sometimes it is erroneously written as such in the literature. Since a molecule that already got degraded cannot be degraded again, the probability of a degradation event happening during the short time period $[t,t+dt]$ definitely depends on how many degradation events happened in times earlier than $t$. However, the Poisson distributions still appear ubiquitously in this class of models as the distributions of the number of mRNA or protein molecules [1, 3, 15, 16, 21, 23, 24], making its origin mysterious. Therefore, I address the following question in this paper: Given that the Poisson distribution appears so ubiquitously in the stochastic dynamics of GRN as the distribution of molecule number of mRNA or protein, is this molecule number equal to the number of certain independent events that happened during a given time interval? As I will show, the answer to this question is affirmative. In fact, it is the number of events where mRNA molecules are created, which are destined to be survive until the end of a given interval. The answer is very simple once stated, almost on the verge of being trivial. The most probable reason that it has not been discussed in the literature is that it may have been counter-intuitive to consider a birth of a particle with a given fate. However, it is important to note that although the fate of a molecule is not determined at the time of its creation, the probability of its given fate at the end of the time interval, death or survival, is already determined at the time of its creation, and it is all that matters in defining a Poisson process.

Using this viewpoint of Poisson processes, I can not only rederive Poisson distribution for models of gene expression with time-dependent rates and a non-Markovian model with time delay, I can also derive analytic forms of general time-dependent distributions by combining the Poisson distribution with other distributions. In fact, time-dependent distributions for arbitrary initial distributions are obtained for the Markovian model by combining the Poisson distribution with the binomial distribution, and distributions for very general class of initial distributions are obtained for a non-Markovian model by combining the Poisson distribution with the multinomial distribution.

The remainder of the paper is organized according to the order of increasing complexity. In section II, I will briefly review Poisson processes, where I will emphasize the importance of inhomogeneous Poisson processes, consisting of independent events that do not necessarily follow identical distributions. In section III, I will consider molecule creations without degradations, a textbook example of a Poisson process. In section IV, I will describe the molecule degradation, which is ${\it not}$ a Poisson process. I show that a general time-dependent probability distribution can be expressed as a superposition of binomial distributions. In section V, I will describe molecule creations with degradations and the underlying Poisson process. I show that any time-dependent distribution can be expressed in terms of the Poisson and the binomial distributions. The result will be generalized to a non-Markovian model with delayed degradation in section VI, where I will not only rederive the Poisson distribution in an earlier work with less effort, but also derive a general time-dependent distribution by combining the Poisson distribution with the multinomial distribution. Finally, the discussions are given in section VII.

The Poisson distribution in Eq.(1) is specified by a single parameter $\mu$, the expected number of independent events in a given time duration. If we partition the time interval $[0,t]$ into $N$ sub-intervals of small size $\Delta t\equiv t/N$, then the probability of events happening more than once in each sub-interval is $O\left((\Delta t)^{2}\right)$, and the probability of single occurrence of the event takes the form $p=\lambda\Delta t+O\left((\Delta t)^{2}\right)$ when such a probability is time-independent. Therefore, the total number of events in $[0,t]$ approximately follows the binomial distribution¹¹1The range of $n$ in the probability distribution $P(n,t)$ will be considered to be all the integers without any restriction in this work, unless specified otherwise. This is acceptable as long as we set $P(n,t)=0$ for illegitimate values of $n$. All the explicit forms of probability distributions considered here such as binomial, multinomial, and Poisson distributions contain factorials of negative integers in the denominators whenever the value of the molecule number is out of legitimate range, and vanish due to the fact that $(j!)^{-1}=[\Gamma(j+1)]^{-1}=0$ whenever $j$ is a negative integer. Taking the range of the molecule number to be the whole integer is especially convenient for summation, where the summation index can be shifted freely, which I will do throughout the text without further elaboration.,

the distribution of the number of successes in $N$ independent and identical trials, with success probability of $p$ at each trial. The Poisson distribution is recovered by taking the limit of $N\to\infty$ with $\mu(t)=Np=\lambda t$ fixed (Appendix A):

where $\mu(t)$ is the expected number of events in the time-interval $[0,t]$.

An important point to note here is that the condition of identical trials is not actually required to obtain the Poisson distribution. We only have to require independent trials [25, 26]. In this case, the binomial distribution in Eq.(4) is generalized to

where the probability of success at $i$-th trial is $p_{i}$, and the summation is over all distinct set of $n$ indices $\{i_{1}<i_{2},\cdots<i_{n}\}\subset\{1,2,\cdots N\}$. Again, the distribution in Eq.(6) becomes the Poisson distribution in the limit of $\Delta t\to 0$, now with the time-dependent function $\lambda(t)$, so that $\mu(t)=\int_{0}^{t}\lambda(t^{\prime})$ [26](Appendix A). This is the process where the probability of an event happening in an infinitesimal time interval $[t,t+dt]$ is $\lambda(t)dt$, called an inhomogeneous Poisson process, to distinguish it from the case with constant $\lambda$, called a homogeneous Poisson process.

This is a simple birth process

where $X$ describes an mRNA or protein molecule, and the creation rate $\alpha(t)$ is time-dependent in general²²2All the time-dependent rates in this work are the rates with predetermined time-dependences, and the corresponding model should not be confused with the one where the rates themselves are stochastic variables [28, 21, 23, 24].. Since molecule creations at distinct time points are independent of each other, these creation events form a Poisson process. If the number of molecules is zero at $t=0$, then the molecule number $n$ at a later time $t>0$ is equal to the number of creation events in the time interval $[0,t]$. Therefore, the probability distribution of $n$ is given by the Poisson distribution, $P(n,t)=P_{\rm Poisson}(n;\mu(t))$, with $\mu(t)=\int^{t}_{0}\alpha(t^{\prime})dt^{\prime}$.

Even when the number of molecules takes a non-zero value $n_{0}$ at $t=0$, the number of molecules $n$ at a later time $t>0$ is simply $n_{0}+n^{\prime}$ where $n^{\prime}$ is the number of creation events in the time interval $[0,t]$. Therefore, $n^{\prime}=n-n_{0}$ still follows the Poisson distribution, leading to a shifted Poisson distribution for $n$,

The dynamics of the probability distribution is Markovian, described by the master equation

It is straightforward to check that Eq.(8) is a solution of the master equation (9) by direct substitution. By the linearity of the master equation (9), an analytic form of the general solution with an arbitrary initial distribution $P(n,0)=v(n)$ can be constructed by superposing the expressions in Eq.(8),

From Eq.(10), we see that the master equation (9) admits stationary solutions

that depend on initial distributions, if and only if $\mu(\infty)=\int_{0}^{\infty}\alpha(s)ds$ exists. In particular, for a constant rate $\alpha$, there is no stationary solution because $\mu(t)=\alpha t$ increases indefinitely with time, and $P(n,t)$ converges to zero pointwise for all $n$. In this case, all the states in the Markov chain are transient states [29].

This is a simple death process

In contrast to the creation process in the previous section, this process is not a Poisson process, because any molecule that has been degraded cannot be degraded again. Therefore, the degradation events are not independent, and the probability of a degradation event depends on the number of molecules available. The master equation describing the dynamics of this model is

In order to get a non-zero number of molecules, the initial number $n_{0}$ of molecules must be non-zero. The number of molecules at later times can be obtained by using the fact that the event of a given molecule surviving at a later time $t$ is independent of what happens to other molecules, whose probability $p_{\rm surv}(t,0)$ is given by

Therefore, the number $n$ of surviving molecules at time $t$ follows the binomial distribution

It can be checked that the expression in Eq.(15) is the solution of the master equation (13) by direct substitution. Again, the general solution with an arbitrary initial distribution $P(n,0)=v(n)$ can be constructed by the superposition

A Poisson distribution appears in the special case where the initial distribution is Poissonian:

in which case we obtain the Poisson distribution $P(n,t)=P_{\rm Poisson}(n;\mu(t))$ with

due to the fact that the superposition of binomial distributions weighted by a Poisson distribution is again a Poisson distribution [26](Appendix B). In the next section, I will show that this distribution can be interpreted as a result of a Poisson process.

When the integral $\int_{0}^{t}\beta(s)ds$ diverges as $t\to\infty$, as in the case of a constant $\beta$, we get $\lim_{t\to\infty}p_{\rm surv}(t,0)=0$. In this case, we see from Eq.(16) that the probability distribution approaches the stationary distribution

regardless of the initial distribution.

This model is described by

which belongs to a class of models called the birth-death models [25, 26, 27]. Let us first consider the case where the number of $X$ molecules is initially zero. Molecules are created in the time interval $[0,t]$, but only a fraction of them survive at $t$, which I will call surviving molecules. When a molecule is created in the short time interval $[t^{\prime},t^{\prime}+dt^{\prime}]$ with $t^{\prime}<t$, its fate at $t$ is undetermined, but the probability of its survival $p_{\rm surv}(t,t^{\prime})$ at time $t$ is already determined to be

and consequently, the probability that a surviving molecule is created in $[t^{\prime},t^{\prime}+dt^{\prime}]$ is given by $\alpha(t^{\prime})p_{\rm surv}(t,t^{\prime})dt^{\prime}$, independent of events that happen in other regions of time. Therefore, the number of molecules at a later time $t$ is equal to the number these independent events in $[0,t]$, following the Poisson distribution $P_{\rm Poisson}(n;\mu(t))$ with

The master equation for this model is

and it is straightforward to check by direct substitution that $P_{\rm Poisson}(n;\mu(t))$ with $\mu(t)$ given by Eq.(22) is a solution of this equation.

We noted in the previous section that the distribution of the molecule number in the degradation-only model is Poissonian for $t>0$ if the initial distribution at $t=0$ is also Poissonian. In fact, this distribution can be considered as a special case of the Poisson distribution in the model with a time-dependent creation rate $\alpha(t)$, by shifting the origin of time: We start from zero molecule at some time $t_{0}<0$ so that we obtain a Poisson distribution for $P(n,0)$. We then require that $\alpha(t)=0$ for $t>0$ so that the model reduces to the degradation-only model for $t>0$.

The probability distribution converges to a stationary Poisson distribution if and only if the expected number of surviving molecules at $t\to\infty$, given by the integral

is finite. For example, if the rates $\alpha$ and $\beta$ are constants, we get

and

Now consider a more general situation where the initial number of molecules is $n_{0}$. The number of molecules at a later time $t$ is the sum of the number $n_{1}$ of surviving molecules among initial $n_{0}$ particles, whose distribution follows the binomial distribution in Eq.(15), and the number $n_{2}$ of surviving molecules created during the interval $[0,t]$, which follows the Poisson distribution. Therefore, the probability distribution for the molecule number $n$ is

where $p_{\rm surv}(t,0)$ and $\mu(t)$ are given by Eq.(14) and Eq.(22), respectively. It is straightforward to check that the solution (27) is a solution of the master equation (23) by direct substitution (Appendix C). By the linearity of the master equation, the solution for an arbitrary initial distribution $P(n,0)=v(n)$ can be constructed by superposing the distributions in Eq.(27),

Due to a complex mechanism of protein degradation, there can be a time-delay in protein degradation. Model with time-delayed degradation attempts to capture such a behavior, and once a molecule enters the degradation process, it gets degraded after a fixed time delay $\tau$ [15]. A more general model where the molecule is allowed to get degraded before $\tau$ was also constructed [16], which I will examine in this section in detail. The model is

Here, $X_{I}$ is an inactive molecule that has entered the degradation process. An active molecule is denoted by $X_{A}$, which can undergo instantaneous degradation with rate $\gamma$. It can enter delayed degradation process with rate $\beta$, at which point it becomes inactive and get degraded after time $\tau$ with certainty, but it can also undergo instantaneous degradation with rate $\zeta$ before the delayed degradation process is complete. Although it is straightforward to allow for time-dependent rates, I will keep them time-independent in order to compare the result with that of ref. [16], as well as for notational simplicity.

This model is non-Markovian, and the master equation for the probability distribution $P(n_{A},n_{I},t)$ takes the form

where $E^{k}_{A,I}f(n_{A,I})\equiv f(n_{A,I}+k)$, $P(n_{A},t)\equiv\sum_{n_{I}}P(n_{A},n_{I},t)$ is the marginal probability for $n_{A}$, and $P^{*}(n_{A},n_{I},t|n^{\prime}_{A}-1)$ is the probability under the initial condition of $n^{\prime}_{A}-1$ active molecules and no inactive molecules, obtained by neglecting the time-delayed degradation [16]:

Eq. (30) along with Eq. (31) admits the Poisson distribution as a time-dependent solution if the particle number is initially zero or the initial distribution is Poissonian [16]. However, it is easy to derive the Poisson distribution without going through the complicated process of solving Eqs.(30) and (31), by considering the underlying Poisson process.

Note that one of three things happen in a given short time interval of $[t^{\prime},t^{\prime}+dt^{\prime}]$: either a molecule that will survive at $t$ as an active one is created, the one that will be survive as an inactive one is created³³3This should not be confused with a birth of an inactive molecule. Every molecule is born active, and only an active molecule can turn into an inactive one. Here I am considering a birth of an active molecule that will eventually survive as an inactive molecule at the later time $t$, or none of these happens. These events are independent of creation events in other regions in time. Also, although the creation of a surviving active molecule and a surviving inactive molecule are exclusive events, they can be treated as independent events because they are rare events: even if we assume they are independent, probability of more than two molecules created in $[t^{\prime},t^{\prime}+dt^{\prime}]$ is negligible. Therefore, exclusive and independent events are indistinguishable if they are rare (Appendix D), and the probability distribution for the numbers $n_{A}$, $n_{I}$ of $X_{A}$, $X_{I}$ are given by the product of Poisson distributions,

It only remains to compute $\mu_{A}(t)$ and $\mu_{I}(t)$.

Suppose that a molecule is created at time $t^{\prime}$. In order for the molecule to remain active at a later time $t$, (i) it should not undergo instantaneous degradation ($X_{A}\xrightarrow{\gamma}\varnothing$), and (ii) should not turn into an inactive molecule $(X_{A}\xrightarrow{\beta}X_{I})$, during the intervening time. Therefore, the conditional probability $p_{A}(t-t^{\prime})$ that this molecule will remain active at $t$ is given by the exponential distribution in $t-t^{\prime}$,

where $a\equiv\beta+\gamma$, and we get

For a molecule created at time $t^{\prime}$ to survive at a later $t$ as an inactive molecule, (i) it must turn into an inactive one within a short time interval $[t_{I},t_{I}+dt_{I}]$ such that $\max(t^{\prime},t-\tau)<t_{I}<t$, (ii) it should not undergo instantaneous degradation ($X_{A}\xrightarrow{\gamma}\varnothing$) in the time between $t^{\prime}$ and $t_{I}$, and (iii) it should not undergo instantaneous degradation ($X_{I}\xrightarrow{\zeta}\varnothing$) in the time between $t_{I}$ and $t$. For a given value of $t_{I}$, such a conditional probability $\rho_{I}(t,t_{I}/t^{\prime})dt_{I}$ is given as

Integrating over $t_{I}$, we get the conditional probability $p_{I}(t-t^{\prime})$ that a molecule created at time $t^{\prime}$ will survive as an inactive one at $t$:

The expected number of surviving inactive molecules, $\mu_{I}(t)$, is then obtained by multiplying $p_{I}(t-t^{\prime})$ by the creation rate $\alpha$ and integrating over $t^{\prime}$:

Eq.(32) along with Eq.(34) and Eq.(41) completely specify the time-dependent Poisson distribution when we start from zero molecule at $t=0$, which agrees with the result in ref. [16]. The total number of molecules, $n=n_{A}+n_{I}$, also follows a Poisson distribution with the expectation value $\mu(t)=\mu_{A}(t)+\mu_{I}(t)$, due to the fact that the sum of two variables that follow Poisson distributions also follows a Poisson distribution [27]. The Poisson process underlying the Poisson distribution for $n=n_{A}+n_{I}$ is the creation of particles destined to survive until time $t$, regardless of being active or inactive.

Finally, as in the case of the Markovian model in the previous section, we can derive general time-dependent distributions by combining the Poisson distribution with the multinomial distribution. However, note that it is almost impossible to allow arbitrary initial distribution at $t=0$. If a non-zero initial value of $n_{I}$ is allowed, there is an ambiguity in the evolution of the system for $t>0$ because the final degradation of the initial inactive molecules depend on the exact time points at $t<0$ that they got inactivated. Therefore, I will only consider the case where only active particles are present at $t=0$, whose number is $n^{0}$. Given an active molecule at $t=0$, the probability that it will survive as an active molecule and the probability that it will survive as an inactive molecule, at a later time $t$, are given by Eqs.(33) and (38), respectively, with $t^{\prime}=0$:

Therefore, the probability that there are $n^{\prime}_{A}$ surviving active molecules and $n^{\prime}_{I}$ surviving active molecules among initial $n_{0}$ active molecules, at time $t$, is given by the multinomial distribution,

The numbers of active and inactive molecules at $t$ are decomposed as $n_{A}=n^{\prime}_{A}+n^{\prime\prime}_{A}$, $n_{I}=n^{\prime}_{I}+n^{\prime\prime}_{I}$, where $n^{\prime\prime}_{A}$and $n^{\prime\prime}_{I}$ are the numbers of the surviving active and inactive molecules created in the time interval $[0,t]$. Therefore, the probability distribution for $n_{A}$ and $n_{I}$ at $t$ is given by the combination of the multinomial and the Poisson distribution:

We can check that the distribution in Eq.(46) is indeed a solution of the master equation Eq.(30) by direct substitution (Appendix E). Again, by the linearity of the master equation, the time-dependent probability distribution for an arbitrary initial distribution of active molecules,

is obtained by the superposition of Eq.(46) weighted by $v(n_{0})$,

Poisson distributions appear ubiquitously in stochastic dynamics of gene expression, and the Poisson noise is considered to be the most basic type of noise when analyzing various components of stochastic fluctuations. However, when gene products are allowed to get degraded, it has not been clear whether the molecule number following a Poisson distribution is equal to the number of certain independent events in time, and if this is the case, what is the corresponding events. I answered this question in this work, by showing that the number of molecules distributed according to the Poisson distribution is actually equal to the number of creations of the molecules that are destined to survive until the end of a given time period, which are indeed independent events in time that form an inhomogeneous Poisson process. Using this viewpoint, I could derive the Poisson distribution not only for the Markovian model with time-dependent-rates, but also for the model with time-delayed degradation, without performing the difficult task of solving the non-Markovian master equation. Furthermore, I could obtain the time-dependent distribution for an arbitrary initial distribution in the case of the Markovian model, by combining the Poisson distribution with the binomial distribution, and the time-dependent distribution for an arbitrary initial distribution of active molecules in the case of the non-Markovian model, by combining the Poisson distribution with the multinomial distribution.

Of course, the molecule number follows the Poisson distribution only under the simplifying assumption of independent creation events. In the case of the protein, the number of protein molecules follows the Poisson distribution only if we approximate the transcription and the translation as a one-step process. In reality, the mRNA molecule has a non-zero life-time, during which protein gets translated with a certain rate, leading to bursty translations [8, 9, 10]. The creations of protein molecules are not a Poisson process in such a model, since the translation rate depends on the mRNA concentration. The Poissonian description also breaks down in most of the models with stochastic rates [28, 21, 23, 24]. In this class of models, the creation rate itself is a random variable distributed according to a probability distribution. A model with a stochastic creation rate can be used as an approximate model to describe the regulation of the expression by the transcription factor binding [28, 23], or to emulate the extrinsic noise due to heterogeneous cellular environments [23]. When the creation rate is a stochastic variable, the creations of the molecules form a Poisson process only if the creation rate has no temporal correlations, which is not valid for most of the non-trivial models. Despite these limitations, general analytic solutions found in this work may be of value as a basis of perturbations to obtain more realistic descriptions of the gene-regulatory network. For example, one may numerically solve the master equation for the stochastic rate, and then use the analytic solution for a given realization of the creation rate, so that the expectation values of various physical quantities are expressed as weighted averages over numerical distributions of the stochastic rates. Similar analyses may be done for other sophisticated models of gene regulatory network by combining the analytic solutions found in the current work with other analytic solutions and/or numerical computations.