Stochastic modeling of gene expression: application of ensembles of trajectories

TL;DR

This paper demonstrates how ensembles of trajectories can be used to analyze stochastic gene expression models, providing insights into dynamics, fluctuations, and rare events beyond static distributions.

Contribution

It applies the ensemble of trajectories approach to various gene expression models, including complex ones with bursts and switches, to analyze dynamics and rare events.

Findings

Generation of cumulant generating functions for protein numbers

Analysis of rare event statistics in gene expression models

Application to simple and complex stochastic models

Abstract

It is well established that gene expression can be modeled as a Markovian stochastic process and hence proper observables might be subjected to large fluctuations and rare events. Since dynamics is often more than statics, one can work with ensembles of trajectories for long but fixed times, instead of states or configurations, to study dynamics of these Markovian stochastic processes and glean more information. In this paper we aim to show that the concept of ensemble of trajectories can be applied to a variety of stochastic models of gene expression ranging from a simple birth-death process to a more sophisticate model containing burst and switch. By considering the protein numbers as a relevant dynamical observable, apart from asymptotic behavior of remote tails of probability distribution, generating function for the cumulants of this observable can also be obtained. We discuss the unconditional stochastic Markov processes which generate the statistics of rare events in these models.

Equations

$$\frac{d}{dt}P(c,t)=\sum _{c^{\prime }\ne c}{w}_{c^{\prime }\to c}P(c^{\prime },t)-\sum _{c^{\prime }\ne c}{w}_{c\to c^{\prime }}P(c,t)$$

$$⟨c|P(t)⟩\equiv P(c,t)$$

$$⟨c|H|c^{\prime }⟩\equiv w_{c^{\prime }\to c}\text{mbox}for{c}^{\prime }\ne c\text{mbox}and⟨c|H|c⟩\equiv -\sum _{c^{\prime }\ne c}{w}_{c\to c^{\prime }}\text{mbox}for{c}^{\prime }=c$$

$$\frac{d}{dt}|P(t)⟩=H|P(t)⟩$$

$$⟨1|H=0$$

$$H|P^{\ast }⟩=0.$$

$$\sum _i(t_{i+1}-t_i){\varphi }_{c(t_i)}$$

$$\sum _i{\theta }_{c(t_i)\to c(t_{i+1})}$$

$$⟨\mathcal{O}{⟩}_s=\frac{1}{Z}\sum _{\mathcal{C}}{\mathcal{O}}_t[\mathcal{C}]P[\mathcal{C}]e^{-s{\mathcal{A}}_t[\mathcal{C}]}.$$

$$Z=⟨e^{-ts{a}_t}{⟩}_s=⟨1|e^{tH(s)}|P(0)⟩$$

$$⟨c|H(s)|c^{\prime }⟩=w_{c^{\prime }\to c}{e}^{-s{\theta }_{c^{\prime }\to c}}\text{mbox}for{c}^{\prime }\ne c\text{mbox}and⟨c|H(s)|c⟩=-\sum _{c^{\prime }\ne c}{w}_{c\to c^{\prime }}-s{\varphi }_c\text{mbox}for{c}^{\prime }=c$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}\ln ⟨e^{-ts{a}_t}{⟩}_s={\Lambda }^{\ast }(s)$$

$$\lim_{t\to\infty}\frac{1}{t}\langle{\cal A}_{t}\rangle_{c}=(-1)^{n}\frac{d^{n}}{ds^{n}}\Lambda^{\ast}(s)\Big{|}_{s=0}$$

$$P(a_t\in da)\asymp e^{-tI(a)}da$$

$$I(a)=-\underset{s\in \mathbb{R}}{\inf }\{sa+{\Lambda }^{\ast }(s)\}$$

$$a=⟨a_t{⟩}_s=-\frac{d}{ds}{\Lambda }^{\ast }(s)$$

$$H(s)|{\Lambda }^{\ast }⟩={\Lambda }^{\ast }(s)|{\Lambda }^{\ast }⟩,$$

$$⟨{\tilde{\Lambda }}^{\ast }|H(s)={\Lambda }^{\ast }(s)⟨{\tilde{\Lambda }}^{\ast }|$$

$$|{\Lambda }^{\ast }⟩\propto \sum _c{P}_f(c,s)|c⟩,$$

$$|{\tilde{\Lambda }}^{\ast }⟩={⟨{\tilde{\Lambda }}^{\ast }|}^T\propto \sum _c{P}_i(c,s)|c⟩$$

$$P_f(c,s)=\frac{⟨c|{\Lambda }^{\ast }⟩}{⟨1|{\Lambda }^{\ast }⟩},$$

$$P_i(c,s)=\frac{⟨{\tilde{\Lambda }}^{\ast }|c⟩⟨c|P(0)⟩}{⟨{\tilde{\Lambda }}^{\ast }|P(0)⟩}.$$

$$P_{\text{mbox}TTI}(c,s)=\frac{⟨{\tilde{\Lambda }}^{\ast }|c⟩⟨c|{\Lambda }^{\ast }⟩}{⟨{\tilde{\Lambda }}^{\ast }|{\Lambda }^{\ast }⟩}.$$

$$H_{\text{mbox}eff}(s)=U(s)H(s)U^{-1}(s)-{\Lambda }^{\ast }(s)$$

$$⟨c^{\prime }|H_{\text{mbox}eff}(s)|c⟩=⟨c^{\prime }|H(s)|c⟩\frac{⟨{\tilde{\Lambda }}^{\ast }|c^{\prime }⟩}{⟨{\tilde{\Lambda }}^{\ast }|c⟩}.$$

$$\frac{d}{dt}P(n,t)=kP(n-1,t)+\mu (n+1)P(n+1,t)-(k+n\mu )P(n,t)$$

$$P^{\ast }(n)=e^{-\gamma }\frac{{\gamma }^n}{n!}\text{mbox}withn=0,1,2,3,\dots$$

$$⟨n⟩=\sum _{n=0}^{\infty }n{P}^{\ast }(n)=\gamma .$$

$$H(s)=\left(\begin{array}[]{ccccccc}-k&\mu&0&0&\cdots\\ k&-k-\mu-s&2\mu&0&\cdots\\ 0&k&-k-2\mu-2s&3\mu&\cdots\\ 0&0&k&-k-3\mu-3s&\cdots\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)\,.$$

$$H(s)|\Lambda ⟩=\Lambda (s)|\Lambda ⟩,$$

$$⟨\tilde{\Lambda }|H(s)=\Lambda (s)⟨\tilde{\Lambda }|$$

$$|{\Lambda }_j⟩=\sum _{n=0}^{\infty }{\varphi }_n^j|n⟩,$$

Full text

1 Introduction

Traditionally, stochastic simulation has been a powerful tool for studying the dynamics of gene regulatory networks [1]. However, this approach is not always efficient specially when rare events occur [2].

The ergodic theorem justifies the replacement of a long-time average of any physical quantity, or equivalently a laboratory measurement, with an ensemble average of that quantity overt configurations. However, when it the comes to the dynamics of the system when a rare value of a time-integrated quantity (time-dependent observables) is observed, the concept of the ensemble of trajectories can serve well. In this case one is interested in studying rare trajectories during which the system behaves in an unusual fashion over an extended period of time and time-averaged quantities remain far from their typical values. While it is hard to observe such rare events experimentally, because of having an exponentially small probability of occurrence, they can be studied using the theory of large deviations.

The theory of large deviations deals with the probabilities of rare events [3]. This theory describes how the probability of occurrence of a time-integrated dynamical observable which is largely deviated from its typical value decays to zero as the length of trajectory goes to infinity. On the other hand, the dynamics of stochastic systems is usually richer than what can be gleaned from their stationary properties. In order to fully understand dynamical behavior of the system, one necessarily has to focus on trajectories, and not just on configurations. A possible efficient approach for achieving that is to study the statistics of dynamical observables whose fluctuation behavior characterizes the dynamics of the system.

The time-integrated/averaged quantities and their fluctuations are of biological significance [4]. Given the fact that gene expression is intrinsically a stochastic process, fluctuations or rare events of proper time-integrated quantities might be observed in this process. Several prominent examples of such rare events have been observed in biological systems which include phenotypic transitions [5]. Cumulative sum of proteins have also been considered in modeling approaches dealing with cancer therapy [6]; therefore, fluctuations of such qualities could be interesting from biological point of view.

Recently, a vast amount of research is devoted to show that the dynamics of stochastic systems can be studied by analyzing the statistical properties of dynamical trajectories (time-realizations) [7, 8, 9, 10, 11, 12, 13, 14, 15]. It is of great importance to characterize fluctuations in the system conditioned on the occurrence of a rare event. Ensembles of trajectories are associated with large deviations of time-integrated quantities. While some trajectories are responsible for creating typical values of the observable, some other trajectories create a specific fluctuation or an atypical value of the observable. By looking at a certain set of trajectories which is responsible for a specific fluctuation, we are imposing constraints on the time-integrated quantities [16]. It has been shown that one can determine the changes in dynamical model parameters so that it reproduces the effects of rare fluctuations [17]. Far from the boundaries of the observation time, there is a temporal regime during which the trajectories of the biased ensemble are close to the steady-state trajectories of an auxiliary model (sometimes called driven or effective model) which is a conditioning-free stochastic Markov process. In other words, the statistics of the conditioning-free Markov stochastic process reproduces the fluctuations of the original Markov process conditioned on the occurrence of a rare event [18].

In comparison with the static ensemble of configurations approach, the above mentioned approach provides us, in principle, with a variety of information. Assuming that the probability distribution function of the dynamical observable has a large deviation form, one can calculate the large deviation function. The scaled cumulant generating function of the time-integrated observable can also be calculated. It is also possible to determine the probability vectors of both the final and initial configuration, knowing that the value of the observable through the evolution of the system is restricted to a certain value. Last but not least, the effective process can also be calculated [3, 15, 17].

In this paper we aim to show that using the ensemble of trajectories approach one can obtain a better understanding of the dynamical behavior of the stochastic models of gene expression. After a brief review of Mathematical preliminaries we start with the simplest model of gene expression consisting of Poissonian creation and degradation of proteins. We will then add burst and switch of promoter to the models and show that these models are still solvable. In comparison with a recent work on modeling of stochastic gene expression conditioned on large deviations [19], we will show that it is possible to add degradation of proteins/mRNA’s to their model and, apart from what have been calculated, one can answer quite interesting questions such as which dynamical trajectories have the most contribution to occurrence of a rare event.

2 Mathematical preliminaries

Let us consider a stochastic process system with discrete configuration space $\{c\}$. The time evolution equation for probability of being in configuration $c$ at time $t$, denoted by $P(c,t)$, is given by

[TABLE]

in which $w_{c\to c^{\prime}}$ is a time-independent rate of jumping from configuration $c$ to $c^{\prime}$. By introducing an orthonormal complete basis $\{|c\rangle\}$ with $\langle c|c^{\prime}\rangle=\delta_{c,c^{\prime}}$ and defining

[TABLE]

and the matrix elements of $H$ as

[TABLE]

we can rewrite the Master equation (1) as follows [20]

[TABLE]

in which $|P(t)\rangle$ is the probability vector at time $t$ and $H$ is the stochastic generator of the process with the following property

[TABLE]

where we have defined a summation vector $\langle 1|$ as $\langle 1|\equiv\sum_{c}\langle c|$. In the long time limit the steady-state probability vector satisfies

[TABLE]

Instead of looking at the time-evolution of the probability vector and calculating the average values of the observables over a static ensemble of configurations, we can look at realizations or trajectories of the process in the space of configurations of the system. By defining dynamical observables over these trajectories, we can then study the dynamical properties of the process such as dynamical phase transitions [16, 21]. There are usually two types of dynamical observables or equivalently time-integrated observables. The first type is a purely spatial observable and depends on both the configurations that the system meets along a trajectory and the time it spends in each configuration. This type of a time-integrated observable which is a functional of the trajectory can be written as

[TABLE]

in which $\phi_{c(t_{i})}$ is the increment of the observable when it meets the configuration $c$ at time $t_{i}$. For example, $\phi_{c}(t_{i})$ can be density of particles or energy of the system associated with configuration $c$ at time $t_{i}$. The second type of a time-integrated observable depends on the transitions among consecutive configurations along a trajectory and can be written as

[TABLE]

in which the increment of this current is denoted by $\theta_{c(t_{i})\to c(t_{i+1})}$. This kind of dynamical observable is sometimes called current-like observable, such as particle current, entropy production rate and activity. For example, taking $\theta_{c(t_{i})\to c(t_{i+1})}=1$ for all transitions led to counting the number of configuration changes during the observation time which is called the activity of a system [3, 17]. We can build an ensemble of trajectories in different ways as we do in traditional classical statistical mechanics for the ensemble of configurations. This will be discussed in the following.

We start with the micro-canonical ensemble of configurations defined for an isolated system with fixed energy $E$. The constraint of fixed energy means that every member of the micro-canonical ensemble has the same energy. It is exactly because of this constraint that calculating the conditional probability of being in $c$ given that $E$ is fixed $P(c|E)$ might be difficult. Similarly, but not on the same footing, we can consider an ensemble of trajectories, defined between an initial time [math] and a final time $t$ (called the observation time which is assumed to be very long), with a constraint on the value of a given dynamical observable over the trajectories. This ensemble contains those trajectories, among all possible trajectories, for which the value of that observable is fixed along the trajectory [15, 16, 22].

We can also consider an ensemble of trajectories assuming that the average value of a given dynamical observable is constant during the observation time [21, 22]. This is similar to the canonical ensemble of configurations in the traditional statistical mechanics where the probability of being in a configuration is calculated at fixed temperature (i.e. the average value of the energy is fixed). This can be achieved by multiplying the probabilities in the micro-canonical ensemble by $e^{-\beta H(c)}$ in which $H(c)$ is the Hamiltonian of the system while being in the configuration $c$. As the temperature is varied, the average energy is changed. Calculating the average value of any observables in this canonical ensemble of energy gives the typical value of those observables under fixed temperature. For the ensemble of trajectories this can be formulated as follows: We multiply the probability of taking the trajectory $\cal C$, denoted by $P[{\cal C}]$, by $e^{-s{\cal A}_{t}[{\cal C}]}$ in which ${\cal A}_{t}$ is a time-extensive dynamical observable explained before. $s$ is a Mathematical biasing field conjugated to $\mathcal{A}_{t}$, analogous to parameter $\beta$ in traditional statistical physics. The ensemble average of a dynamical observable, say ${\cal O}_{t}[{\cal C}]$, which is a functional of $\cal C$ is given by

[TABLE]

The dynamical partition function $Z$ appeared in (3) is a normalization factor whose logarithm is called the dynamical free energy. Note that, as in the canonical ensemble, the field $s$ conjugated to ${\cal A}_{t}$ can fix the average value of this observable during the observation time. This means that the averages are now calculated in an ensemble of trajectories under fixed $s$. The above statistics is constructed over a biased ensemble of trajectories which is sometimes called the $s$-ensemble. Averages in the $s$-ensemble with $s=0$ correspond to the steady-state averages of ${\cal O}_{t}$ which is the only physically accessible ensemble. These averages are clearly time-translational invariant; however, this invariance is broken for the averages over the $s$-ensemble with $s\neq 0$. If the observation time is very long and we are far from the boundaries $t^{\prime}=0$ and $t^{\prime}=t$, there is a temporal regime during which the time translational invariance (TTI) is recovered. In other words, the $s$-ensemble averages during the TTI regime do not depend on the instance of time at which the averaging in performed. However, before the system relaxes into the TTI regime and after it goes out of it, these averages depend on the time of averaging [15, 17].

Let us look at the large deviations of the dynamical observable ${\cal A}_{t}$ in the long time limit. The generating function of the moments of the time-integrated quantity $a_{t}=\frac{1}{t}{\cal A}_{t}$ can be written as [15]

[TABLE]

in which the average $\langle e^{-tsa_{t}}\rangle_{s}$ is taken over all possible trajectories during the time interval $[0,t]$ and $H(s)$ is a modified generator which is equal to the stochastic generator of the process $H$ at $s=0$. The matrix elements of $H(s)$ are

[TABLE]

in which $\phi_{c}$ and $\theta_{c^{\prime}\to c}$ are the increments of the time-dependent observables of the first and the second type respectively from which we construct the $s$-ensemble [3, 16]. In the long time limit one can write

[TABLE]

where $\Lambda^{\ast}(s)$, which is called the Scaled Cumulant Generating Function (SCGF) of the observable, is the largest eigenvalue of the modified generator $H(s)$. Note that the derivatives of $\Lambda^{\ast}(s)$ evaluated at $s=0$ give the cumulants of ${\cal A}_{t}$ scaled by time i.e.

[TABLE]

from which one can characterize the fluctuations of the observable. According to the Gärtner–Ellis Theorem if $\Lambda^{\ast}(s)$ exists and is differentiable for all $s\in\mathbb{R}$, then $a_{t}$ satisfies a large deviation principle namely

[TABLE]

with the rate function $I(a)$ given by the Legendre–Fenchel transform of $\Lambda^{\ast}(s)$

[TABLE]

in which the symbol “inf” above stands for “infimum of”, which can be taken to mean the same as “minimum of”. Expanding $I(a)$ beyond quadratic order gives information about non- Gaussian fluctuations of ${\cal A}_{t}/t$, which are referred to as large deviations. Moreover, from (7) we find

[TABLE]

a relation which connects the average value of the observable $a_{t}$ to the conjugate field $s$ [3]. The right and the left eigenvectors of $H(s)$ corresponding to the largest eigenvalue of the modified generator $\Lambda^{\ast}(s)$ are quite meaningful in the context of large deviations in terms of trajectories. Considering the following eigenvalue relations

[TABLE]

it has been shown that

[TABLE]

in which $P_{f}(c,s)$ and $P_{i}(c,s)$ are in fact the probability of the final and initial configuration respectively, knowing that the value of $a$ is fixed and related to the observable $s$ through (8) [7, 15, 17]. It is worth mentioning that these normalized probabilities are defined as follows

[TABLE]

We could also look at the probability of observing a given configuration $c_{t^{\prime}}$ at a time $t^{\prime}$ during the evolution of the system far from the initial and final configuration, i.e. $0\ll t^{\prime}\ll t$, conditioned on fixed $s$. Starting from $t^{\prime}=0$ the system relaxes, on a time scale $\tau$, into the TTI regime during which the probability of being in the configuration $c$ at fixed $s$ is given by

[TABLE]

During the TTI regime defined as $\tau\ll t^{\prime},\tau\ll t-t^{\prime}$ in which $\tau$ is a relaxation time into this temporal regime, the steady-state trajectories of the effective process, discussed below, are those of the biased ensemble of trajectories [16, 17, 18].

As we mentioned before, in the long time limit each specific fluctuation in the system can be described by a stochastic Markov process called the effective process. This process is equivalent to the conditioning of the original process on seeing a certain fluctuation. The stochastic generator of this effective stochastic process is given by

[TABLE]

which is a generalization of Doob’s h-transform [16, 17]. $U(s)$ in (16) is a diagonal matrix with the matrix element $\langle c|U(s)|c\rangle=\langle\tilde{\Lambda}^{\ast}|c\rangle$. The off-diagonal matrix elements of (16) are given by

[TABLE]

and the diagonal elements can be obtained using the fact that (16) is a stochastic matrix [17].

Finally, there is an important family of stochastic processes for which the effective process is identical to the original process except that the effective rates are just rescaled values of those in the original process [23]. This property has an important consequence. If we can calculate the the effective rates (by calculating the left eigenvector $\langle\tilde{\Lambda}^{\ast}|$) and the steady-state probability vector $|P^{\ast}\rangle$, we can substitute the original rates in the steady-sate probability vector with the effective rates and using (15) obtain the right eigenvector $|\Lambda^{\ast}\rangle$ of the modified generator associated with $\Lambda^{\ast}(s)$. We will use this property later in the paper.

Considering this short introduction to the large deviations, we apply the above results to different well-known models of gene expression. We will start with the simplest model and then cover the more detailed ones. We will show how the powerful toolbox of large deviations can serve us find considerable amount of information about the dynamics of these models.

3 The models

As we mentioned earlier, gene expression is intrinsically a stochastic process. After modeling of gene expression as a Markovian stochastic process one can, in principle, calculate the probability distribution of many relevant observables, such as number of proteins or mRNA’s, and also their moments [24, 25]. However, until recently, not much is done on the dynamical properties of this stochastic process when these observables are considered as a dynamical observables.

The toolbox of large deviations, as explained briefly in the previous section, allows one to study the dynamical properties of the stochastic processes as well as the large deviations of these observables when considered as dynamical variables defined on trajectories or time-realizations of the process. This approach, as we will see, provides us with a deeper understanding of the dynamical properties of gene expression as a Markovian stochastic process.

In what follows, we will consider stochastic models of gene expression at mRNA or protein level and investigate the dynamical properties of these models using the concept of $s$-ensemble. We will show that the large deviation form of the probability distribution of dynamical observables and also the cumulants of these observables can be calculated exactly.

3.1 Model 1

As the first model we consider a birth-death model at protein level in which the proteins are created with the rate $k$ and degradation is occurred with the rate $n\mu$ in which $n$ is the number of the proteins [26]. A simple sketch of this model is given in Fig. (1). The Master equation for this process can be simply written as

[TABLE]

in which $P(n,t)$ is the probability of existing $n$ protein at time $t$. Note that here the configuration of the system is determined by the number of proteins and that the space of configurations is infinite dimensional. The steady-state probability distribution $P^{\ast}(n)$ (which turns to be an equilibrium steady-state in this case) can easily be obtained by letting $t\to\infty$ and solving (18). The result is given by

[TABLE]

in which $\gamma\equiv\frac{k}{\mu}$. The average number of proteins in the steady-state (which is sometimes called the typical value) can be easily calculated

[TABLE]

Since mean and noise of protein numbers are quite relevant quantities in the realm of gene expression, let us consider the number of proteins as a dynamical observable (the first type) defined on the trajectories or time-realizations of the process during an observation time of length $t$. More precisely, we consider the first type time-integrated quantity $\mathcal{A}_{t}$ as $N_{t}=\int_{0}^{t}n(t)dt$ and define the time-averaged number of protein over observation time $t$, i.e. $n_{t}=N_{t}/t$, as the dynamical observable. The modified generator $H(s)$ in the complete orthonormal basis $\{|0\rangle,|1\rangle,|2\rangle,|3\rangle,\dots\}$ is given by

[TABLE]

It turns out that the full spectrum of the modified generator (20) can be calculated exactly. Considering the following eigenvalue relations

[TABLE]

we have found

[TABLE]

and

[TABLE]

where

[TABLE]

$C_{n}(j;x)$ is the Poisson-Charlier polynomial defined explicitly by [27]

[TABLE]

for $x>0$. These polynomials are a family of orthogonal polynomials satisfying the following relation

[TABLE]

Note that in (23) $s$ is restricted to $s\in]-\mu,+\infty]$. On the other hand, $s\in]-\mu,0[$ ($s\in]0,+\infty]$) correspond to the ensemble average of the observable larger (smaller) than its typical value while $s=0$ gives the typical value of observable. Note that using $\sum_{j=0}^{\infty}\phi_{m}^{j}\tilde{\phi}_{n}^{j}=\delta_{m,n}$, the left and right eigenvectors of $H(s)$ satisfy

[TABLE]

The largest eigenvalue of (20) is given by $j=0$ in (23) i.e. $\Lambda^{\ast}(s)=\frac{-ks}{\mu+s}$. As we mentioned earlier, derivatives of $\Lambda^{\ast}(s)$ with respect to $s$ generate the cumulants of $n_{t}$ in the $s$-ensemble. The first derivative of the largest eigenvalue at $s=0$ gives

[TABLE]

According to (6), the higher-order derivatives of $\Lambda^{\ast}(s)$ at $s=0$ give higher-order cumulants of $N_{t}$ scaled with time $t$. It is worth mentioning that while the first cumulant of the observable in the steady-state is equal to the first cumulant of the dynamical observable at $s=0$, the higher-order cumulants obtained from steady-state distribution of the observable are not equal to the dynamical ones at $s=0$. The reason is that the definition of these observables are basically different.

The right and the left eigenvectors of (20) corresponding to the largest eigenvalue $\Lambda^{\ast}(s)$ (or $j=0$ in (23)) are now given by

[TABLE]

and

[TABLE]

respectively.

Let us now investigate the effective dynamics of the process. As we have already explained, an unconditional stochastic process for which the typical value of the dynamical observable is an atypical value of the same observable in the original process, is called the effective process. Considering the number of proteins as a dynamical observable, the matrix elements of the stochastic generator of this effective process are found to be

[TABLE]

We remind the reader that the configuration of the system $c$ is now defined by the number of proteins in the system $n$. As can be seen these effective rates of birth and death are just rescaled values of those in the original process.

It turns out that in the long time limit, the probability distribution for the number of proteins $n_{t}$ has a large deviation form $P(n_{t}=\bar{n})\asymp e^{-tI(\bar{n})}$ with the following rate function

[TABLE]

Let us have a look at another interesting quantity that is the $s$-ensemble average of protein numbers at time $t^{\prime}$, a time between [math] and $t$. This quantity is given by

[TABLE]

in which the modified generator $H(s)$ is given by (20) and $\tilde{P}(n,s,t^{\prime})\equiv\langle n|\tilde{P}(s,t^{\prime})\rangle$ is the solution of the following Master equation

[TABLE]

The diagonal matrix $\hat{n}$ has the diagonal elements $\langle n|\hat{n}|n\rangle=n$. Straightforward calculations result in

[TABLE]

in which

[TABLE]

Finally, after some calculations, we find

[TABLE]

It is easy to see that

[TABLE]

which comes from the fact that $s$-ensemble is time-reversal symmetric. The result (34) is also interesting from the point that in the long time limit the $s$-ensemble average of the number of proteins at the initial and final times are given by

[TABLE]

After a relaxation time $\tau=1/(\mu+s)$ the system falls into the TTI regime during which we have

[TABLE]

As we mentioned, this model is in equilibrium in the long time limit i.e. the steady-state probability distribution $P(n)$ given by (19), besides the birth and death rates, satisfy the detailed balance condition. On the other hand, the average of the current of the proteins in the system (the current of birth minus the current of death) is zero because of the detailed balance. However, in the $s$-ensemble of trajectories the average of this current is only zero during the TTI regime while it is non-zero during the initial and final transient regimes characterized by the relaxation time $\tau$. In order to clarify this, we have studied that average of the net current of proteins in the $s$-ensemble of number of proteins. The $s$-ensemble average of the net current is given by

[TABLE]

in which $J$ is a square matrix with the following elements

[TABLE]

After some straightforward calculations one finds

[TABLE]

At three different times, including the TTI regime, the asymptotic behavior of the current is given by

[TABLE]

As can be seen the ensemble average of the net current is non-zero during the transient regimes while it is zero during the TTI regime as expected. In Fig. 2 we have plotted both $\langle n(t^{\prime})\rangle_{s}$ and $\langle J(t^{\prime})\rangle_{s}$ as a function of $t^{\prime}$ for different values of $s$. As can be seen during the transient regimes both quantities are functions of time except $s=0$; however, there is a temporal regime where $\langle n(t^{\prime})\rangle_{s}$ is a constant while $\langle J(t^{\prime})\rangle_{s}$ is zero.

We conclude the analysis of this model in the $s$-ensemble of trajectories by discussing the probability of observing a given configuration $n$ at time $t^{\prime}\in[0,t]$ given that $s$ is fixed. This quantity is given by

[TABLE]

Straightforward calculations result in

[TABLE]

The fact that the $P_{s}(n,t^{\prime})=P_{s}(n,t-t^{\prime})$, which reveals the time-reversal symmetry of $s$-ensemble, can easily be seen. Now one can take the limit of this probability in the long time limit to find that

[TABLE]

These can also be obtained from the properties of the left and the right eigenvectors of $H(s)$ associated with its largest eigenvalue, as already explained in (13), (14) and (15). More specifically, $\lim_{t\to\infty}P_{s}(n,0\ll t^{\prime}\ll t)$ can be obtained by substituting the effective rates in the steady-state probability distribution function (19).

In Fig. 3 we have plotted (3.1) as a function of $t^{\prime}$ for three different values of $s$. The diagrams show which configuration is more probable in each case as $s$ is varied. As can be seen in (3.1), the probability distributions at the beginning and the end of the trajectory are the same because of time-reversal symmetry. The average value of the number of proteins during the TTI regime is given by the first derivative of $\Lambda^{\ast}(s)$ with respect to $s$.

3.2 Model 2

Recent experiments have observed a feature of gene regulation which can be captured by defining stochastic processes in which protein production often occurs in bursts resulting from many factors [25, 28, 29]. In what follows we consider an effective model for protein production. We assume that protein production occurs with rates $k$, and burst sizes $n$ drawn from a state-dependent geometric distribution $b_{n}$ given by

[TABLE]

The Master equation for the probability distribution $P(n,t)$ of having $n$ protein at time $t$ is

[TABLE]

The steady-state distribution is given by

[TABLE]

in which the symbol $(x)_{n}\equiv x(x+1)\cdots(x+n-1)$ is the ordinary Pochhammer symbol. The average number of proteins in the steady-state is also given by

[TABLE]

Let us now consider the number of proteins as a dynamical observable defines on the time-realizations of the process and investigate the fluctuations of this quantity. The modified generator for this process is given by

[TABLE]

It turns out that the largest eigenvalue of $H(s)$ is given by

[TABLE]

and the right and left eigenvectors corresponding to this eigenvalue are also as follows

[TABLE]

Using the left eigenvector $\langle\tilde{\Lambda}^{\ast}|$ the effective rates can be calculated using (17)

[TABLE]

for $n,r=0,1,2,\dots$. Considering the number of proteins as a dynamical observable, its probability distribution in the long time limit has a large deviation form which is given by

[TABLE]

in which the rate function $I(\bar{n})$ can easily be calculated using (7)

[TABLE]

3.3 Model 3

In the third example, we consider a stochastic process of gene expression from a promoter with $2$ internal states $i=0,1$. The promoter makes random transition from [math] to $1$ with rate $\alpha$ and from $1$ to [math] with rate $\beta$. In each state, a single mRNA is generated with rate $k_{0}$ ($k_{1}$) when the system is at the state [math] ($1$) [26]. This model is schematically drawn in Fig. (4)

In what follows we will show that it is possible to present a full description of the dynamical properties of the process. On the other hand, we can discuss the fluctuations of the number of proteins in the long time limit. Considering the number of proteins as the dynamical observable the modified generator of the process in the complete orthonormal basis $\{|i,n\rangle\}=\{|0,0\rangle,|1,0\rangle,|0,1\rangle,|1,1\rangle,\dots\}$ is given by

[TABLE]

in which

[TABLE]

for $n=0,1,2,\dots$ and also

[TABLE]

It turns out that the largest eigenvalue $\Lambda^{\ast}(s)$ and its corresponding left and right eigenvectors of $H(s)$ can be calculated exactly. Assuming the following form for the left eigenvector

[TABLE]

in which

[TABLE]

one obtains the largest eigenvalue of the modified generator

[TABLE]

As we mentioned earlier, the derivatives of $\Lambda^{\ast}(s)$ at $s=0$ generate the cumulants of the observable in the steady-state of the process. It is easy to check the the average of the proteins in this model is given by

[TABLE]

The right eigenvector corresponding to $\Lambda^{\ast}$ can also be calculated with more efforts. Let us consider

[TABLE]

Straightforward calculations show that the generating functions of $\phi^{(0)}_{n}$ and $\phi^{(1)}_{n}$ defined as

[TABLE]

satisfy

[TABLE]

We will show how to calculate $f_{0}(x)$, because $f_{1}(x)$ can easily be obtained by applying the following transformation to $f_{0}(x)$

[TABLE]

By defining a new variable $z$ as

[TABLE]

the equation governing $f_{0}(z)$ becomes

[TABLE]

in which

[TABLE]

Now by choosing

[TABLE]

in which

[TABLE]

we obtained the equation governing $\mathcal{F}_{0}(z)$ as follows

[TABLE]

in which

[TABLE]

The final result is now given by a Hypergeometric function as follows

[TABLE]

Having the left eigenvector associated with the largest eigenvalue, we can calculate the effective rates as follows

[TABLE]

for $n=0,1,2,3,\dots$. As can be seen, the effective rates are rescaled when compared to those in the original process.

3.4 Model 4

The fourth example is similar to the third model except we have added burst. As in the previous model we consider a stochastic process of gene expression from a promoter with $2$ internal states $i=0,1$. The promoter makes random transition from [math] to $1$ with rate $\alpha$ and from $1$ to [math] with rate $\beta$. In each state, bursts of gene expression leading to the production of mRNAs occur with rate $k_{0}$ ($k_{1}$) when the system is at the state [math] ($1$), and burst sizes $n$ drawn from a state-dependent distribution $b_{n}$ as introduced in (44). This model can also be used to represent gene expression at the level of proteins. As in the previous models we are going to show that the largest eigenvalue of the modified generator besides the left eigenvector associated with this eigenvalue can be calculated exactly. These quantities are enough to provide us with a full description of the dynamical properties of the process.

By considering the complete orthonormal basis $\{|i,n\rangle\}=\{|0,0\rangle,|1,0\rangle,|0,1\rangle,|1,1\rangle,\dots\}$ the modified generator of the process in the presence of burst is given by

[TABLE]

in which

[TABLE]

and also

[TABLE]

for $n=0,1,2,\dots$. It turns out that the largest eigenvalue and also the corresponding left eigenvector can be calculated as in the third model. Assuming the following form for the left eigenvector

[TABLE]

in which

[TABLE]

one obtains

[TABLE]

The first derivative of the largest eigenvalue at $s=0$ gives the average number of proteins in the steady-state

[TABLE]

Note that for $k_{0}=k_{1}=k$, we recover the results of the second model i.e.

[TABLE]

In Fig. 5 we have plotted the largest eigenvalue and the rate function calculated using the Legendre–Fenchel transformation (7). Note that as $\beta$ approaches to zero the first derivative of the largest eigenvalue changes discontinuously at certain $s$. This indicates a first-order dynamical phase transition.

In contrast to the previous cases, we have found that calculating the right eigenvector of $H(s)$ for an arbitrary $s$ is a formidable task. The effective rates of the effective dynamics can be calculated using (17) as for the previous models. The results are

[TABLE]

for $n,r=0,1,2,3,\dots$. As for the previous models, the effective rates of the effective dynamics are just rescaled in comparison to the the reaction rates of the original process.

3.4.1 Generalization

The previous model can be generalized to the one in which the promoter has more than two internal states. We assume that the promoter changes its internal state from $i$ to $j$ with the rate $\alpha_{i\to j}$ with $i,j=1,2,3,\dots,N$. The Master equation governing the time evolution of the probability distribution function for having $n$ proteins at time $t$ when the promoter is in the state $i$ is given by

[TABLE]

Considering the $s$-ensemble of number of proteins, the modified generator is given by

[TABLE]

in which the matrix elements of $N\times N$ matrices $A_{n}$, $B_{n}$ and $D_{n}$ are given by

[TABLE]

for $i,j=1,2,3,\dots,N$. It turns out that, because of the special structure of $H(s)$, the largest eigenvalue and its corresponding left eigenvalue of $H(s)$, as an infinite-dimensional matrix, reduce to those of an $N\times N$ matrix. Considering the following left eigenvector

[TABLE]

with

[TABLE]

we only need to solve

[TABLE]

in which $\tilde{H}(s)$ is defined as follows

[TABLE]

Now the effective rates can be calculated, after calculating $\phi_{n}$’s, as follows

[TABLE]

with $i,j=1,2,3,\dots,N$ and $n,r=1,2,3,\dots$.

Comparing with a recent work on stochastic gene expression conditioned on large deviations [19], this generalization adds degradation to the process while remaining integrable. As we saw, the conditioning-free effective process is represented by a process similar to the original process except that it has renormalized parameters [23].

4 Concluding remarks

In this paper, we studied four Markovian stochastic models of gene expression conditioned on large deviations of the population of proteins. Starting with a simple birth-death model, which is a continuous-time Markov process often used to study how the number of individuals in a population change through time, we showed that the model if fully integrable in the sense that the SCGF and also the corresponding left and right eigenvectors of the associated modified generator could be calculated exactly. The SCGF gives us, through the Gärtner–Ellis Theorem, the large deviation form of the probability distribution function of the number of proteins. Dynamical behavior of the system also obtained using the probabilistic concepts of the left and right eigenvectors associated with the SCGF.

We gradually generalized the simple birth-death process to the stochastic processes which were more realistic in the realm of gene expression. We showed that they were still integrable and fluctuations of the dynamical observable could be traced using the large deviation theory. Finally, we showed that the modeling of gene expression as a Batch Markovian Arrival Process (BMAP) studied in [19] could be generalized to a model which includes degradation and that this generalization does not alter the integrability of the system.

As it has been explained in [30], the SCGF is not always given by the largest eigenvalue of the modified generator. This might happen in the systems with infinite-dimensional configuration space. We have checked that, despite the configuration space is infinite-dimensional, the SCGF is still given by the largest eigenvalue of the modified generator.

Acknowledgment

We would like to thank Prof. Rahul Kulkarni for his suggestions through an email correspondence.