Analytical results for a stochastic model of gene expression with arbitrary partitioning of proteins

TL;DR

This paper provides an exact analytical solution for a stochastic gene expression model with arbitrary protein partitioning, utilizing the PPA-mapping to transform the problem into a hierarchy of simpler models, enabling explicit calculations of key statistics.

Contribution

It introduces a novel application of the PPA-mapping to a variation of the two-stage gene expression model with time-dependent rates, deriving explicit solutions and extending to models with parameter fluctuations.

Findings

Derived an integral expression for the time-dependent generating function.

Obtained explicit formulas for mean, variance, and correlation functions.

Extended results to models with stochastic parameter fluctuations.

Abstract

In biophysics, the search for analytical solutions of stochastic models of cellular processes is often a challenging task. In recent work on models of gene expression, it was shown that a mapping based on partitioning of Poisson arrivals (PPA-mapping) can lead to exact solutions for previously unsolved problems. While the approach can be used in general when the model involves Poisson processes corresponding to creation or degradation, current applications of the method and new results derived using it have been limited to date. In this paper, we present the exact solution of a variation of the two-stage model of gene expression (with time dependent transition rates) describing the arbitrary partitioning of proteins. The methodology proposed makes full use of the the PPA-mapping by transforming the original problem into a new process describing the evolution of three biological switches. Based on a succession of transformations, the method leads to a hierarchy of reduced models. We give an integral expression of the time dependent generating function as well as explicit results for the mean, variance, and correlation function. Finally, we discuss how results for time dependent parameters can be extended to the three-stage model and used to make inferences about models with parameter fluctuations induced by hidden stochastic variables.

Tables

Table 1. Table 1: Transitions and associated rates for the original model (model-0).

Event	Update	Transition rates
$𝒜$-production	$n\to n+1$	$k\left(t\right)$
$𝒜$-degradation	$n\to n-1$	$n\mu \left(t\right)$
${𝒫}_j$-production	$m_j\to m_j+1$	$n{q}_j\left(t\right)$
${𝒫}_j$-degradation	$m_j\to m_j-1$	$m_j{\gamma }_j\left(t\right)$

Table 2. Table 2: Three different types of histories need to be considered. Each path starts with θ = 0 𝜃 0 \theta=0 , and transits no more than once from OFF to ON.

Representation	Description
	$\theta$ is constantly in the OFF-state
	$\Theta \left(\tau \right)=0,\forall \tau \in [0,t[$
	$\theta$ is switching state at time $s$
	$\Theta \left(\tau \right)=1$, if $\tau \in [s,t[$ and $\Theta \left(\tau \right)=0$ otherwise
	$\theta$ is switching state at time $s$ and $s^{\prime }$
	$\Theta \left(\tau \right)=1$, if $\tau \in [s,s^{\prime }[$ and $\Theta \left(\tau \right)=0$ otherwise

Equations

$$\frac{d}{dt}{P}_{n,m_1,m_2}$$

$$\frac{d⟨n⟩}{dt}$$

$$\frac{d⟨m_j⟩}{dt}$$

$$\frac{⟨n⟩(t)}{⟨n{⟩}^{\ast }}$$

$$\frac{\langle m_{j}\rangle(t)}{\langle m_{j}\rangle^{*}}=\left\{\begin{array}[]{cc}1-\frac{1}{\gamma_{j}-\mu}(\gamma_{j}e^{-\mu t}-\mu e^{-\gamma_{j}t}),&\gamma_{j}/\mu\neq 1\\ 1-(1+\mu t)e^{-\mu t},\ \gamma_{j}/\mu=1,\end{array}\right.$$

$$⟨n⟩(t)={\int }_0^t\mathrm{d}sk(s)e^{W_{\mu }(s)-W_{\mu }(t)},$$

$$⟨m_j⟩(t)$$

$${\mathcal{K}}_j(s,s^{\prime },t)=k(s)q_j(s^{\prime })e^{W_{\mu }(s)-W_{\mu }(s^{\prime })+W_{{\gamma }_j}(s^{\prime })-W_{{\gamma }_j}(t)},$$

$$\frac{d⟨n^2⟩}{dt}$$

$$⟨n^2⟩(t)={\int }_0^t\mathrm{d}s[k(s)+(2k(s)+\mu (s))⟨n⟩(s)]e^{2W_{\mu }(s)-2W_{\mu }(t)}.$$

$$\frac{d⟨n{m}_j⟩}{dt}$$

$$\frac{d⟨m_1{m}_2⟩}{dt}$$

$$G^{(0)}(x,z_1,z_2,t)=\sum _{n,m_1,m_2}{x}^n{z}_1^{m_1}{z}_2^{m_2}{P}_{n,m_1,m_2}(t),$$

$$\frac{d{G}^{(0)}}{dt}$$

$$P_{m_1,m_2}(t)=\sum _{n=0}^{\infty }{P}_{n,m_1,m_2}(t),$$

$$\displaystyle G^{(0)}(z_{1},z_{2},t)=\lim_{N\rightarrow\infty}{\Big{[}}G^{(1)}_{N}(z_{1},z_{2},t){\Big{]}}^{N}.$$

$$G_N^{(1)}(z_1,z_2,t)=\sum _{\Theta }{\Psi }_N(\Theta )G_{\Theta }^{(2)}(z_1,z_2,t).$$

$$G_{\Theta }^{(2)}(z_1,z_2,t)=\prod _{j=1,2}{G}_{j|\Theta }^{(2)}(z_j,t).$$

$$G_{j|\Theta }^{(2)}=\underset{M\to \infty }{\lim }{\left[G_{M;j|\Theta }^{(3)}\right]}^M.$$

$$\displaystyle G^{(0)}(z_{1},z_{2},t)=\lim_{N\rightarrow\infty}{\Bigg{[}}\underbrace{\sum_{\Theta}\Psi_{N}(\Theta)\underbrace{\prod_{j=1,2}\underbrace{\lim_{M\rightarrow\infty}{\Big{[}}G^{(3)}_{M;j|\Theta}(z_{j},t){\Big{]}}^{M}}_{G^{(2)}_{j|\Theta}(z_{j},t)}}_{G^{(2)}_{\Theta}(z_{1},z_{2},t)}}_{G^{(1)}_{N}(z_{1},z_{2},t)}{\Bigg{]}}^{N}.$$

$$⟨m_j⟩$$

$$⟨m_j{⟩}_N^{(1)}$$

$$⟨m_j{⟩}_{\Theta }^{(2)}$$

$$⟨m_j{⟩}_{M|\Theta }^{(3)}=\frac{1}{M}{\int }_0^t\mathrm{d}\lambda \Theta (\lambda )q_j(\lambda )e^{W_{{\gamma }_j}(\lambda )-W_{{\gamma }_j}(t)}.$$

$$C_{1,2}=⟨m_1{m}_2⟩-⟨m_1⟩⟨m_2⟩=\underset{N\to \infty }{\lim }N⟨m_1{m}_2{⟩}_N^{(1)}.$$

$$⟨m_1{m}_2{⟩}_N^{(1)}=\sum _{\Theta }{\Psi }_N(\Theta )⟨m_1{⟩}_{\Theta }^{(2)}⟨m_2{⟩}_{\Theta }^{(2)}.$$

$$G_N^{(1)}(z_1,z_2,t)=\sum _{\Theta ,a,b}{z}_1^a{z}_2^b{\phi }_{N;(\Theta ,a,b)}(t).$$

$$G_N^{(1)}(z_1,z_2,t)=\sum _{\Theta }{\Psi }_N(\Theta )G_{\Theta }^{(2)}(z_1,z_2,t),$$

$$G_{\Theta }^{(2)}(z_1,z_2,t)=\sum _{a,b}{z}_1^a{z}_2^b{\psi }_{a,b|\Theta }(t).$$

$$\frac{d{G}_{\Theta }^{(2)}}{dt}$$

$$G_{\Theta }^{(2)}(z_1,z_2,t)=\prod _{j=1,2}{G}_{j|\Theta }^{(2)}(z_j,t),$$

$$\frac{d{G}_{j|\Theta }^{(2)}}{dt}=(z_j-1)(q_j(t)\Theta (t)-{\gamma }_j(t){\partial }_{z_j})G_{j|\Theta }^{(2)}.$$

Full text

1 Introduction

Gene expression is the biological process by which information from a gene is used to synthesize RNA macromolecules and proteins. With a few exceptions, until the 1990s, this process was commonly understood from “a deterministic viewpoint” [1, 2]. Since then, the combination of experimental and theoretical approaches has clarified that gene expression is often stochastic in nature (see [3, 4, 5, 6, 7] for review articles). The effect of fluctuations (noise) is usually limited when we are dealing with large numbers of molecules [8]. In cells however, wherein genes and mRNAs are often present in low numbers, stochasticity has an important role on cellular functions. The importance of fluctuations [9, 10, 11, 12, 13, 14, 15] can be illustrated by the observation that, amongst a genetically identical population in a homogenous environment, cell-to-cell variations in gene expression can result in phenotypic heterogeneity.

There exists various mechanisms, some more complex than others, allowing cells to tame and exploit randomness [16]. In order to unveil those processes, research efforts are directed on both experimental [17, 18, 19, 20] and theoretical fronts [21, 22, 23, 24, 25, 26, 27, 28, 29]. Collaborations between biologists, physicists and mathematicians aim to reveal the conditions under which transcriptional noise may, or may not, cascade to affect downstream genetic products.

The two stage and three stage models [30, 31, 32, 33, 34] give a minimalist description of the simplest yet non-trivial biological processes leading to gene expression. The two-stage model includes only transcription and translation processes, while the three-stage model also incorporates free and repressed states of the DNA promoter region. Analytical techniques and results [35, 36, 37] for the previously mentioned processes are the cornerstone for further theoretical developments. These models are the elementary bricks allowing for the construction of more complex reaction networks including non-exponential waiting times, transcriptional burst, feedback loops et cetera [38, 39, 40, 41, 42]. Even when it is possible to derive the exact mean and variance of protein and mRNA numbers, obtaining an exact closed-form expression for the generating function is often a challenging problem. The two-stage model is a perfect example. It has been the subject of numerous studies since the paper of Thattai and Oudenaarden [40] in 2001. The model presents linear propensities so that all moments can be derived exactly. Such problems (like the one considered in this paper) are said to be “exactly solvable”. Nevertheless, the exact generating function for the two stage model [37] was obtained only after ten years of extensive theoretical and experimental studies.

The search for exact solutions is often challenging because a small variation of a model’s definition can make analytical results unattainable. Typically methods aiming for a full characterisation of a given process (beyond results for the mean and variance), focus on the master equation approach and its partner equation for the generating function. Once the generating function is obtained, all moments are in principle known: given by successive derivatives. This approach can provide insights into different limiting cases and into the behaviour of the distribution in different regions of parameter space. It is important to mention that analytical results (for the probability distribution) have been obtained for a class of models such as monomolecular reactions systems [61] or deficiency zero networks [62, 63]. Unfortunately, once outside these classes there exists no systematic analytical recipe applicable independently of a model’s structure. Research efforts are naturally turning towards numerical simulations which though powerful “bring no intuition into the underlying […] interactions” [40]. To reach a better understanding one needs to investigate the joint distribution of mRNA and proteins, as well as temporal data, beyond the two-time autocorrelation function. As research progresses, emphasis is given to real time measurements with the hope to “expose the true cell dynamics buried in the average” [23]. Nowadays experimental advances allow for the count of individual molecules over time [43, 44, 45, 46] highlighting the need for both time-dependent and steady-state theoretical results .

In recent work, the partitioning of Poisson arrivals [47] was invoked to map Poisson processes to simple biological switches. This method is based on the separation of creation events (mRNA creation or protein creation) into independent groups. When applicable, this procedure leads to a mapping between creation/degradation process and a simple two-states biological switch. Applied to the two-stage model, this method led to the time dependent protein distributions [47] using already known results [21, 18, 36] for mRNA distributions in models with promoter-based regulation. The PPA-mapping needs, however, to be applied with care. It is important to warn the reader that, in a given model, not all creation/degradation process can be mapped onto a biological switch. For the PPA-mapping to apply, one needs to be able to partition a given creation event into independent processes. And for a given model, this will depend on upstream regulation of the creation/degradation process under consideration. This restriction is strong and appears as a serious limitation of the mapping applicability. We therefore need a more systematic way to use the PPA-mapping. As it is, the PPA mapping, can only be applied on models presenting a mixture of zero and first order reactions. It is unclear as if and how this method can be used or adapted to study models in the presence of feedback. In the simplest model describing bursty mRNA production, a variation of PPA mapping leads to an alternative derivation of the mRNA generating function [48]. But so far, this method has not been the subject of much attention and few models have been solved using this approach.

In direct connection with the applicability question of the PPA mapping is the inverse problem: Assuming the arbitrary partition of a creation event into two ‘types’ (type $1$ and type $2$), the latter process being itself regulated by an upstream mechanism, what correlation (between $1$ and $2$) should we expect? Is the correlation bounded? Does it vanish under particular conditions? Also, it is of experimental interest, to search for ways to infer the protein levels of a given type using measurement data on the other.

The model and methodology proposed in this paper were designed to (1) study correlations induced by the arbitrary partition of proteins arrival and (2) obtain the generating function making full use of the idea implicit in the PPA-mapping. The process we consider appears to be a simplified version of mechanisms involved in alternating splicing processes allowing a single gene to code for multiple proteins. With alternating splicing, a particular pre-mRNA can lead to different messenger RNAs, each being responsible for the production of isoform proteins (differing in their amino acid sequence). Results recently published in [49] focused on both bursty and constitutive pre-mRNA creation. One should mention that alternating splicing is far from being rare. Many genes have multiple splicing patterns [50, 51, 52] and numerous examples confirm that alternating splicing contributes to the development of cancer (see [53] and [54] for review articles). Our goal is to obtain the time dependent solution of the proposed model in term of the generating function. Our method is based on the construction of different mappings. Each transformation aims to reduce the study of a given model to the analysis of a simpler one. After a succession of transformations the problem is condensed to the study of two-state biological switches. The nesting between models and reduced models is reflected in a set of relations between generating functions. This hierarchy allows us to derive relations between mean numbers and higher order moments. We show that the PPA mapping allows us to consider arbitrary time dependent transition rates. Other studies such as [58, 59] and [60] have considered explicit time dependent parameters to investigate the effect of upstream hidden dynamics on downstream populations. Here, we obtain the time dependent generating function without solving any complicated differential equation but rather a simple first order equation (with time dependent coefficients). Accessing analytical results for time dependent coefficient provides a way to tackle models with noisy transition rates induced by hidden stochastic variables. In this paper, we show how results for the mean and correlations for time dependent model can be used to access the solution for the three-stage model.

The paper is organised as follows. In section (2.1), we start with the presentation of the model under consideration. We give the master equation governing the evolution of the probability distribution and present the solution of the first order moments. In section (2.2), we define the generating function, and outline the three different steps defining our method. Step $1$ and $3$ both describe the transformation under the PPA-mapping at different stages of the derivation. The intermediate step $2$ defines the decomposition of a given process over all possible histories. Each step is detailed in sections (2.4), (2.5) and (2.6). The succession of transformations takes us relatively far away from the solution of the original problem. To proceed further, we give in section (2.7), the probability associated to each relevant histories. Finally, the time dependent generating function is presented in section (3), where our result is generalised to arbitrary partition numbers. We show how for some particular cases our results match the known solution for the two-stage model (see reference [37] and [47]). Finally we discuss how to use results for time dependent parameters to extend this work to the three-stage model and others processes including additional random variables.

2 Theory

2.1 The model

The model under consideration describes the stochastic evolution of protein numbers in a variation of the two-stage model, for which proteins are arbitrarily separated into two groups (${\cal P}_{1}$ and ${\cal P}_{2}$). Note that our results will be easily generalized to an arbitrary number of protein types. We denote by ${\cal A}$ the upstream molecule regulating proteins production (see Figure (1)). Time dependent transition rates for proteins production and degradation are respectively written $q_{j}(t)$ and $\gamma_{j}(t)$ (with $j=1,2$). The level of regulator ${\cal A}$ is itself governed by the transition rates $k(t)$ and $\mu(t)$, respectively associated to creation and degradation (see table (1)). Because the method presented in this paper invokes mappings to other processes, it is convenient to refer to the original model as model-0.

The state of the system is, at any time, characterised by the numbers $n$, $m_{1}$ and $m_{2}$, of molecules ${\cal A}$ and proteins ${\cal P}_{1}$ and ${\cal P}_{2}$ respectively. We write $P_{n,m_{1},m_{2}}(t)$ the probability distribution of state $(n,m_{1},m_{2})$. We should keep in mind that the latter quantity is conditional on the initial state. In particular, we will consider $P_{0,0,0}(t=0)=1$. Further down in this paper we will explain how the latter initial condition is imposed by the initial state of the reduced model resulting from successive mappings. In order to consider other initial states, an extension of the proposed method is needed. This procedure would add an extra layer of complexity to the work presented here and is not discussed further in this paper. The probability distribution is governed by the master equation

[TABLE]

At this point, it is premature to solve the full master equation, and we start by deriving equations for mean population numbers. Let us denote $\langle{\cal O}\rangle$ the average of the observable ${\cal O}$, given by $\sum_{n,m_{1},m_{2}}{\cal O}(n,m_{1},m_{2})P_{n,m_{1},m_{2}}$. For mean numbers we derive the equations:

[TABLE]

If we restrict ourself to the well known case of constant reaction rates, with $\langle n\rangle(t=0)=\langle m_{j}\rangle(t=0)=0$, the solutions of these equations are:

[TABLE]

and

[TABLE]

with the stationary values $\langle n\rangle^{*}={k}/{\mu}$ and $\langle m_{j}\rangle^{*}={kq_{j}}/{\mu\gamma_{j}}$. At this stage, it is not hard to extend these results to time-dependent coefficients. With a bit of work, one can show that the solution of equation $\eqref{eq_2}$ is

[TABLE]

where we define $W_{\mu}(t)=\int_{0}^{t}\mathrm{d}\lambda\ \mu(\lambda)$. For arbitrary functions $k(t)$ and $\mu(t)$ it is however impossible to comment on the existence of a stationary state unless we assume the existence of the limits $k(t\rightarrow\infty)=k^{*}$ and $\mu(t\rightarrow\infty)=\mu^{*}$. It is here important to note that the solution of equation $\eqref{eq_3}$ can be expressed as a double integral. To proceed we use equation \eqref{solution_n_time_dep} with the substitution $k(t)\rightarrow q_{j}(t)\langle n\rangle(t)$ and $\mu(t)\rightarrow\gamma_{j}(t)$. A few lines of calculation leads to

[TABLE]

with kernel

[TABLE]

where $W_{\gamma_{j}}(t)=\int_{0}^{t}\mathrm{d}\lambda\ \gamma_{j}(\lambda)$. Interestingly, it is this quantity ${\cal K}_{j}(s,s^{\prime},t)$ which will reappear explicitely in the final expression for the generating function. Along the same lines, it is possible to push further, writing equations for second order moments such as $\langle n^{2}\rangle$

[TABLE]

Once again, the solution for constant coefficients is easy to derive and can be expressed as a function of $\langle n\rangle(s)$:

[TABLE]

To evaluate correlations of the form $\langle m_{1}m_{2}\rangle(t)$ and $\langle nm_{1}\rangle(t)$ we write

[TABLE]

Together with $d\langle m_{j}\rangle/dt$ ($j=1,2$), $d\langle n\rangle/dt$ and $d\langle n^{2}\rangle/dt$, Eq. (11) and (12) define a system of seven equations. Importantly, equations governing the evolution of correlators do not involve higher order moments. As a consequence, correlations at any order can be obtained by solving a finite set of equations [55, 56, 57]. Even if one considers constant reaction rates, the generalisation of the solution to three or more protein types is not trivial. It requires the solution of a new and bigger set of equations. One possible avenue is to pursue with approximation of the “mean field” type, which consists in assuming $\langle m_{1}m_{2}\rangle\simeq\langle m_{1}\rangle\langle m_{2}\rangle$. A priori, the later approximation holds for weakly correlated systems only. Hence we need to quantify correlation numbers in order to select the appropriate analytical methods.

2.2 The generating function

Let us start by defining the generating function of the original model

[TABLE]

which obeys the differential equation

[TABLE]

Focusing our attention on the numbers of proteins only, we define the marginal probability

[TABLE]

for which the generating function is $G^{(0)}(z_{1},z_{2},t)=G^{(0)}(1,z_{1},z_{2},t)$. In order to attain an analytical expression we will successively reduce the original model into simpler ones. To avoid confusion we choose to denote as model-1, model-2 and model-3, the processes which will be emerging from these successive mappings. We write $G^{(1)}$, $G^{(2)}$ and $G^{(3)}$ the generating functions for each model respectively. The following gives a short description of the steps taken in this paper, while each of them is further developed in sections (2.4), (2.5) and (2.6).

1. Step 1:

The PPA-mapping [47] is based on the partitioning of Poisson processes (see Figure (2)). It allows for simplification of the original problem to $N$ independent processes all identical to model-1 (see Figure (3) and section (2.4)). In the reduced model, the production of protein is regulated by a biological switch taking values $\theta=0$ (OFF) and $\theta=1$ (ON). Since $N$ appears as a parameter of the reduced model, we write $G^{(1)}_{N}(z_{1},z_{2},t)$ the generating function of model-1. The latter is related to the original generating function via:

[TABLE] 2. Step 2:

Denoting by $\Theta$ a particular history (or path) generated by the time evolution of the variable $\theta$, we define $\Psi_{N}(\Theta)$ to be the probability of a given path. Model-2 is defined for one particular history as if frozen (Figure (4)). We write $G^{(2)}_{\Theta}$ as the associated generating function and express $G^{(1)}_{N}$ as an average over all possible histories (see section (2.5))

[TABLE]

Once the differential equation for $G^{(2)}_{\Theta}$ has been derived, we will be able to show that protein numbers are uncorrelated in model-2. It follows that $G^{(2)}_{\Theta}$ can be expressed as the product of two functions, each associated to a given protein type:

[TABLE] 3. Step 3:

To access the solution of model-2, we exploit the PPA-mapping one more time. Splitting the creation process into $M$ independent processes, it ultimately reduces to the study of biological switches (see model-3 in Figure (5)). Writing $G^{(3)}_{M;j|\Theta}(z_{j},t)$ as the generating function of the switch $j$ ($j=1,2$), we show the relation (see section (2.6))

[TABLE]

Finally, nesting all steps together, the original generating function is given by

[TABLE]

2.3 Consequences: hierarchy in mean and correlation numbers

Before entering the heart of the subject with the application of the PPA-mapping, one can investigate consequences of these successive transformations. The nesting of generating functions allows us to derive direct relations between mean numbers in the different models. We write $\langle m_{j}\rangle_{N}^{(1)}$, $\langle m_{j}\rangle_{\Theta}^{(2)}$ and $\langle m_{j}\rangle_{M|\Theta}^{(3)}$ the mean numbers of $j$-proteins in model-1, 2 and 3 respectively. For simplicity, we choose to omit the superscript [math] so that $\langle m_{j}\rangle$ denotes the average number of proteins in the original model. To ease the notations further we choose not to make the time dependance explicit, since the relations derived bellow are true for all time $t$. Equations (16), (17) and (19) bring us to

[TABLE]

The calculation of $\langle m_{j}\rangle^{(3)}$ is a pretty simple affair. Each protein being reduced to a biological switch, $m_{j}^{(3)}$ is restricted to the value [math] and $1$. We give here, the expression of $\langle m_{j}\rangle^{(3)}$, for which the derivation is presented in section (2.6):

[TABLE]

To continue further, eq. (22) requires knowledge of the probability $\Psi_{N}(\Theta)$ for a given path. This is not particularly difficult as one only needs to consider paths with probability up to the order $1/N$ (see section (2.7)). Without further knowledge of the generating function, once $\Psi_{N}(\Theta)$ and $\langle m_{j}\rangle^{(3)}$ given, the reader can derive the time evolution for mean number of proteins using equations (21), (22) and (23). Practically, those steps give a convoluted way to reach the result already presented in (5). It however reflects on the strategy adopted here to access the generating function.

Considering the correlation function, with the help of Eq. (16), we can show

[TABLE]

The protein number being uncorrelated in model-2 we have $\langle m_{1}m_{2}\rangle^{(2)}_{\Theta}=\langle m_{1}\rangle_{\Theta}^{(2)}\langle m_{2}\rangle_{\Theta}^{(2)}$, which leads us to

[TABLE]

From the latter two equations, we conclude that $C_{1,2}>0$ unless at least one of $\langle m_{j}\rangle=0$ ($j=1,2$). Hence, there is no non-trivial point in parameter space such that the correlation between protein number vanishes. As a consequence, there is no region of the parameter space in which the mean field approach is valid. In [49], the authors focus on alternative splicing mechanism, investigating the stationary state of a slightly different model from the one presented here. This study considers the transition from a pre-mRNA to two different mature mRNAs. For constitutive expression (no bursty pre-mRNA creation), they show that (in the stationary state) the mRNA numbers (of type $1$ and $2$) are independent. They however observe, for bursty pre-mRNA production, the emergence of correlations between the two mature mRNA types.

2.4 First transformation: from model-0 to model-1

The PPA-mapping is based on the partitioning property of Poisson processes. Without entering into technical details, the mapping can be understood as follow:

Consider the creation/degradation process of regulator ${\cal A}$ (with rates $k(t)$ and $\mu(t)$). 2. 2.

Partition every creation events into $N$ ’types’ (Figure (2)). The partition is homogeneous so that each ${\cal A}$ molecule is equally likely to be assigned to a given type. It follows that the creation rate associated to a particular type is given by $k(t)/N$. 3. 3.

Take the limit $N\gg 1$. As a consequence, the probability of observing more than one ${\cal A}$ molecule of a particular type can be neglected. It follows that the random variable describing the number of molecules ${\cal A}$ (of a given type) is restricted to the value [math] or $1$.

Model-1 as defined under this procedure is illustrated on Figure (3). Note that $N$ appears as a parameter in the reduced model. Along the lines presented in [47] we write $G^{(0)}=[G^{(1)}_{N}]^{N}$. Equation (2.2) shows that $G^{(1)}_{N}$ obeys the same differential equation under the transformation $k(t)\rightarrow k(t)/N$. As a consequence the probability of observing (in the reduced model) more than one ${\cal A}$ molecule is of order $1/N^{2}$ and can be neglected as $N\rightarrow\infty$. While the previous logical argument shows how model-1 is emerging from model-0, an alternative derivation, based on the probability distribution instead of the generating function, allows for the reversed construction: building model-0 starting with $N$ independent model-1. This derivation, not presented in the literature so far, is presented in an appendix.

2.5 Second transformation: from model-1 to model-2

Let us remind the reader that $\theta$ is the new stochastic variable (taking value in $\{0,1\}$) emerging in model-1. The decomposition over all possible histories, generated by the variable $\theta$, emerges from the use of conditional probabilities. To be more explicit we write $\Theta$ as a particular history associated to the variable $\theta$. For a given path, we write $\Theta(t)$ as the value taken by the random variable $\theta$ at time $t$. We continue further by writing $\varphi_{N;(\Theta,a,b)}(t)$ as the probability associated to a particular history $\Theta$ and protein numbers $a$ and $b$. The generating function $G^{(1)}_{N}$ can be rewritten as

[TABLE]

Defining $\psi_{a,b|\Theta}(t)$ as the conditional probability on $\Theta$ while $\Psi_{N}(\Theta)$ is the probability of a given history, the equality $\varphi_{N;(\Theta,a,b)}(t)=\Psi_{N}(\Theta)\psi_{a,b|\Theta}(t)$ leads to

[TABLE]

with

[TABLE]

For a known history $\Theta(t)$, we have

[TABLE]

Note that in the last equation $G^{(2)}_{\Theta}$ is clearly independent of $N$ and so is the conditional probability $\psi_{a,b|\Theta}$. The dependence in $N$ is now carried by the probability $\Psi_{N}$. At this point we see that $G^{(2)}_{\Theta}$ can be written as

[TABLE]

where each generating function is governed by

[TABLE]

Thereupon the two protein numbers are uncorrelated in model-2.

2.6 Third transformation: from model-2 to model-3

To reach the expression of $G^{(2)}_{\Theta}$, one applies the PPA-mapping one more time. This will reduce the original model to the study of biological switches (Figure (5)). For each protein type ${\cal P}_{j}$ ($j=1,2$), we once again, choose to partition every creation event into $M$ groups. The partition being homogeneous, each protein is equally likely to be assigned to a given group. The creation rate for a particular group is given by $q_{j}(t)/M$. Taking the limit $M\gg 1$ allows us to neglect the creation of more than one protein in each group. To be explicit one writes $G^{(2)}_{j|\Theta}=[G^{(3)}_{M;j|\Theta}]^{M}$ in equation (32). This procedure leads to the same differential equation with the transformation $q_{j}(t)\Theta(t)\rightarrow q_{j}(t)\Theta(t)/M$. Hence, in the limit $M\rightarrow\infty$, the number of proteins of type $j$ are restricted to [math] and $1$. It follows that

[TABLE]

The function $G^{(3)}_{M;j|\Theta}$ describes the dynamics of a two-state model and can be written has

[TABLE]

where $f_{M;j|\Theta}$ is the probability to find the switch $j$ in the ON-state, knowing the history $\Theta$. The latter is the solution of the following equation

[TABLE]

We now have reached the point where one needs to define the initial state. We choose to consider $f_{M;j|\Theta}(t=0)=0$. Let us remind the reader that, in model-3, the total number of switches $j$ in the ON-state equals the number of proteins ${\cal P}_{j}$ in model-2. The hierarchy builds up to the number of proteins in model-0. As we look at equation (21), (22) and (23), we see that choosing (at time $t=0$) all switches ($j=1,2$) in the OFF-state, imposes the following initial state on to the original model

[TABLE]

A simple calculation gives

[TABLE]

which, to the first order in $1/M$, simplifies to

[TABLE]

with

[TABLE]

Nesting equation (38) into (34) leads to:

[TABLE]

With Eq. (33) the latter result allows us to write

[TABLE]

2.7 Summing over all histories

In order to derive the generating function $G^{(1)}_{N}$ from $G^{(2)}_{\Theta}$, using equation (17), we focus on the expression of the probability $\Psi_{N}(\Theta)$ for all relevant histories $\Theta$. As mentioned earlier, we simply need to evaluate $\Psi_{N}(\Theta)$ up to the order $1/N$. We choose to consider the initial state $\Theta(t=0)=0$. Amongst $N$ identical models, the number of switches in the ON-states defines the number $n$ of molecules ${\cal A}$. It follows that the initial state must satisfy:

[TABLE]

Together, equations (36) and (42) fully specify the initial state, so that $P_{0,0,0}(t=0)=1$. We remind the reader that the probability of the transition form $ON\rightarrow OFF$ (between time $t$ and $t+\delta t$) is given by $\mu(t)\delta t$. In addition, the probability of observing the transition $OFF\rightarrow ON$ is given by $k(t)\delta t/N$. At the first order in $1/N$, three different types of histories are relevant. They are symbolically represented by
,
and
and detailed in table (2). The probability associated to each path is

[TABLE]

It is particularly useful to rewrite the last equation as

[TABLE]

which can be used to verify the conservation of probability:

[TABLE]

The latter relation confirms that all relevant paths have been taken into consideration. The expression of $G^{(1)}_{N}$ is symbolically given by

[TABLE]

The explicit calculation (using $G^{(2)}_{{\makebox[8.0pt][l]{ \leavevmode\hbox to5.52pt{\vbox to0.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}{{}{}}{{}}{} {}{}{}{}{}{}{{}}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.66365pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{4.66365pt}{0.0pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} }}}=1$) points us to

[TABLE]

with

[TABLE]

3 Result: final expression of $G^{(0)}(z_{1},z_{2},t)$

The methodology presented in the previous section leads us to the following generating function (obtained by taking the limit $N\rightarrow\infty$ in equation (49)):

[TABLE]

where

[TABLE]

with the same kernel ${\cal K}_{j}$ defined in eq. $\eqref{eq_kernel}$).Finally, $G^{(2)}_{\Theta}$ for the path
is explicitly given by

[TABLE]

To keep notation as compact as possible, we will simply write ${\cal K}_{j}$, omitting the variables $s,s^{\prime}$ and $t$. The relation between $G^{(0)}$ and $G^{(2)}_{\Theta}$ is particularly interesting and allows the kernel ${\cal K}_{j}$ to play a key role in a new set of relations between model-2 and the original model. In model-0, the mean $\langle m_{j}\rangle$ and correlation $\langle m_{j}(m_{j}-1)\rangle$ are obtained using $\langle m_{j}\rangle(t)=\partial_{z_{j}}G^{(0)}(z_{1},z_{2},t)|_{z_{1},z_{2}\rightarrow 1}$ and $\langle m_{j}(m_{j}-1)\rangle(t)=\partial_{z_{j}}^{2}G^{(0)}(z_{1},z_{2},t)|_{z_{1},z_{2}\rightarrow 1}$. For the mean, a simple calculation leads to $\langle m_{j}\rangle(t)=\int_{0}^{t}\mathrm{d}s\ \int_{s}^{t}\mathrm{d}s^{\prime}\ {\cal K}_{j}(s,s^{\prime},t)$ (identical to eq. (7)). For the variance, defined by ${\rm{Var}}[m_{j}](t)=\langle m_{j}^{2}\rangle(t)-[\langle m_{j}\rangle(t)]^{2}$, writing $\langle m_{j}\rangle^{(2)}_{{\makebox[15.0pt][l]{ \leavevmode\hbox to8.09pt{\vbox to2.96pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}{}{{}}{} {}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{2.56181pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{2.56181pt}{0.0pt}\pgfsys@lineto{2.56181pt}{2.56181pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{2.56181pt}{2.56181pt}\pgfsys@lineto{5.12364pt}{2.56181pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{5.12364pt}{2.56181pt}\pgfsys@lineto{5.12364pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}{}{}{}{{}}\pgfsys@moveto{5.12364pt}{0.0pt}\pgfsys@lineto{7.22546pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{7.22546pt}{0.0pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} }}}=\partial_{z_{j}}G^{(2)}_{{\makebox[15.0pt][l]{ \leavevmode\hbox to8.09pt{\vbox to2.96pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}{}{{}}{} {}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{2.56181pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{2.56181pt}{0.0pt}\pgfsys@lineto{2.56181pt}{2.56181pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{2.56181pt}{2.56181pt}\pgfsys@lineto{5.12364pt}{2.56181pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}\pgfsys@moveto{5.12364pt}{2.56181pt}\pgfsys@lineto{5.12364pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } {}{{}}{} {}{}{}{}{}{}{{}}\pgfsys@moveto{5.12364pt}{0.0pt}\pgfsys@lineto{7.22546pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{7.22546pt}{0.0pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} }}}|_{z_{1},z_{2}\rightarrow 1}$, we obtain :

[TABLE]

while the correlation function $C_{i,j}(t)=\langle m_{i}m_{j}\rangle(t)-\langle m_{i}\rangle(t)\langle m_{j}\rangle(t)$ becomes

[TABLE]

With the mean and correlation numbers in hand, the variance can be reached easily using ${\rm{Var}}[m_{j}](t)=\langle m_{j}\rangle(t)+C_{j,j}(t)$. It is clear that the generating function can be generalized to an arbitrary number $J$ of proteins by replacing $\sum_{j=1,2}\rightarrow\sum_{j=1}^{J}$ in equation (51) and (53). In this situation, the generating function depends of $J$ variables: $z_{1}$, $z_{2}$, …, $z_{J}$.

3.1 Constant reaction rates

Let us first focus on simplifications occurring when all transition rates are constant. We have $W_{\mu}(t)=\mu t$, $W_{\gamma_{j}}(t)=\gamma_{j}t$ and ${\cal K}_{j}(s,s^{\prime},t)=kq_{j}e^{\mu(s-s^{\prime})+\gamma_{j}(s^{\prime}-t)}$. It is then convenient to define $u=e^{\gamma_{j}(s-t)}$ and $v=e^{\gamma_{j}(s^{\prime}-t)}$ so that $A_{j}$ takes the form

[TABLE]

with

[TABLE]

At this stage it is not hard to show that the mean protein number is given by (5). Figure (6) confirms the validity of our results. For correlation numbers, it is useful to define the normalised function $\tilde{C}_{i,j}(t)=C_{i,j}(t)/{\langle m_{i}\rangle^{*}\langle m_{j}\rangle^{*}}$. For $\gamma_{i}+\gamma_{j}\neq\mu$, $\gamma_{i}\neq\mu$ and $\gamma_{j}\neq\mu$ the latter quantity is given by

[TABLE]

The case $\gamma_{i}=\mu$ (or $\gamma_{j}=\mu$) has to be treated separately. To proceed one can (1) set $\gamma_{i}=\mu$ in the kernel $K_{i}$, or alternatively (2) write $\gamma_{i}=\mu+\epsilon$ and take the limit $\epsilon\rightarrow 0$. The case $\gamma_{i}+\gamma_{j}=\mu$ (or $2\gamma_{i}=\mu$ when considering the variance) has to be treated similarly. Those limits lead to relatively more compact expressions, for example when $\gamma_{i}=\gamma_{j}=\mu$, we have:

[TABLE]

The agreement (for all time $t$) between analytical expressions and numerical simulations can be seen in Figure (7). In the limit $t\rightarrow\infty$ all expressions of $\tilde{C}_{i,j}$ converge to a single form. In the stationary state, the correlation function has a unique expression

[TABLE]

with $\eta_{i}=\gamma_{i}/\mu$ and $\eta_{j}=\gamma_{j}/\mu$. So that

[TABLE]

For homogeneous degradation rates ($\gamma_{i}=\gamma_{j}=\gamma$), the correlation is invariant under the exchange $\gamma\leftrightarrow\mu$. The last equation clearly shows that the correlation function does not vanish (unless one out of $\langle m_{i}\rangle^{*}$ and $\langle m_{j}\rangle^{*}$ vanishes). We note that the correlation $C^{*}_{i,j}$ is strictly monotonic (decreasing) in terms of $\gamma_{i}$ and $\gamma_{j}$ (keeping all other parameters constant). As a consequence, if one can estimate lower and upper bounds of both $\gamma_{i}$ and $\gamma_{j}$ it is, in principle, possible to restrain the range of correlation values to an interval: $[C^{*}_{min},C^{*}_{max}]$. In addition, we observe, for a fixed value of $\eta_{j}$, that the correlation $\tilde{C}^{*}$ presents a maximum at $(\eta_{i})_{max}=\eta_{j}+\sqrt{2\eta_{j}(\eta_{j}+1)}$. If $C^{*}$ is strictly monotonic, it is when varying $\gamma_{i}$ while keeping $\langle m_{i}\rangle^{*}$ constant that a non monotonic behaviour is observed. In this case, $C^{*}_{i,j}$ can be rewritten as

[TABLE]

and presents a maximum in $(\eta_{i})_{max}$. Keeping both protein levels $\langle m_{i}\rangle^{*}$ and $\langle m_{j}\rangle^{*}$ constant, the correlation function becomes $C^{*}_{i,j}=\langle m_{i}\rangle^{*}\langle m_{j}\rangle^{*}\tilde{C}^{*}_{i,j}$. Looking for an upper bound into the $(\eta_{i},\eta_{j})$-plane, one needs to solve $\partial_{\eta_{i}}C^{*}=0$ and $\partial_{\eta_{j}}C^{*}=0$ simultaneously. However there are no strictly positive solutions to the latter system of equations. Hence, under this constrain, $C^{*}_{i,j}$ does not present a maximum when varying both $\gamma_{i}$ and $\gamma_{j}$.

3.1.1 The $2$-stage model: $J=1$

In the case $J=1$, the model with constant transition rates, reduces to the conventional two-stage model. To pursue, we define $r=\mu/\gamma$ and $\delta(z)=q(z-1)/\gamma$. We can show that our result leads to the solution first presented in [37] and later in [47]:

[TABLE]

where $\hphantom{.}{}_{1}F_{1}$ is the confluent hypergeometric function. As it is, the identity between equation (51) (for one protein type only) and equation (63) is not obvious. To proceed, we use the Taylor expansion of $e^{\delta v}$ and $e^{-\delta u}$ and write $\epsilon=e^{-\gamma t}$ which we assume small compared to one. Considering $r=\mu/\gamma\neq 1$ and keeping the lowest order in $\epsilon$ (see appendix) we show that

[TABLE]

with $G^{*}(z)$ given by equation (63) and

[TABLE]

such that $\lim_{z\rightarrow 1}H(z,t)=\lim_{t\rightarrow\infty}H(z,t)=1$. The latter approximation leads to:

[TABLE]

in agreement with equation (5). The case $\gamma=\mu$ is treated separately in appendix.

3.1.2 Homogeneous degradation rates: $\gamma_{j}=\gamma$, $\forall j$

When dealing with $J$ protein types ($J>1$) and homogeneous degradation rates ($\gamma_{j}=\gamma$ $\forall j$), the generating function reduces to a form close to the one previously obtained for the two-stage model. Defining $\Delta(\{z_{j}\})=\sum_{j}q_{j}(z_{j}-1)/\gamma$, we can show that $G^{(0)}(\{z_{j}\},t)$ is given by equation (64) under the substitution $\delta(\{z_{j}\})\rightarrow\Delta(\{z_{j}\})$. It follows that

[TABLE]

The generating function ${\cal G}$, associated to the total number of proteins ($M=\sum_{j}m_{j}$), is defined by ${\cal G}(z)=\sum_{M}P_{M}z^{M}$, with

[TABLE]

We see that ${\cal G}$ is given by ${\cal G}(z,t)=G^{(0)}(\{z_{j}=z\},t)$:

[TABLE]

with the average creation rate defined by $J\bar{q}=\sum_{j}q_{j}$. The mean of total protein number ($M=\sum_{j}m_{j}$) satisfies ${\langle M\rangle}/{J}={k}{\bar{q}}/({\mu}{\gamma})$.

3.2 Time dependent transition rates: a bridge towards other models

Results for time dependent parameters allow for the study of fluctuations (induced by hidden variables) in production and/or degradation rates. In a recent paper Dattani and Barahona [60] proposed a framework to model gene expression with stochastic or deterministic transcription and degradation rates. Along the same lines, let us start by defining random variables $x_{\phi}$, for all parameters $\phi$ of the model ($\phi=k,\mu,q_{1},q_{2},...,\gamma_{1},\gamma_{2},...$). We choose to write $\phi(t)=\phi_{0}+\phi_{1}x_{\phi}(t)$ with $\phi_{0},\phi_{1}\in{\mathbb{R}}$, but other functional forms could be chosen. More explicitly, we have for $k(t)$: $k(t)=k_{0}+k_{1}x_{k}(t)$. A particular time history of the random variable $x_{k}(t)$ is written $X_{k}=(x_{k}(t)|\forall t)$. Identically, we write $X_{\phi}=(x_{\phi}(t)|\forall t)$ the history for random variable $x_{\phi}$. Finally, we define ${\mathbb{X}}$ as the set of histories ${\mathbb{X}}=\{X_{k},X_{\mu},X_{q_{1}},...X_{\gamma_{1}},...\}$ so that the generating function is now explicitly dependent on ${\mathbb{X}}$: we write $G^{(0)}_{{\mathbb{X}}}$. Writing ${\cal P}({\mathbb{X}})$ as the probability of the set of histories, the new generating function is

[TABLE]

where $\sum_{{\mathbb{X}}}$ symbolically represents the sum over all possible histories of all parameters. The mean protein number is given by

[TABLE]

with $\langle m_{j}\rangle_{{\mathbb{X}}}=\partial_{z_{j}}G^{(0)}_{{\mathbb{X}}}|_{\forall z_{i}=1}$. To give a concrete example, we will restrict ourself to time dependent production rate $k(t)$ while all other transition rates are constant. We show how our results bridge towards the three-stage model, allowing us to access the exact mean and correlation functions. In this situation ${\cal K}_{j}(s,s^{\prime},t)=k(s)q_{j}e^{\mu(s-s^{\prime})+\gamma_{j}(s^{\prime}-t)}$. Choosing the appropriate function $k(t)$ can give information on the behaviour induced by state fluctuation of the DNA operational site (Figure (8)). We will write $\langle m_{j}\rangle_{X_{k}}$ the mean protein numbers, for a particular history ${X_{k}}$. For $\gamma_{j}\neq\mu$ and without restriction on $k(t)$, we show that

[TABLE]

with the convolution product

[TABLE]

For $\gamma_{j}=\mu$ we write $\gamma_{j}=\mu+\epsilon$ in (72) together with the limit $\epsilon\rightarrow 0$. 111Along the same line, when considering time dependent production rate $q_{j}(t)$ while keeping all other parameters constant, we get an equation similar to (72): $\langle m_{j}\rangle_{X_{q_{j}}}(t)=({k}/{\mu})\left[q_{j}-q_{j}e^{-}_{\mu}\right]\star e^{-}_{\gamma_{j}}.$

The previous equation becomes $\langle m_{j}\rangle_{X_{k}}(t)=q_{j}\left(-\frac{\partial}{\partial\mu}\right)(k\star e^{-}_{\mu})(t)$. For all values of $\gamma_{j}$ and $\mu$, the Laplace transform of the mean number of protein $L[\langle m_{j}\rangle]$ simplifies to a single expression:

[TABLE]

with $L[k]$ as the Laplace transform of $k(t)$. Assuming the limit $k(t\rightarrow\infty)=k^{*}$ exists, the final value theorem leads to $\langle m_{j}\rangle^{*}=k^{*}q_{j}/(\gamma_{j}\mu)$. With the three-stage model in mind, we set $k_{0}=0$ and write $k(t)=k_{1}x_{k}(t)$, where the random variable $x_{k}(t)$ takes value in $\{0,1\}$. It describes the possible states, active ($x_{k}=1$) or inactive ($x_{k}=0$), of the promoter region. Governed by a simple two state dynamics (with transition rates $W_{0\rightarrow 1}=\alpha$ and $W_{1\rightarrow 0}=\beta$) the variable $x_{k}(t)$ “oscillates” between those states (see illustration (8)). This motivates in [59] the choice of a sinusoidal function: $k(t)=c_{1}\sin(\omega t+\phi)+c_{2}$. However, the time evolution of the variable $x_{k}(t)$ is stochastic and, starting from initial condition $\langle x_{k}\rangle(t=0)=0$, it satisfies

[TABLE]

with stationnary state $\langle x_{k}\rangle^{*}=\alpha/(\alpha+\beta)$. It follows that the mean for the three stage model is given by $\langle\langle m_{j}\rangle\rangle$ representing the average over the history of the variable $x_{k}$. For $\gamma_{j}\neq\mu$ we have to evaluate

[TABLE]

with $\langle\langle m_{j}\rangle\rangle^{*}=\langle x_{k}\rangle^{*}k_{1}q_{j}/(\mu\gamma_{j})$. A simple calculation leads to the exact expression:

[TABLE]

as long as $\gamma_{j}\neq\mu$, $\mu\neq\alpha+\beta$ and $\gamma_{j}\neq\alpha+\beta$. Once again, the Laplace transform gives one single expression valid in all parameter space:

[TABLE]

This result is not new and could have alternatively been obtained by writing $d\langle\langle m_{j}\rangle\rangle/dt=q_{j}\langle\langle n\rangle\rangle-\gamma_{j}\langle\langle m_{j}\rangle\rangle$, which solution is $\langle\langle m_{j}\rangle\rangle=q\int_{0}^{t}\mathrm{d}s\ \langle\langle n\rangle\rangle(s)e^{-\gamma_{j}(t-s)}$ and using the time evolution of mRNA level $\langle\langle n\rangle\rangle$ (presented in [21]). Figure (9) shows agreement between analytical predictions and numerical simulations. If it is mathematically convenient to consider the time evolution starting from an “empty” initial state (all stochastic variable to zero), this situation does not seem to be biologically relevant. One could however, consider the similar scenario starting from the state $x_{k}=0$, with initial numbers $n,m_{1},m_{2}$ of $A$ macromolecules and proteins. Solving this new problem requires a different approach based on an variation of the PPA mapping, which is not considered in this paper.

Finally, let us discuss how to infer on correlation numbers between protein types. First, with the help of

[TABLE]

we can use equation (55) to express the correlation function $C_{i,j|{X_{k}}}=\langle m_{i}m_{j}\rangle_{X_{k}}(t)-\langle m_{i}\rangle_{X_{k}}(t)\langle m_{j}\rangle_{X_{k}}(t)$:

[TABLE]

One should note that the variance is given by $C_{j,j|{X_{k}}}(t)+\langle m_{j}\rangle_{X_{k}}(t)$ and can be evaluated using equations (72) and (80). For singular cases $\gamma_{i}=\mu$, $\gamma_{j}=\mu$ or $\gamma_{i}+\gamma_{j}=\mu$ a similar expression can be derived from the previous equation taking the limit appropriately. To continue further one has to proceed more carefully. In fact the correlations in the model presented in figure 6 are defined by ${\cal C}_{i,j}=\langle\langle m_{i}m_{j}\rangle\rangle-\langle\langle m_{i}\rangle\rangle\langle\langle m_{j}\rangle\rangle$, which can be expressed using $\langle C_{i,j}\rangle$ (the average of $C_{i,j|X_{k}}$ over the history $X_{k}$):

[TABLE]

One can notice that $\langle C_{i,j}\rangle$ (just like $\langle\langle m_{i}\rangle\rangle$ and $\langle\langle m_{j}\rangle\rangle$) is a functional of the mean $\langle x_{k}\rangle(t)$, and can be evaluated easily. The challenge comes from the term $\langle\langle m_{i}\rangle\langle m_{j}\rangle\rangle$ as it requires knowledge of correlators $\langle x_{k}(s)x_{k}(s^{\prime})\rangle$:

[TABLE]

In fact the $x_{k}$-correlation can be calculated exactly. For $s<s^{\prime}$, it is given by

[TABLE]

This last expression, together with the help of (80) and (77) allow for the evaluation of ${\cal C}_{i,j}$. To compare correlations ${\cal C}_{i,j}$ for the three-stage model with $C_{i,j}$ (61) for the two-stage model, we impose the equality $\alpha k_{1}=(\alpha+\beta)k$ which assures, in both models, identical regulator, and proteins levels. It follows that ${{\cal C}_{i,j}^{*}}/{C_{i,j}^{*}}=1+R$, where $R$ has a cumbersome expression, dependent on all parameters but $q_{i}$ and $q_{j}$. It vanishes for $\beta=0$ (as well as $k=0$) and satisfies $R>0$ for all other finite parameter values. Hence we conclude that, in the stationary state, for identical regulator and protein levels, correlations between protein numbers are higher in the model with promoter-based regulation: ${\cal C}_{i,j}^{*}>C_{i,j}^{*}$. Restraining ourself to homogeneous degradation rates ($\gamma_{i}=\gamma_{j}=\gamma$) the expression for $R$ is more manageable :

[TABLE]

Once again, we note that $R$ and ${\cal C}_{i,j}^{*}$ are invariant under the exchange $\mu\leftrightarrow\gamma$. Figure (9) shows a comparison between the time evolution of protein number in the two and three stage model. Our data validate the equality ${{\cal C}_{i,j}^{*}}/{C_{i,j}^{*}}=1+R$ and seem to indicate that ${\cal C}_{i,j}(t)/C_{i,j}(t)\simeq 1+R$ is an acceptable approximation, at least for the set of parameter selected.

4 Conclusion

In this paper, we present a variation on the two-stage model including the arbitrary partition of protein and arbitrary time dependent parameters. The mechanism considered is similar to the one involved in alternative splicing [49]. Our aim was to develop an analytical approach leading to the time dependent generating function, based on the PPA-mapping, and which did not require knowledge of already known results. Constructed on a succession of transformations, our work simplifies the original problem to the analysis of two-states biological switches. As a consequence, a series of different reduced models emerges, with clear relations linking their respective generating functions. In particular we show how the kernel ${\cal K}_{j}$ plays a important role in the final expression of the generating function. We show that the emerging hierarchy allows us to connect mean numbers and higher order moments between models. This leads to an explicit relation between correlation $C_{i,j}$ in the original model and mean protein number in model-2. For constant parameter values, we derived the exact time dependent expression of correlation $C_{i,j}$. We note that particular cases such as $\gamma_{j}=\mu$, $\gamma_{j}=\mu/2$ or $\gamma_{1}+\gamma_{2}=\mu$ have to be treated separately. However, in the stationary state this distinction vanishes as each function converges towards a common asymptotic expression. Considering constant parameters, for $J=1$ or for homogeneous degradation rates ($\gamma_{j}=\gamma$ $\forall j$), we show how our results reduce to the solution of two stage model. Finally, we show how results for arbitrary time dependent transition rates can be used to study models presenting parameter fluctuations induced by hidden random variable. We give an explicit example by extending our results for the mean and correlation to the three state model. In particular we show that promoter based regulation leads to higher protein-protein correlations ${\cal C}^{*}_{i,j}>C_{i,j}^{*}$. The method presented here is applicable in other scenarios (with zero and first order reactions) and may be use as a guide in the study of other biological systems. We hope this methodology will contribute to the development of new analytical avenues for future research.

Acknowledgements

The authors would like to thanks the AMRC in Coventry and the Stat. Phys. Group in Nancy for their support and constant efforts to protect curiosity driven research. In particular, T. Platini extends his acknowledgements to K. Aujogue, R. Low, S. Mukherjee, S. Vantieghem and R. Kulkarni.

5 Appendix

5.1 From model-0 to model-1: a detailed derivation

In the following we aim to construct the original problem (model-0) starting with $N$ independent models, all identical to model-1. The derivation is presented for constant parameter, but could easily be extended to time dependent parameter. We will start with the definition of $\phi_{N;(\theta,a,b)}(t)$: the probability of finding, in model-1, the biological switch in the state $\theta$ together with $a$ and $b$ proteins of type ${\cal P}_{1}$ and ${\cal P}_{2}$ respectively. We continue by defining $\Phi_{N}(n,m_{1},m_{2})$ as the probability of finding amongst $N$ independent identical models (each labeled with subscript $\nu$) the total numbers of active switches $n=\sum_{\nu}\theta_{\nu}$ and proteins $m_{1}=\sum_{\nu}a_{\nu}$, $m_{2}=\sum_{\nu}b_{\nu}$. Finally when, taking the limit $N\rightarrow\infty$, we will be able to show that $\lim_{N\rightarrow\infty}\Phi_{N}(t)$ obeys the master equation $\eqref{eq_Master_eq}$. In other words $P_{n,m_{1},m_{2}}=\lim_{N\rightarrow\infty}\Phi_{N;(n,m_{1},m_{2})}$.

To start let us try to keep notation as compact as possible by defining the flip operator ${\hat{F}}$ as the operator acting on the triplet $(\theta,a,b)$ and such that ${\hat{F}}(\theta,a,b)=(1-\theta,a,b)$. In addition, we write ${\hat{A}}^{\pm}$ and ${\hat{B}}^{\pm}$ as the operators defined by

[TABLE]

Let $\phi_{N;(\theta,a,b)}(t)$ be the probability distribution associated to model-1. It obeys the following equation

[TABLE]

Considering $N$ independent but identical models, we define ${\mathbb{S}}$ as the set of triplets: ${\mathbb{S}}=\{(\theta_{\nu},a_{\nu},b_{\nu});\nu=1,2,\ldots,N\}$. Note that ${\mathbb{S}}$ gives a full description of the state of all $N$ independent reduced models. In addition, we will write $\theta_{\nu}({\mathbb{S}})$, $a_{\nu}({\mathbb{S}})$ and $b_{\nu}({\mathbb{S}})$ the variable $\theta$, $a$ and $b$ in the $\nu^{th}$ triplet of the set ${\mathbb{S}}$. The definitions of the operator ${\hat{F}}$, ${\hat{A}}^{\pm}$ and ${\hat{B}}^{\pm}$ are extended onto the set ${\mathbb{S}}$. So that ${\hat{F}}_{\nu}$, ${\hat{A}}^{\pm}_{\nu}$ and ${\hat{B}}^{\pm}_{\nu}$ act on the $\nu^{th}$ triplet of the set ${\mathbb{S}}$, leaving all others unchanged. We can now define the overall probability $\phi_{N,{\mathbb{S}}}=\prod_{\nu=1}^{N}\phi_{N,(\theta_{\nu},a_{\nu},b_{\nu})}$, which evolution is governed by

[TABLE]

With a little bit of effort, the latter equation leads to

[TABLE]

Moving forward, we define the probability distribution $\Phi_{N}$ as

[TABLE]

where the constraint $\Delta$ is defined by

[TABLE]

where $\delta$ is the Kronecker symbol. To keep notations as compact as possible we write $\Delta^{(3)}=\Delta^{\{\theta_{\nu}\}}_{n}\Delta^{\{a_{\nu}\}}_{m_{1}}\Delta^{\{b_{\nu}\}}_{m_{2}}$. The master equation for the distribution $\Phi_{N}$ is

[TABLE]

Every sum, for which $\phi_{N;{\mathbb{S}}}$ appears explicitly can be easily evaluated. As an example we give here the second term of the first line in the previous equation. Using the constraint $\Delta^{\{\theta_{\nu}\}}$, which impose $\sum_{\nu}\theta_{\nu}=n$, we have

[TABLE]

When $\phi_{N,{\mathbb{S}}}$ does not appear explicitly we need to re-labeled the sum over ${\mathbb{S}}$. As an example we present details to the calculation of the first term of the first line in which we have $\phi_{N,{\hat{F}}_{i}{\mathbb{S}}}$. Defining $\tilde{\mathbb{S}}={\hat{F}}_{i}{\mathbb{S}}$ we have ${\hat{F}}_{i}\tilde{\mathbb{S}}={\mathbb{S}}$ so that

[TABLE]

In relabelling ${\mathbb{S}}$ to $\tilde{\mathbb{S}}$ the expression of $\Delta^{\{\theta_{\nu}\}}$ has changed. To keep track of this change we will replace $\Delta^{\{\theta_{\nu}\}}\rightarrow\tilde{\Delta}^{\{\theta_{\nu}\}}$ where

[TABLE]

It follows that

[TABLE]

Note that the only elements which will contribute are such that $\theta_{i}(\tilde{\mathbb{S}})=0$, which allows us to write

[TABLE]

Finally, we are able to express the first summation in term of $\Phi_{N,(n,m_{1},m_{2})}$:

[TABLE]

Proceeding along the same line for every summation symbol we have

[TABLE]

Taking the limit $N\rightarrow\infty$, we see that the latter equation converges towards the master equation (2.1). In other words $P_{n,m_{1},m_{2}}=\lim_{N\rightarrow\infty}\Phi_{N;(n,m_{1},m_{2})}$, from which it naturally follow $G^{(0)}=\lim_{N\rightarrow\infty}\left(G^{(1)}_{N}\right)^{N}$.

5.2 Large time approximation

The large time approximation, for one protein type only, is obtained by writting $r=\mu/\gamma$, $\delta=q(z-1)/\gamma$ and $\epsilon=e^{-\gamma t}$ so that

[TABLE]

with

[TABLE]

Using the Taylor expansion leads to

[TABLE]

for which there is no particular problem unless in the last integral we have $n-r=-1$ for some value of $n$. Avoiding this situation, by choosing $r\notin\mathbb{N}$, leads to

[TABLE]

with

[TABLE]

Keeping the lowest order in $\epsilon$ we get $(z-1)A\simeq I_{F}+h$ with

[TABLE]

It follows that $G^{(0)}(z,t)\simeq G^{*}(z)H(z,t)$, with $H=e^{h}$ as presented in equation (65). Going back to equation (102), we can work with $r=1$ ($\mu=\gamma$). In this case, when keeping terms of order $\epsilon$ and $\epsilon\ln(\epsilon)$ we get :

[TABLE]

Under the following approximation we get $\langle m\rangle(t)/\langle m\rangle^{*}\simeq 1-\mu te^{-\mu t}$ in agreement with the long time limit of equation (5).

6 Bibliography