Authors: Chen Jia, Ramon Grima
Categories: Article, Biological sciences, Mathematical biosciences, Systems biology
Source: iScience
The standard model describing the fluctuations of mRNA numbers in single cells is the telegraph model which includes synthesis and degradation of mRNA, and switching of the gene between active and inactive states. While commonly used, this model does not describe how fluctuations are influenced by the cell cycle phase, cellular growth and division, and other crucial aspects of cellular biology. Here, we derive the analytical time-dependent solution of an extended telegraph model that explicitly considers the doubling of gene copy numbers upon DNA replication, dependence of the mRNA synthesis rate on cellular volume, gene dosage compensation, partitioning of molecules during cell division, cell-cycle duration variability, and cell-size control strategies. Based on the time-dependent solution, we obtain the analytical distributions of transcript numbers for lineage and population measurements in steady-state growth and also find a linear relation between the Fano factor of mRNA fluctuations and cell volume fluctuations. We show that generally the lineage and population distributions in steady-state growth cannot be accurately approximated by the steady-state solution of extrinsic noise models, i.e. a telegraph model with parameters drawn from probability distributions. This is because the mRNA lifetime is often not small enough compared to the cell cycle duration to erase the memory of division and replication. Accurate approximations are possible when this memory is weak, e.g. for genes with bursty expression and for which there is sufficient gene dosage compensation when replication occurs.
**Subject ** Biological sciences, Mathematical biosciences, Systems biology
Experiments have revealed a large cell-to-cell variation in the number of mRNA molecules in isogenic populations.^1^^,^^2^^,^^3^ This can in part be explained by stochastic effects in gene expression due to the low copy numbers of many components, including DNA and important regulatory molecules.^4^ Live-cell imaging approaches allow a direct visualization of stochastic bursts of gene expression in living cells.^5^ However, these experiments are challenging and hence more commonly one measures the mRNA expression per cell from single-molecule fluorescence in situ hybridization^5^ or single-cell RNA-sequencing (scRNA-seq) experiments.^6^
The experimental distributions of mRNA numbers are fitted to the predictions of mathematical models, by which one can obtain estimates of the rates of several important transcriptional processes.^7^^,^^8^^,^^9^^,^^10^ The most common model of this type is the so-called two-state or random telegraph model of gene expression.^11^^,^^12^ This is composed of four (effective) reactions
where the first two reactions describe the switching of the gene between an active state G∗ and an inactive state G, the third reaction describes transcription while the gene is in the active state, and the fourth reaction describes the degradation of the mRNA M. The chemical master equation (CME) describing the telegraph model can be exactly solved in steady state, as well as in time.^12^^,^^13^^,^^14^^,^^15^ Extensions of this model to include more than two gene states have also been considered.^16^^,^^17^^,^^18^
A substantial number of genes are inactive most of the time and in the brief time that they are active, a large number of mRNA molecules are transcribed but not degraded.^19^ This leads to bursty expression. The probability of r new mRNA molecules being transcribed before the gene switches off, i.e. a burst of size r, is P(r)=pr(1−p), where p=ρ/(ρ+σ0) is the probability that the gene synthesizes an mRNA molecule, conditional on it being in the active state.^20^ This distribution is geometric with mean ρ/σ0. The average time between two consecutive bursts is 1/σ0+1/σ1≈1/σ1 since the gene spends most of its time off (σ0≫σ1); in other words, the rate of burst production is approximately σ1. It follows that the reaction scheme given in Equation 1 can be reduced to an effective one-state model composed of only two reactions
where k is the transcriptional burst size which is geometrically distributed with mean ρ/σ0. The geometric burst size distribution has been validated experimentally.^1^ The CME for this model can be solved exactly in steady state leading to the well-known negative binomial distribution of mRNA numbers,^21^^,^^22^ which is also widely used in scRNA-seq analysis.^23^ Because of the unimodality of this distribution, this simplified model cannot explain bimodality in gene expression,^24^^,^^25^ a feature that can be explained by the two-state model.
However, the conventional one-state and two-state models are very limited in their predictive power because they lack a description of many cellular processes that are known to have a profound impact on the distribution of mRNA numbers in single cells, e.g. the doubling of gene copy numbers upon DNA replication,^26^ partitioning of molecules during cell division,^27^ scaling of the mRNA synthesis rate with cell volume,^28^^,^^29^^,^^30^^,^^31^^,^^32^ and stochasticity in the cell cycle duration and growth rate that is related to cell-size control strategies.^33^^,^^34^^,^^35^^,^^36^^,^^37^^,^^38^^,^^39^ Recently, numerous efforts have been made to extend the conventional one-state and two-state models to include some description of these processes and yet retain analytical tractability. Some studies focused on the moment statistics (mean and variance) of mRNA and protein numbers,^40^^,^^41^^,^^42^^,^^43^^,^^44^ while other studies additionally obtain the analytical distributions of molecule numbers.^22^^,^^45^^,^^46^^,^^47^^,^^48^^,^^49^ Please refer to Table 1 for a summary of exactly solvable extensions of the one-state and two-state models that explicitly capture cell birth, growth, and division.
Due to mathematical complexity, most previous work is limited to the effective one-state model with the gene product (mRNA or protein) produced in a constitutive or bursty manner.^22^^,^^45^^,^^46^^,^^47^^,^^48^ Some of these models incorporate the scaling of transcription activity with cell volume,^22^^,^^48^ while the rest do not. We note that the latter case is not to be seen as unphysical since while the scaling of transcription with volume is commonly observed, it is by no means a universal phenomenon (in both prokaryotic^50^^,^^51^ and eukaryotic cells,^52^^,^^53^^,^^54^^,^^55^ there are examples where there is no such scaling). As for the conventional one-state model shown in Equation 2, the main limitation is the assumption of instantaneous bursts, while in reality there is a finite time for the bursts to occur. A distinct advantage of the extended one-state models over the conventional one is that those which describe gene replication^47^ are able to produce bimodal distributions.
The exact solution of extended two-state models that incorporate cell birth, growth, and division has not received much attention. A recent study^49^ made progress in this direction. In particular, the two-state telegraph model in growing and dividing cells was shown to be exactly solvable when (i) the mRNA synthesis rate scales linearly with cell volume and (ii) there is no variation of gene copy numbers across the cell cycle, i.e. gene replication is not taken into account. As previously mentioned, while (i) is common, it is not universal. The assumption behind (ii) is of course a means to simplify the model but clearly unrealistic. Relaxing any one of these two properties means that within the theoretical framework presented in Ref. 49, it is not possible to obtain an exact solution for the distribution of mRNA numbers.
While the aforementioned literature summarized in Table 1 has sought to fix the biological limitations of the conventional one-state and two-state models by directly introducing more processes and solving the master equation of the resulting complex models, a different indirect approach has also been proposed. This approach takes the point of view that biological processes not explicitly modeled by the conventional models can be incorporated by considering the model parameters themselves to vary between cells, and therefore to be drawn from probability distributions;^4^^,^^56^^,^^57^^,^^58^ we call this an extrinsic noise model (ENM). This model can be solved exactly in steady state for various distributions of parameter values (see Table I of Ref. ^58^). It is expected that such an approach produces meaningful results provided the parameters controlling cell-to-cell variability change very slowly. Under certain conditions, the solution of the ENM might even exactly match that of complex models of stochastic gene expression. For example, it has recently been shown that the exact solution of the two-state telegraph model in growing and dividing cells where gene replication is ignored and where the mRNA synthesis rate scales with cell volume is precisely the same as that of the ENM with the mRNA synthesis rate sampled from the distribution of cellular volume and with the mRNA degradation rate being replaced by an effective rate that also incorporates the dilution of molecules at cell division.^49^ A natural question is, if in a two-state telegraph model we introduce gene replication and allow the mRNA synthesis rate potentially to be volume dependent, then does the ENM still provide an exact or at least an accurate approximation of this model?
In this paper, we first exactly solve an extension of the telegraph model that explicitly describes cell birth, growth, division, replication, and an mRNA synthesis rate that can be either independent of cell volume or else that linearly scales with it. Many of the known exact solutions of the one-state and two-state models to-date can be shown to be special cases of the present theory. The analytical distribution of transcript numbers is subsequently used to study the accuracy of the ENM. We show that the transcript number distribution in steady-state growth is generally not well approximated by the steady-state distribution of the ENM. Conditions under which the ENM provides an accurate approximation are derived and verified using simulations.
We consider an extension of the telegraph model which takes into account cell growth, cell division, gene replication, gene dosage compensation, and volume-dependent transcription (see Figure 1 for an illustration). The specific meaning of all model parameters can be found in Table 2. The model has the following properties.
where σ0 and σ1 are the switching rates between the two gene states, and d is the mRNA degradation rate. For many genes in fission yeast,^28^^,^^29^ mammalian cells,^30^^,^^31^ and plant cells,^32^ there is evidence that the mRNA number scales linearly with cell volume in order to maintain approximately constant concentrations (concentration homeostasis; for a recent review see^61^). This is due to a coordination of the mRNA synthesis rate with cell volume—we shall refer to this mechanism as balanced mRNA synthesis. However, in both prokaryotic^50^^,^^51^ and eukaryotic cells,^52^^,^^53^^,^^54^^,^^55^ there are examples where there is no such scaling. Since each cell has a different volume, the mechanism of volume-dependent transcription is a source of extrinsic noise,^57^ potentially accounting for a significant amount of the observed cell-to-cell variation in mRNA numbers. To unify non-balanced and balanced mRNA synthesis, we assume that the mRNA synthesis rate depends on cell volume V(t) via a power law form with proportionality constant ρ and power β∈[0,1]. Then β=1 (β=0) corresponds to the situation where the mRNA synthesis rate scales linearly with cell volume (does not depend on cell volume). It has recently been postulated that the nonlinear scaling between gene expression levels and cellular volume is due to the heterogeneous recruitment abilities of promoters to RNA polymerases.^62^
Figure 1 ModelSchematic of an extension of the telegraph model of gene expression in growing and dividing cells. The volume V(t) of a cell grows exponentially with constant growth rate g and doubling time T. The gene expression dynamics is characterized by a two-state model with volume-dependent transcription and volume-independent degradation. Specifically, the gene can switch between an active state G∗ and an inactive state G. Transcription occurs when the gene is active. The synthesis rate of mRNA depends on cell volume V(t) via a power law form with power β∈[0,1], and the degradation rate of mRNA is a constant. Gene replication occurs at a time T0 where w=T0/T∈(0,1) is some fixed proportion of the cell cycle. Upon replication, the activation rate for each gene copy decreases from σ1 to σ1′ due to gene dosage compensation.
Here, we compute the time-dependent distribution of the mRNA number within a cell cycle under arbitrary initial conditions. We first consider the dynamics before replication for haploid cells. The microstate of the gene of interest can be described by an ordered pair (i,n), where i denotes the state of the gene with i=0,1 corresponding to the inactive and active states, respectively, and n denotes the number of mRNA molecules. Let pi,n(t) denote the probability of having n transcripts at time t∈[0,wT] when the gene is in state i. Note that t=0 corresponds to cell birth. Then, the stochastic gene expression dynamics before replication is governed by the coupled set of CMEs
where p1,−1=0 by default, the term involving ρ represents mRNA synthesis, the terms involving d represent mRNA degradation, and the terms involving σ0 and σ1 represent gene state switching, which will be referred to simply as gene switching in what follows. To solve them, we define a pair of generating functions Fi(t,z)=∑n=0∞pi,n(t)(z+1)n for i=0,1. Note that here we use (z+1)n rather than the conventional zn in the definition of the generating function—with this choice, the formulas given below are much more concise. In addition, let pn(t)=p0,n(t)+p1,n(t) denote the probability of having n transcripts at time t and let F(t,z)=F0(t,z)+F1(t,z) be the corresponding generating function. In terms of the generating functions, Equation 3 can be converted into the first-order linear partial differential equations (PDEs)
To solve them, we first convert them into a second-order parabolic PDE and then transform the second-order PDE into a hypergeometric differential equation through a change of variables. Complicated computations show that for each t∈[0,wT], the generating functions Fi, i=0,1 can be computed in closed form as (see STAR Methods for the proof)
Here Fi(0,z) and i=0,1 are the generating functions at t=0 which can be determined by the initial conditions, and the functions Kij, i,j=0,1 are given by
where the parameters r, a, b, and u are given by
Adding the two identities in Equation 5 gives the explicit expression of the generating function F before replication, i.e.
where the functions Li, i=0,1 are given by
When b=1, the term b−1 appears in the dominator of these equations and the equalities should be understood in the limiting sense. Note that when the mRNA synthesis rate is volume independent (β=0), the expression of F given in Equation 8 coincides with the time-dependent solution of the standard telegraph model.^14^
We next focus on the dynamics after replication for haploid cells. Since there are two daughter gene copies after replication, to distinguish them, we call them daughter copy A and daughter copy B. The dynamics of each gene copy is governed by the CMEs given in Equation 3 with σ1 being replaced by σ1′. Let pn(t) denote the probability of having n transcripts at time t∈[wT,T] and let F(t,z)=∑n=0∞pn(t)(z+1)n be the corresponding generating function. In STAR Methods, we prove that the generating function F after replication can be computed in closed form as
where L0′ and L1′ are functions obtained from L0 and L1 by replacing the parameters r, a, b, and u with
In summary, we have derived the analytical expression of the generating function F at any time t∈[0,T] within a cell cycle, which is given by
where Fi(wT,z) and i=0,1 are determined by Equation 5. The time-dependent distribution of the mRNA number can be recovered by taking the derivatives of the generating function F at z=−1, i.e.
Our analytical expression of the transient mRNA distribution is rather complicated. However, it can be greatly simplified in some special cases. In STAR Methods, we show how the analytical solution can be simplified for two non-trivial special (i) the gene switches rapidly between the active and inactive states (σ0,σ1≫g); (ii) the mRNA is produced in a bursty manner (σ0≫σ1), i.e. the gene is mostly inactive but transcribes a large number of mRNA when it becomes active.^72^^,^^73^^,^^74^^,^^75^ In the latter case, the burst frequency is σ1 before replication and the total burst frequency for the two gene copies is 2σ1′ after replication.
Thus far, we have obtained the transient mRNA distribution for haploid cells. For diploid cells, since the two alleles act independently and since each allele has the mRNA distribution given in Equation 11, the generating function for the total number of transcripts at any time t∈[0,T] is given by Fdiploid(t,z)=F(t,z)2, where F(t,z) is given by Equation 10. Here, we have used the fact that the generating function of two independent random variables is the product of their respective generating functions. Due to independence of the two alleles, when the rate parameters for each allele are fixed, the gene expression noise (measured by the coefficient of variation squared of mRNA numbers) in diploid cells is one half that in haploid cells. Without loss of generality, we always focus on haploid cells in what follows.
Thus far, we have derived the exact mRNA distribution at any time within a cell cycle. Here, we focus on the full time-dependence of the mRNA distribution across cell cycles under arbitrary initial conditions. To this end, we not only need the expression of F at any time t∈[0,T] but also need the expressions of Fi, i=0,1.
Recall that Equation 5 gives the analytical expressions of the generating functions Fi, i=0,1 before replication under any initial conditions. In particular, at replication, we have
Now, we focus on the dynamics of daughter copy A after replication. Let pi,n(t) denote the probability of having n transcripts at time t∈[wT,T] when the daughter copy A is in state i and let Fi(t,z)=∑n=0∞pi,n(t)(z+1)n be the corresponding generating function. In STAR Methods, we prove that the generating functions Fi, i=0,1 after replication can be computed exactly as
where Kij′ and i,j=0,1 are functions obtained from Kij by replacing the parameters r, a, b, and u with the parameters r′, a′, b′, and u′, respectively. Inserting Equation 12 into Equation 13 and taking t=T, we obtain
where
Suppose that the daughter cell with daughter copy A is tracked after division. Since we have assumed binomial partitioning of molecules at division, the probability pi,nnext(0) at birth in the next generation is given by
In terms of the generating function, the above relation can be written as
This gives the initial conditions for the next generation and the time-dependent mRNA distribution within the next cell cycle can be computed via Equation 10. Applying Equations 14 and 17 repeatedly, we are able to compute the full time-dependence of the Fi functions across cell cycles; substituting these in Equation 10 gives the full time-dependence of the mRNA distribution across cell cycles.
As a check of our analytical solutions, we compare the exact distributions of the mRNA number with the numerical ones obtained from a modified version of the finite-state projection (FSP) algorithm^76^ at three different time points (birth, replication, and division) across four cell cycles (Figure 2). In this algorithm, we couple the standard FSP with cell cycle events; for details see STAR Methods. Here, we assume that initially there are no mRNA molecules in the cell and the gene is off. This mimics the situation where the gene has been silenced by some repressor over a period of time such that all transcripts have been removed via degradation (while after silencing there may be some background level of mRNA, for simplicity we assume that all transcripts have been degraded). At time t=0, the repressor is removed and we study how gene expression recovers. When using FSP, we truncate the state space (to exclude states that are visited very rarely) and solve the associated truncated master equation numerically using the MATLAB function ODE45 with the dynamics before and after replication solved separately. Note that while the FSP and the stochastic simulation algorithm (SSA) yield comparable distributions of molecule numbers, the computational time of the former is much less than of the latter, provided the biochemical reaction networks are small enough—hence here we used the FSP. As expected, the analytical and simulated solutions coincide with each other completely at all times, and the mRNA distributions at birth, replication, and division reach a steady state within a few cell cycles. This can be also seen from Figure S1, where we illustrate the time-dependent mean and Fano factor of the mRNA number across four cell cycles.
Figure 2 Time-dependent mRNA distributions at birth, replication, and division across four cell cyclesThe blue curves show the analytical distributions computed by applying Equations 10, 14, and 17 repeatedly, and the red circles show the numerical ones obtained from FSP. The model parameters are chosen as Vb=1,g=1,β=1,w=0.4,d=5,ρ=20deff,σ0=1.5,σ1=3,σ1′=2.4 .
Another interesting observation is that the time-dependent mRNA distributions for our detailed telegraph model may exhibit three modes (Figure 2)—this is the combined effect of gene replication and slow switching between gene states. The zero mode is since there is no transcription when the gene is off while the two non-zero modes are due to transcription when the gene is turned on during the pre-replication and post-replication phases of the cell cycle. According to simulations, distributions with more than three modes are not observed in our model. This is different from the prediction of the conventional telegraph model^12^ whose distribution has at most two modes.
Thus far, we have obtained the full time-dependence of the mRNA distribution across cell cycles under arbitrary initial conditions. After several generations, the distribution at any fixed time within a cell cycle (such as the distributions at birth, replication, and division) becomes independent of the generation number. This is also called the cyclo-stationary condition in the literature^46^ or steady-state growth.^20^ Next, we compute the time-dependent mRNA distribution within a cell cycle under cyclo-stationary conditions.
Before computing the mRNA distribution, we first derive the probabilities of the gene being in the active and inactive states at any time within a cell cycle under cyclo-stationary conditions. Let pon(t) denote the probability of each gene copy being in the active state at time t∈[0,T]. Before replication, the dynamics of the active probability satisfies the differential equation p˙on=σ1(1−pon)−σ0pon. Solving this equation gives rise to
where we have used the fact that a/b=σ1/(σ0+σ1) and r+βg=σ0+σ1 (see Equation 7). Recall that the gene activation rate decreases from σ1 to σ1′ upon replication. After replication, the dynamics of the active probability satisfies the differential equation p˙on=σ1′(1−pon)−σ0pon. Solving this equation yields
where pon(wT) is determined by Equation 18. Combining Equations 18 and 19, we obtain the active probability of the gene at division, i.e.
Under cyclo-stationary conditions, the active probabilities at cell birth in two successive generations must be the same, i.e. ponss(0)=ponss(T). Then, the steady-state active probability of the gene at birth is given by
and thus the steady-state inactive probability at birth is given by poffb=1−ponb. It then follows from Equation 18 that the steady-state active probability of the gene at replication is given by
and thus the steady-state inactive probability at replication is given by poffr=1−ponr.
Next, we focus on the time-dependent mRNA distributions under cyclo-stationary conditions. Recall that we have obtained the time-dependent mRNA distributions within a cell cycle, whose generating function F(t,z) is given by Equation 10, provided that the initial conditions Fi(0,z), i=0,1 are known. Under cyclo-stationary conditions, the values of Fi(0,z) in two successive generations must be the same, i.e. Fi(0,z)=Finext(0,z), where Finext(0,z) has been derived in Equations 14 and 17. It then follows that the steady-state values of Fi(0,z) should satisfy
where
is a matrix-valued function with K˜ij, i,j=0,1 being given in Equation 15. Applying Equation 22 repeatedly, we obtain
Taking n→∞ in the above equation yields
where we have used the fact that
Once we have derived the steady-state values of Fi(0,z), i=0,1, it immediately follows from Equation 10 that the time-dependent generating function F under cyclo-stationary conditions is given by
Consider the case where gene replication is not taken into account (w=1) and when the mRNA synthesis rate scales with cell volume (β=1).^49^ In this case, the functions K˜ij(z) given in Equation 15 reduce to K˜ij(z)=Kij(T,z), and it is not difficult to see that Equation 22 can be solved analytically as
Inserting these equations into Equation 24 yields
For a given cell of volume V, its age is given by t=log(V/Vb)/g. Substituting t=log(V/Vb)/g in the above equation shows that the steady-state generating function for a cell of constant volume V is given by FV(z)=M(a;b;u˜z), where a=σ1/(d+g), b=(σ0+σ1)/(d+g), and u˜=ρV/(d+g). We make a crucial observation that this is exactly the steady-state generating function of the mRNA distribution for the conventional telegraph model.^12^
This result has been found in Ref. ^49^, which states that when w=β=1, the steady-state mRNA distribution for a cell of constant volume V of the detailed telegraph model is the same as that of the conventional telegraph model with effective decay rate deff=d+g. Note that the two terms in this rate capture the fact that transcripts are lost both by active degradation (with rate d) and by dilution at cell division (with rate g)—hence a model of this type is known as an effective dilution model (EDM).^77^ Intuitively, the EDM considers a population of cells with synchronized cell cycles so that at each time, all cells have the same volume.
Experiments have shown that in bacteria, most mRNAs have a half-life that is much shorter than the cell cycle duration, i.e. d≫g (see Table S1 for the typical values of d and g in various cell types), and thus are very unstable. The value of η=d/g can be used to measure the stability of mRNA. For unstable mRNAs (η≫1), the terms e−dt and e−d(t−wT) in Equation 10 are very small and thus can be approximated by zero (whenever t is not very close to 0 and wT). In this case, the time-dependent generating function F under cyclo-stationary conditions reduces to
where we have used the fact that Fi(0,0)=∑n=0∞pi,n(0) and Fi(wT,0)=∑n=0∞pi,n(wT) are the probabilities of the gene being in state i at birth and at replication, respectively. Imposing the term e−dt as zero in Equation 9 yields
When one of the gene switching rates σ0 and σ1 is very large, we have r=σ0+σ1−βg≫g and thus the second term on the right-hand side of Equation 27 can be neglected. This may occur when (i) the gene switches rapidly between the two states (σ0,σ1≫g), or (ii) the mRNA is produced in a constitutive manner (σ1≫σ0,g), or (iii) the mRNA is produced in a bursty manner (σ0≫σ1,g). In this case, the cyclo-stationary generating function Fss can be simplified significantly as
This contains much information. For a given cell of volume V<2wVb, its age is given by t=log(V/Vb)/g<wT and hence there is only one gene copy in the cell. Substituting t=log(V/Vb)/g in the above equation shows that the steady-state generating function for a cell of constant volume V is given by FV(z)≈M(a;b;u˜z), where
Here, we have used the fact that deff/(d+g)≈1 when mRNA is very unstable. Note that FV(z) is exactly the steady-state generating function of the mRNA distribution for the EDM
On the other hand, for a given cell of volume V>2wVb, its age is given by t=log(V/Vb)/g>wT and hence there are two gene copies in the cell. In this case, the EDM should be modified as
where GA and GB denote the two daughter copies whose dynamics are both governed by the conventional telegraph model. Substituting t=log(V/Vb)/g in Equation 28 shows that the steady-state generating function for a cell of constant volume V is given by FV(z)=M(a;b;u˜z)2, where a′≈σ1′/deff and b′≈(σ0+σ1′)/deff. Note that FV(z) is exactly the steady-state generating function of the mRNA distribution for the EDM given in Equation 30 since the two gene copies are independent of each other.
In summary, our analysis shows that for mRNAs with short lifetimes, the EDM makes a good approximation when one of the gene switching rates σ0 and σ1 is large (here the cell age t cannot be very close to 0 and wT, i.e. newborn cells and cells that have just finished gene replication should be excluded). This can be understood as follows. Previous studies^78^ have shown that the relaxation speed of the EDM to the steady state is governed by both the mRNA degradation rate d and the total gene switching rate σtot=σ0+σ1. When d and σtot are both large, any memory at birth from the previous cycle (due to binomial partitioning of molecules at division and to the gene state prior to division) and any memory at replication (due to gene state copying of the two daughter copies) will be rapidly erased. Each time that the volume changes, the mRNA distribution instantaneously equilibrates and hence the EDM works. Note that when the cell age t is close to 0 and wT, the memory at birth and at replication cannot be erased, which leads to the failure of the EDM. Relatively slow mRNA degradation and relative slow gene switching will both result in a deviation of the EDM from the full model.
In Figure 3, we compare the exact mRNA distributions with the numerical ones obtained from FSP at three different time points (birth, replication, and division) across the cell cycle under cyclo-stationary conditions. The truncated master equations are solved across several (usually less than five) cell cycles until the Hellinger distance between mRNA distributions at birth in two successive generations is less than 10−4. This guarantees that cyclo-stationary conditions are reached. When gene replication is not taken into account (w=1) and when the mRNA synthesis rate scales with cell size (β=1), the distributions of the full model agree perfectly with those of the EDM given in Equation 25 (Figure 3A). This coincides with our theoretical predictions. When gene replication is taken into account, the EDMs before and after replication are given by Equations 29 and 30, respectively. In this case, the EDM may deviate remarkably from the full model with the deviation being much larger at early stages of the cell cycle (Figure 3B), especially when mRNA degradation and gene switching are relatively slow. This can be understood as follows. According to the steady-state properties of the conventional telegraph model, in the presence of gene replication, the mean and the Fano factor of the mRNA number at birth for the EDM are given by
Figure 3 Comparison between the full model and the EDM(A) Steady-state mRNA distributions at birth, replication, and division for the full model and the EDM when gene replication is not taken into account. The blue curves show the analytical distributions given in Equations 23 and 24, the red circles show the numerical ones obtained from FSP, and the gray regions show the distributions of the EDM.(B) Same as (A) but when gene replication is taken into account. In (A) and (B), the model parameters are chosen as Vb=1,g=1,β=1,d=4,ρ=20deff,σ0=1.5,σ1=3,σ1′=2.4 . The parameter w is chosen as w=1 in (A) and w=0.4 in (B).(C) Same as (B) but in the special case where mRNA synthesis is balanced and bursty, and dosage compensation is perfect. The model parameters are chosen as Vb=1,g=1,β=1,w=0.4,d=4,ρ=200deff,σ0=300,σ1=30,σ1′=15.
and the mean and the Fano factor of the mRNA number at division are given by
Under cyclo-stationary conditions, it follows from Equation 16 that the mean mRNA numbers at birth and at division for the full model should satisfy ⟨n⟩(T)=2⟨n⟩(0) and FanoT=2Fano0−1. However, these two restrictions in general do not hold for the EDM—the EDM satisfies these two restrictions only when
Note that when mRNA synthesis is balanced (β=1) and bursty (σ0≫σ1), the above restrictions are satisfied when dosage compensation is perfect (σ1′=σ1/2), i.e. when the total burst frequency does not change when replication occurs. When these three conditions are satisfied, the EDM makes accurate predictions and the mRNA number follows a negative binomial distribution (Figure 3C). The breakdown of the above restrictions will give rise to the deviation of the EDM from the full model, as observed in Figure 3B. Intuitively, this is because the mRNA distribution at birth is affected by the fluctuations of the two gene copies at division and thus in general it cannot be captured solely by an EDM with only one gene copy. Note that special case 2 discussed above may not satisfy the above moment equalities since in this special case, the EDM fails for newborn cells.
We next compute the steady-state distributions of transcript numbers measured over a cell lineage or from a population snapshot . In lineage measurements, the mRNA number from an individual cell is tracked at any point in time, i.e. once the cell divides, only one of the two daughter cells is tracked. Clearly, the probability of observing a cell of age t∈[0,T] is 1/T for lineage measurements. As a result, the generating function of the steady-state distribution along a cell lineage is given by
In contrast, in population measurements, the mRNA numbers in a population of isogenic cells are observed at a particular time. Previous studies^20^ have shown that the probability of observing a cell of age t∈[0,T] is 2(1−t/T)(log2)/T=2ge−gt for population measurements. Thus, the generating function of the steady-state distribution in a population of cells is given by
Our analytical expression of the steady-state distribution is rather complicated since we have to integrate the time-dependent distribution over time which involves complex confluent hypergeometric functions. However, it can be simplified to a large extent in some special cases. In STAR Methods, we show how the analytical solution can be simplified in two non-trivial special (i) the mRNA is unstable and the gene switches rapidly between the two states; (ii) the mRNA is unstable and the gene switches slowly between the two states. In particular, in the latter case, the steady-state distribution for lineage measurements is given by
and the steady-state distribution for population measurements is given by
where Γ(n,λ)=∫λ∞tn−1e−tdt is the incomplete gamma function and δ0(n) is the Kronecker delta which takes the value of 1 when n=0 and takes the value of 0 otherwise. Interestingly, for the two types of measurements, the mRNA distribution has a zero-inflated part (the first part). Indeed, numerous biostatistical papers^79^^,^^80^^,^^81^ have used zero-inflated models to characterize mRNA distributions in scRNA-seq data analysis. For population data such as in scRNA-seq, our theory predicts that the probability of true zeros, i.e. dropout events that are not due to purely technical reasons, is p0pop=(2−21−w)poffb+(21−w−1)poffr. This is important for interpreting scRNA-seq data and may be potentially useful for data imputation.
Figures 4A and 4B compare the population mRNA distribution given in Equation 34 and the numerical one obtained from FSP. Clearly, the two distributions agree perfectly when the transcripts are highly unstable and when the gene switching rates are very small. Particularly, we find that the mRNA distribution is capable of exhibiting three modes, corresponding to the three terms in Equation 34 (left panel of Figure 4A). This shows that a trimodal distribution may occur in the special case of unstable mRNA and slow gene switching. When the transcripts are relatively stable, trimodality becomes less apparent and Equation 34 fails to capture the real mRNA distribution, while it can still well capture the probability of zero observations (right panel of Figure 4A).
Figure 4 Distribution and moment analysis for population measurements(A) Steady-state mRNA distribution in a cell population when gene switching is very slow. The blue curves show the analytical distributions given in Equation 34 and the red circles show the numerical ones obtained from FSP. The model parameters are chosen as Vb=1,g=1,β=1,w=0.3,ρ=8deff,σ0=2×10−4,σ1=10−3,σ1′=9×10−4 . The mRNA degradation rate is chosen as d=50 in the left panel and d=10 in the right panel.(B) Fano factor of the mRNA number versus the Fano factor of cell volume in a cell population. The blue line is computed from the analytical expressions given in Equations 35 and 37 and the red circles are obtained from FSP. The model parameters are chosen to be the same as in the left panel of (A).
We next analyze the moments of transcript numbers. Here, we consider the general case without making any simplifying assumptions. The mean and the second factorial moment of the mRNA number at any time t∈[0,T] within a cell cycle can be recovered by taking the derivatives of the generating function F at z=0, i.e.
Straightforward computations show that
where ⟨n⟩(0) is the mean of the initial mRNA number, λ=au/b, λ′=a′u′/b′, and
Under cyclo-stationary conditions, the mean at division should be twice that at birth, i.e. ⟨n⟩(T)=2⟨n⟩(0). This shows that the steady-state mean at birth is given by
where η=d/g. Inserting this equation into Equation 35 gives the steady-state mean at any time within a cell cycle.
The explicit expression for the second moment is extremely complicated since we need to take the second derivative of a complex generating function. However, when the mRNA is unstable, the generating function has a relative simple expression (see Equation 26) and taking the second derivative of this function yields the second factorial moment of the mRNA number at any time within a cell
where
The steady-state moments of the mRNA number for lineage and population measurements can then be obtained by integrating Equations 35 and 37 over time. For example, the mean and second factorial moment of the mRNA number for population measurements are given by
These explicit expressions can be obtained but are omitted here since they are too complicated. In STAR Methods and Figure S2, we find that the lineage mean is always greater than the population mean, and the difference between them is at most 10%.
A crucial observation made from the analytical results is that for both types of measurements, the Fano factor of mRNA number fluctuations, FanomRNA, and the Fano factor of cell volume fluctuations, Fanovolume, must satisfy the following relation when mRNA synthesis is balanced (see STAR Methods for the proof):
where C is a constant independent of the birth volume Vb and the growth rate g. In Figure 4B, we validate this relation using both the exact solution and FSP. Our result shows that the fluctuations in gene expression and cell volume, characterized by the Fano factors, are linearly correlated when the mRNA synthesis rate scales with cell size. This may be potentially useful for checking whether mRNA synthesis is balanced in living organisms.
In particular, when the mRNA is unstable and when gene expression is bursty, the Fano factor of mRNA number fluctuations can be computed explicitly by integrating Equations 35 and 37 over
When dosage compensation is perfect, we have 2a′=a and thus the above equation reduces to
where we have used the fact that 2(log2)2≈1. Note that when gene expression is bursty, the mean burst size at time t is ρV(t)/σ0 and hence the mean burst size over the whole cell cycle is given by B=(ρVb)/[(log2)σ0], which can be obtained by averaging ρV(t)/σ0 over time. In this case, we have FanomRNA≈1+B, which reduces to the well-known result for the conventional telegraph model.^72^
Our detailed telegraph model involves the coupling between gene expression dynamics, cell volume dynamics, and cell cycle events. In STAR Methods and Figure S3, we show that the steady-state distribution of the detailed model cannot be captured by the steady-state solution of the conventional telegraph model given in Equation 1 with volume-independent rates, even when gene replication is not taken into account (w=1). In previous studies, the lineage and population distributions for the detailed model have often been approximated by the distributions for the ENM.^49^ In the ENM, the mRNA distribution for a cell of constant volume V is exactly the one predicted by the EDM, and the fluctuations of cell volume V are regarded as extrinsic noise.^57^^,^^58^ In other words, the mRNA distribution for the ENM is given by
where Π(V) is the distribution of cell volume. We emphasize here that the EDM varies depending on the number of gene copies and thus also depending on cell volume. For a cell of volume V<2wVb, there is only one gene copy and the EDM is given by Equation 29; for a cell of volume V≥2wVb, there are two gene copies and the EDM is given by Equation 30. In addition, note that the distribution of cell volume is different for lineage and population measurements. Since cell volume V(t) and cell age t are related by V(t)=Vbegt, the cell volume distribution can be obtained from the cell age distribution which has already been given above (see the paragraphs before Equations 31 and Equation 32). Specifically, the volume distribution for lineage measurements is given by^70^
and the volume distribution for population measurements is given by
Inserting the above two equations into Equation 39 gives the mRNA distribution for the ENM.
To evaluate the performance of the ENM approximation, we first illustrate the Hellinger distance D between the lineage distributions of the full model and the ENM as a function of σ0/σ1 and σ1′/σ1 when mRNA synthesis is balanced, i.e. β=1 (Figure 5A). It can be seen that the ENM serves as a good approximation when gene expression is bursty (σ0≫σ1) and when dosage compensation is perfect (σ1′=σ1/2). This is indeed a sufficient condition for mRNA to display concentration homeostasis when gene replication is taken into account.^48^ A proof of this condition can be found in STAR Methods. The breaking of either dosage compensation or bursty expression will give rise to a significant deviation of the ENM from the full model (Figure 5B). In particular, the distribution of the ENM can show bimodality whereas that of the full model is unimodal.
Figure 5 Comparison between the full model and the ENM(A) Heat plot of the Hellinger distance D between lineage distributions for the full model and the ENM as σ0/σ1 and σ1′/σ1 vary. The model parameters are chosen as Vb=1,g=1,β=1,w=0.5,d=5,σ1=30 .(B) Comparison between the lineage distributions for the full model and the ENM as σ0/σ1 and σ1′/σ1 vary. The blue curves show the analytical distributions for the full model given in Equation 31, the red circles show the numerical ones obtained from FSP, and the gray regions show the distributions for the ENM. The model parameters are chosen as in (A). The parameter σ0 is chosen as σ0=10σ1 (bursty case) and σ0=0.5σ1 (non-bursty case). The parameter σ1′ is chosen as σ1′=σ1/2 (perfect dosage compensation) and σ0=σ1 (no dosage compensation). The parameters associated with the four panels are marked in (A) by stars.(C) Heat plot of D as β and σ1′/σ1 vary. The model parameters are chosen as Vb=1,g=1,w=0.5,d=5,σ0=300,σ1=30 .(D) Heat plot of D as η and σtot/g vary. The model parameters are chosen as Vb=1,g=1,β=1,w=0.5,σ0=2.5σ1,σ1′=σ1.(E) The model parameters are chosen to be the same as in the third panel of (B) but η is varied. In (A)–(E) the parameter ρ is chosen so that ⟨n⟩lin=30.
It is still unclear how the ENM performs when mRNA synthesis is not balanced (β<1). To see this, we further illustrate D as a function of β and σ1′/σ1 when gene expression is bursty (Figure 5C). Interestingly, there is a region of parameter space (shown in dark blue) where D is minimized. In particular, when the mRNA synthesis rate is volume independent (β=0), the ENM works well when σ1′/σ1 is between 0.65 and 0.8. This shows that to maintain the effectiveness of the ENM approximation, a lack of balanced mRNA synthesis requires also a lower degree of dosage compensation. Recent studies have shown that even when β<1, strong concentration homeostasis (characterized by a small coefficient of variation of the mean concentrations across the cell cycle) can still be obtained when σ1′/σ1≈1/2β+1 (shown by the yellow dashed line) and when replication occurs halfway through the cell cycle (w=0.5).^48^ Note that the region where D is minimized is exactly around the yellow dashed line. This shows that the effectiveness of the ENM approximation is closely related to concentration homeostasis even when β<1.
To further confirm our results, we use the transcriptional parameters inferred in Ref. 26. In this case, the mRNA distributions for two bursty genes Oct4 and Nanog in mouse embryonic stem cells were measured as a function of time in the cell cycle from which all the rate parameters involved in our model were estimated. Since the cell-to-cell variability in volume within each cell cycle phase was quite small, it was assumed that β=0, i.e. the mRNA synthesis rate is volume independent. Dosage compensation was found to be apparent for both genes, with σ1′/σ1 estimated to be 0.63 for Oct4 and 0.71 for Nanog. Based on the inferred parameters, we compare the mRNA distributions of population measurements for the full model and the ENM (Figure S4A). We find that the ENM performs well for both genes. This agrees with our prediction that the ENM is valid when β=0 and when σ1′/σ1 is around 0.7 (Figure 5C). However, if we keep all rate parameters the same but reset σ1′/σ1 to 1 (no dosage compensation), then the EDM approximation will become significantly less satisfying (Figure S4B). This also coincides with the simulations shown in Figure 5C.
When mRNA synthesis is balanced and bursty, we have seen that the ENM approximation is accurate when dosage compensation is strong. However, in bacteria and budding yeast, there has been some evidence that dosage compensation is not widespread.^50^^,^^82^ It is unclear under what conditions the ENM is still valid when dosage compensation is weak. To see this, we also depict D as a function of η=d/g and σtot/g when there is no dosage compensation, i.e. σ1′=σ1 (Figure 5D). In this case, we find that the ENM still works well when the mRNA is very unstable (d≫g) and when the total gene switching rate is very large (σtot≫g). This is fully consistent with our earlier theoretical predictions for the accuracy of the EDM, on which the ENM depends. In particular, when gene expression is bursty (σ0≫σ1,g), increasing the mRNA degradation rate will give rise to a better ENM approximation (Figure 5E). Note that this is not true when the total gene switching rate is slow (Figure 5D). We emphasize that while Figures 5B and 5E show the mRNA distributions for lineage measurements, the same results are applicable for population measurements (Figure S5).
The value of η=d/g can be determined experimentally since both d and g can be measured. In bacteria, η is typically between 6−30, depending strongly on the strain and the growth condition; in yeast, it is typically between 3−8; and in mammalian cells, it is typically between 2−4 (see Table S1 for the median and range of η in various cell types). This suggests that the ENM approximation may be generally most useful in bacteria and less useful in yeast and mammalian cells.
Thus far, we have considered a detailed telegraph model of gene expression with a cell cycle description when the cell volume dynamics and the cell cycle duration are deterministic. However, in naturally occurring systems, the cell cycle duration is appreciably stochastic (see Figure 1C of Ref. ^47^ for experimental distributions of cell cycle durations in eight different cell types). Moreover, there has been ample evidence^33^^,^^34^^,^^35^^,^^36^^,^^37^^,^^38^^,^^39^ that the amount of growth produced during the cell cycle must be controlled such that, on average, larger cells at birth have shorter cell cycle durations than smaller ones. This mechanism maintains size homeostasis.
To model cell-cycle duration variability and size homeostasis, we use the size-additive autoregressive model of stochastic cell volume dynamics.^36^^,^^83^ The model assumes that the volume at birth Vb and the volume at division Vd are connected by the relation
where 0≤α≤2 is the strength of size control, v¯>0 is the typical (average over generations) birth volume which is a time-independent constant, and ϵ∼N(0,σϵ2) is a Gaussian noise term independent of Vb. The idea behind the model is as upon being born with volume Vb, the cell attempts to grow for a period of time such that its target volume at division is f(Vb)=αVb+(2−α)v¯, but due to stochasticity, the actual volume at division may deviate from the target volume. Due to exponential cell growth, the cell cycle duration T is given by
where for simplicity we have assumed constant growth rate g across generations. This implies that on average, larger cells at birth have shorter cell cycle durations than smaller ones. Different size control strategies correspond to different values of α. When α=0, the target division volume f(Vb)=2v¯ is constant; this corresponds to the sizer strategy, where cells have to reach a certain size before division occurs. When α=1, the cell attempts to add a constant volume f(Vb)−Vb=v¯ to its newborn size; this corresponds to the adder strategy. Since the growth is exponential, attempting to grow for a constant time is the same as having f(Vb)=2Vb; hence α=2 corresponds to the timer strategy. The adder or near-adder behavior has been observed in bacteria, budding yeast, and mammalian cells,^35^^,^^37^^,^^39^ while fission yeast exhibits a near-sizer behavior.^33^
When σϵ=0, the model reduces to deterministic (previously considered) cell volume dynamics, in which case the timer, adder, and sizer strategies are exactly the same since Vb=v¯ is a constant. As σϵ increases, the time series of cell volume becomes much more noisy; however, it is difficult to identify whether there is a change in the magnitude of fluctuations solely from the time series of the mRNA number (Figure 6A). Note that when σϵ is small, the model produces a steady-state cell size distribution (from lineage simulations) characterized by three a fast increase in the size count for small cells, a slow decay for moderately large cells, and a fast decay for large cells (Figure 6B). This is consistent with the cell size distribution in E. coli.^70^ A natural question is what are the values of σϵ in naturally occurring systems. To see this, we examined the publicly available lineage data of cell size in E. coli and fission yeast^36^^,^^84^ and found that the typical value of σϵ is between 0.2v¯ and 0.3v¯ (see STAR Methods for a discussion about the inference of σϵ and the estimated values of σϵ in E. coli and fission yeast under different growth conditions).
Figure 6 Effects of stochastic cell volume dynamics on mRNA fluctuations(A) Typical trajectories of cell size and mRNA number as σϵ increases.(B) Cell volume distribution of lineage measurements as σϵ increases. In (A) and (B), the model parameters are chosen as v¯=1,g=1,β=1,w=0.4,d=4,ρ=20deff,σ0=1.5,σ1=3,σ1′=2.4,α=1.(C) Comparison between the steady-state mRNA distributions of lineage measurements for deterministic and stochastic cell size dynamics under different size control strategies. The blue curves show the analytical distributions for stochastic cell size dynamics, the red circles show the numerical ones obtained from the SSA, and the gray regions show the distributions for deterministic cell size dynamics. The model parameters are chosen as v¯=1,g=1,β=1,w=0.5,d=5,σ0=σ1=100,σ1′=50,σϵ=0.4. The parameter ρ is chosen so that ⟨n⟩lin=10 for deterministic cell size dynamics. Previous studies^93^ have shown that the timer strategy with α=2 is not stable since it cannot produce a finite and nonzero mean of cell volume. Hence, we choose α=1.8 for the timer strategy here.(D) Steady-state mRNA distributions at birth, replication, and division as σϵ increases. The model parameters are the same as in (A) and (B).(E) Comparison between the lineage distributions of the full model and the ENM for stochastic cell size dynamics. The gray regions show the distributions for the ENM. The model parameters are chosen to be the same as in Figure 5B.
To compute the mRNA distribution for stochastic cell volume dynamics, note that the evolution of the system within a cell cycle is controlled by four random (i) the gene state at birth αb, (ii) the mRNA number at birth Nb, (iii) the birth volume Vb, and (iv) the cell cycle duration T. Once the values of the four variables are fixed, the generating function F at any time t∈[0,T] within a cell cycle is given by Equation 10, i.e.
Here, the initial conditions Fi(0,z), i=0,1 are determined by αb and Nb as
The functions Li and Li′, i=0,1 given in Equation 9 depend on Vb since the parameters u and u′ are functions of Vb; the replication time wT depends on T. Hence, the generating function F depends on all the four variables. Once we know the joint distribution of the four variables in some generation, we can use Equation 42 to compute their joint distribution in the next generation. In this way, we obtain the full time-dependence of the mRNA distribution cross cell cycles. In STAR Methods, we have generalized the analytical results obtained previously to the model with stochastic cell volume dynamics. Specifically, we have derived the exact time-dependent mRNA distribution for a cell of any age in any generation, as well as the exact steady-state distribution for lineage measurements.
To reveal the influence of cell-cycle duration variability and size homeostasis on gene expression, we compare the lineage distributions for the model with deterministic cell size dynamics and the model with stochastic cell size dynamics under different size control strategies (Figure 6C). The two distributions deviate remarkably from each other for the timer strategy, but the deviation is much smaller for the adder and sizer strategies. This demonstrates the advantage of the adder and sizer strategies in reducing gene expression noise. In addition, Figure 6D illustrates the steady-state mRNA distribution at three different points (birth, replication, and division) across the cell cycle as noise in cell size dynamics, characterized by σϵ, varies. Clearly, larger noise in cell size results in larger noise in gene expression, as expected. A sufficiently large σϵ may even change the number of modes of the mRNA distribution. Interestingly, we find that when the mRNA distribution exhibits multimodality, increasing σϵ will not change the height of the zero peak but may affect the height and position of non-zero peaks (Figure 6D).
Finally, we investigate the accuracy of the ENM approximation for stochastic cell volume dynamics. Note that we can no longer use the EDM to approximate the mRNA distributions at birth, replication, and division, since the cell volumes are stochastic. We compare the steady-state mRNA distributions at birth, replication, and division for the full model with their ENM approximations in Figure S6 and also compare the lineage distribution for the full model with its ENM approximation in Figure 6E (see STAR Methods for the analytical expressions of the ENM approximations). The model parameters in the two figures are chosen to be the same as in Figures 3B and 3C and 5B, respectively. We can see that in the presence of fluctuations in cell volume, the results of the present paper are still valid—the ENM does not work in general but performs well when mRNA synthesis is balanced and bursty and when dosage compensation is perfect. Comparing Figure S6 with Figures 3B and 3C and comparing Figure 6E with Figure 5B, we also find that the differences between the mRNA distributions for the full model and the ENM are slightly diminished when the cell volume dynamics is stochastic.
In this work, we analytically solved a detailed model of stochastic gene expression with cell cycle and cell volume descriptions including gene switching, cell growth, cell division, volume-dependent transcription, gene replication, and gene dosage compensation. We first considered the case where the cell volume dynamics is deterministic and then generalized the results to include cell-cycle duration variability and cell-size control strategies. Previous models of stochastic mRNA dynamics in growing and dividing cells^22^^,^^47^ can be seen as special cases of the present modeling framework. For example, when mRNA synthesis scales with cell volume and when the gene inactivation rate is much higher than the gene activation rate, our model reduces to the one studied in Ref. 22. Under this timescale separation assumption, there is essentially only one gene state and the computation is much easier than the one given in the present paper. If the intrinsic noise due to the random birth-death of transcripts is ignored, then our model reduces to the one-state model studied in Ref. 40. In addition, we emphasize that our model not only characterizes the mRNA dynamics but can also be used to describe the protein dynamics. For example, when gene expression is bursty and when the degradation rate is taken to be zero, our model reduces to the effective one-state model of the protein dynamics proposed in Refs. ^41^^,^^42^^,^^46^.
Our work is also distinctive from recent related work^49^ since our derivations of the distributions of mRNA numbers as a function of cell age and generation number, and of the distributions in steady-state growth do not need the assumption of stochastic concentration homeostasis (SCH); the relaxation of this assumption is crucial to model the variation of gene copy numbers across a cell cycle due to DNA replication. We have also investigated how well can the model be approximated by the effective dilution and extrinsic noise models (EDM/ENM). When gene replication is taken into account, we showed that the mRNA distributions of the full model may differ significantly from the predictions of the EDM/ENM. We elucidated three cases where the EDM/ENM makes accurate approximations.
The first case occurs when the mRNA is very unstable and the total gene switching rate (the sum of the gene activation and inactivation rates) is very large such that on the timescale of volume change, the mRNA distribution instantaneously equilibrates. This condition is intuitive and has been discussed in earlier work.^57^ However, as we showed using data from various cell types, the typical mRNA lifetime in eukaryotes (especially mammalian cells) is generally not small enough compared to the cell cycle duration to enforce instantaneous equilibrium; rather, the fluctuations have memory of birth and replication events.
The second case takes place when mRNA synthesis is balanced and bursty and when dosage compensation is perfect. While our model does not generally obey SCH due to gene copy number variation upon replication, however, in this case, parameter conditions effectively enforce SCH. Note that if expression is balanced and it is bursty with weak dosage compensation or else it is constitutive with perfect dosage compensation, there is an apparent breakdown of the EDM/ENM’s ability to accurately approximate the full model. This is since in these cases the dependence of the mean mRNA numbers with cellular volume is significantly influenced by the doubling of gene copy numbers at replication. Examples where expression is balanced but the effects of replication are not completely buffered by dosage compensation are starting to be uncovered, e.g. in human cells while the overall mRNA synthesis rates increase with cell volume, however, S/G2-phase cells show increased synthesis rates compared to G1-phase cells of the same volume.^85^ As pointed out in Ref. ^61^, this is reminiscent of a step-increase in RNA production during or after S phase which was previously observed in synchronized HeLa cell populations and other organisms^86^—this suggests that perfect dosage compensation in mammalian cells may not be common.
The third case is when mRNA synthesis is non-balanced and bursty, and when dosage compensation is of an intermediate strength such that concentration homeostasis is approximately maintained, i.e. there is only a small variation of the mean mRNA concentration throughout the cell cycle—note that this is a much weaker condition than SCH. We showed that this is indeed the case for two genes Oct4 and Nanog in mouse embryonic stem cells, whose parameters have been previously estimated before and after gene replication.^26^ Recent studies^48^ have shown, using both theory and data, that when gene expression is bursty, deviations from SCH show up as deviations from the gamma distribution in the mRNA concentration. This can be used to test whether SCH is approximately valid in vivo.
Our model is complex due to the coupling between gene expression dynamics, cell volume dynamics, and cell cycle events. A natural question is whether all the parameters involved in the model can be inferred accurately. In fact, parameter inference for models that are more complex than the telegraph model but simpler than our model has been made in our previous papers using the method of distribution matching^48^ or power spectrum matching.^66^ Whether accurate parameter estimation is possible by fitting mRNA distributions from population snapshot data to the analytical expression given by our calculations remains an open question.
In summary, our work shows that caution is needed when the ENM is applied to explain data collected in growing and dividing cells and that the accuracy of this reduced model of gene expression cannot be a priori assumed genome-wide. Our model, though detailed, has some limitations. We have focused on models that explain cell-to-cell variability in the synthesis rates due to their dependence on cell size. However, likely other descriptors of cell state (such as shape, local cell crowding, mitochondrial abundance, and capacity to respond to Ca^2+^) can explain a higher degree of cell-to-cell variability than cell size alone.^87^^,^^88^ In addition, here we have not considered the G0 phase, where cells are not growing and are outside of the replicative cell cycle. For the subpopulation of cells permanently in G0 such as senescent and many differentiated cells, since cells do not grow and divide, we can always use the ENM to characterize their gene expression dynamics.
Last but not least, here we have considered the expression of unregulated genes but it is well known that many genes regulate each other resulting in complex gene regulatory networks.^89^ Overcoming the last limitation is particularly pressing but it is analytically challenging because such models have nonlinear propensities stemming from the modeling of bimolecular interactions between transcriptional factors and genes.^90^ Progress in this direction will be reported in a separate paper. We also anticipate that the results of the present paper can be generalized to include more than two gene states.^16^^,^^17^^,^^18^
Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Ramon Grima (ramon.grima@ed.ac.uk).
This study did not generate new unique reagents.
Here we derive the analytical expression of the generating function F before replication. Adding the two identities in Equation 4 shows that F1 can be represented by F as
Inserting this equation into the second equation of (4) shows that F satisfies the second-order parabolic PDE
where r=σ0+σ1−βg. Following Ref. ^14^, we introduce a new variable τ=logz−dt. Let F˜(τ,z) and F˜i(τ,z) be the functions with variables τ and z that are associated with F(t,z) and Fi(t,z), respectively, i.e.
Then Equation 47 can be simplified to a large extent as
If we fix the variable τ, this is an ordinary differential equation (ODE) with respect to the variable z. Note that
Inserting this equation into Equation 48 yields
This is a modified version of the confluent hypergeometric differential function and its solution can be written in general form as
where
with a=σ1/(d+βg), b=(σ0+σ1)/(d+βg), u=ρVbβ/(d+βg), and with M(a;b;x) being the confluent hypergeometric function.
The remaining question is how to determine the functions φ0 and φ1 based on the initial conditions. By the definition of F˜i, it is easy to see that
Taking τ=logz in Equation 50 yields
where
Moreover, it follows from Equation 46 that
This shows that
It follows from Equation 50 that
Combining the above two equations yields
where
Combining Equations 52 and 53, we obtain
This shows that
By means of the Wronskian of confluent hypergeometric functions, it is easy to check that
Moreover, straightforward computations show that
Inserting the above two equations into Equation 54, we obtain
In terms of the original variables t and z, it follows from Equation 50 that the generating function F is given by
Interesting Equations 51 and 55 into the above equation, we finally obtain the analytical expression of the generating function F, which is given by
Here F0(0,z) and F1(0,z) are the generating functions at t=0 which can be determined by the initial conditions, and the functions L0 and L1 are given by
Here we derive the analytical expression of the generating functions F0 and F1 before replication. From the first equation of (4), F1 can be represented by F0 as
Inserting into the second equation of (4) shows that F0 satisfies the second-order parabolic PDE
where R=σ0+σ1+d. Using the new variables τ and z, Equation 58 can be simplified to
If we fix the variable τ, this is exactly an ODE with respect to the variable z. Inserting Equation 49 into Equation 59 yields
This is a modified version of the confluent hypergeometric differential function and its solution can be written in general form as
where
The remaining question is how to determine the functions φ0 and φ1 based on the initial conditions. By the definition of F˜i, it is easy to see that
Taking τ=logz in Equation 60 yields
where
Moreover, it follows from Equation 57 that
This shows that
It follows from Equation 60 that
Combining the above two equations yields
where
Combining Equations 62 and 63, we obtain
This shows that
By means of the Wronskian of confluent hypergeometric functions, it is easy to check that
Moreover, straightforward computations show that
Inserting the above two equations into (64), we obtain
In terms of the original variables t and z, it follows from Equation 60 that the generating function F0 is given by
Inserting Equations 61 and 65 into the above equation, we finally obtain the analytical expression of the generating function F0, which is given by
Here F0(0,z) and F1(0,z) are the generating functions at t=0 which can be determined by the initial conditions, and the functions K00 and K01 are given by
where we have used the fact that R=r+d+βg. Since we have derived both F0 and F, we finally obtain the explicit expression of F1=F−F0, which is given by
where
We next focus on the dynamics after replication for haploid cells. Since there are two daughter gene copies after replication, to distinguish them, we call them daughter copy A and daughter copy B. For convenience, we assume that all mRNA molecules are allocated to daughter copy A when replication occurs, while no molecules are allocated to daughter copy B. Similarly, the microstate of each daughter copy can be described by an ordered pair (i,n), where i=0,1 denotes the gene state and n denotes the number of mRNA molecules belonging to that daughter copy, i.e. the molecules allocated to that daughter copy at replication and the molecules produced by that daughter copy after replication. The total number of molecules after replication is then the sum of the numbers of molecules that belong to the two daughter copies. Note that other choices for how the mRNA molecules are allocated to each of the gene copies have no effect on the calculation of the statistics of the total number of molecules after replication.
Let pi,nA(t) denote the probability of having n transcripts that belong to daughter copy A at time t∈[wT,T] when daughter copy A is in state i. Similarly, let pi,nB(t) denote the same quantity for daughter copy B. Recall that the gene activation rate decreases from σ1 to σ1′ upon replication. Then the stochastic gene expression dynamics for each gene copy after replication is governed by the CMEs
where l=A,B indicates which daughter copy is considered. To solve these, for each l=A,B, we define a pair of generating functions
In addition, let pnl(t)=p0,nl(t)+p1,nl(t) denote the probability of having n transcripts that belong to daughter copy l at time t and let Fl(t,z)=F0l(t,z)+F1l(t,z) be the corresponding generating function. Then Equation 68 can be converted into the PDEs
For any t∈[0,wT], let α(t) denote the state of the mother copy and let X(t) denote the number of mRNA molecules at time t. For any t∈[wT,T] and j=A,B, let αj(t) denote the state of daughter copy j and let Xj(t) denote the number of mRNA molecules that belong to daughter copy j at time t. Since the two daughter copies inherit the gene state of the mother copy at replication, we have
Since all molecules are allocated to daughter copy A and no molecules are allocated to daughter copy B at replication, we have
Since αA(wT)=αB(wT), the initial distributions for the two daughter copies are correlated and thus the dynamics for the two daughter copies after replication are not independent of each other. However, once the gene state of the mother copy at replication is fixed (conditioned on α(wT)=k with k=0,1), the initial distributions for the two daughter copies are (conditionally) independent of each other, and hence the dynamics for the two daughter copies after replication are also (conditionally) independent of each other. We now use the conditional independence of the two daughter copies to compute the generating function F after replication.
Recall that F(t,z) is the generating function of
which denotes the probability of having n mRNA molecules in the cell at time t. Moreover, recall that FA(t,z) is the generating function of
which denotes the probability of having n mRNA molecules that belong to daughter copy A at time t. In addition, recall that FB(t,z) is the generating function of
which denotes the probability of having n mRNA molecules that belong to daughter copy B at time t. Using the probabilistic notation, the generating functions F(t,z), FA(t,z), and FB(t,z) can be represented by
where 1A denotes the indicator function of the set A.
We now make a crucial observation that conditioned on α(wT)=j, i.e. the gene state of the mother copy is j at replication, the dynamics {αA(t),XA(t)}wT≤t≤T for daughter copy A and the dynamics {αB(t),XB(t)}wT≤t≤T for daughter copy B are independent of each other. This shows that
where we have used the fact that Fk(wT,0)=P(α(wT)=k). Note that E[(z+1)XA(t)|α(wT)=k] is the generating function of pnA(t) conditioned on α(wT)=k. Similarly to the proof of Equation 56, it is easy to prove that
where Lj′, i,j=0,1 are functions obtained from Lj by substituting the parameters r, a, b, and u with the parameters r′, a′, b′, and u′, respectively,
and
Note that E[(z+1)XB(t)|α(wT)=k] is the generating function of pnB(t) conditioned on α(wT)=k. Similarly to the proof of Equation 72, we have
where
and
Inserting Equations 72 and 73 into Equation 71, we finally obtain the generating function F after replication, i.e.
In summary, we have derived the analytical expression of the generating function F at any time t∈[0,T] within a cell cycle, which is given by
where Fi(wT,z), i=0,1 are determined by Equations 66 and 67. The time-dependent distribution of the mRNA number can be recovered by taking the derivatives of the generating function F at z=−1.
We next compute the generating functions Fi, i=0,1 after replication. Recall that Fi(t,z) is the generating function of
which denotes the probability of having n mRNA molecules in the cell at time t when the daughter copy A is in state i. In addition, recall that FiA(t,z) is the generating function of
which denotes the probability of having n mRNA molecules that belong to daughter copy A at time t when daughter copy A is in state i. It is easy to see that the generating functions Fi(t,z) and FiA(t,z) can be represented by
where 1A denotes the indicator function of the set A.
We now make a crucial observation that conditioned on α(wT)=j, i.e. the gene state of the mother copy is j at replication, the dynamics {αA(t),XA(t)}wT≤t≤T for daughter copy A and the dynamics {αB(t),XB(t)}wT≤t≤T for daughter copy B are independent of each other. This shows that
where we have used the fact that Fk(wT,0)=P(α(wT)=k). Note that E[(z+1)XA(t)1{αA(t)=i}|α(wT)=k] is the generating function of pi,nA(t) conditioned on α(wT)=k. Similarly to the proof of Equations 66 and 67, it is easy to prove that
where Kij′, i,j=0,1 are functions obtained from Kij by substituting the parameters r, a, b, u, and v with the parameters r′, a′, b′, u′, and v′, respectively,
and
Inserting Equations 73 and 76 into Equation 75, we obtain
Next we focus on two non-trivial special cases where the time-dependent distribution given in Equation 11 can be greatly simplified. We assume the same setup as Figure 2, i.e. initially there is no mRNA molecules in the cell and the gene is in the inactive state. In this case, we have F0(0,z)=1 and F1(0,z)=0.
The first special case occurs when the gene switches rapidly between the two states, i.e. σ0,σ1≫d,g. In this case, we have a,b,r≫1 and thus the confluent hypergeometric function term in Equation 9 reduces to
Direct computations show that the generating function given in Equation 10 reduces to F(t,z)=eμ(t)z, where
with λ=au/b and λ′=a′u′/b′. Then the time-dependent mRNA distribution is given by
which is a Poisson distribution with mean μ(t).
The second special case occurs when the mRNA is produced in a bursty manner, i.e. σ0≫σ1 and ρ/σ0 is finite.^74^ In this case, the gene is mostly in the inactive state, but when it becomes active, it produces a large number of mRNA molecules. Clearly, we have b≫a, u/b is finite, and thus the confluent hypergeometric function terms in Equation 9 reduce to
where we have used Kummer’s transformation in the first equation. Straightforward computations show that the generating function given in Equation 10 can be simplified as
where A(t)=(u/b)e−dt is a decay term due to mRNA degradation and B(t)=(u/b)eβgt is the mean burst size at time t∈[0,wT]. In the bursty limit, the burst frequency for each gene copy decreases from σ1 to σ1′ upon replication. When σ1′=σ1/2, the total burst frequency does not change when replication occurs and dosage compensation is perfect. In this case, we have 2a′=a and the generating function F reduces to
Then the time-dependent mRNA distribution is given by
If the mRNA has a much shorter lifetime compared to the cell cycle duration, i.e. d≫g, then A(t)≈0 and thus the mRNA number has the negative binomial distribution
In some special cases, the integrals in Equations 31 and 32 can be computed explicitly. The first case occurs when mRNA is unstable and the gene switches rapidly between the two states. In this case, we have a,b≫1 and thus M(a;b;z)≈eaz/b. It then follows from Equations 26 and 27 that the generating function Fss(t,z) at any time within a cell cycle is given by
where λ=au/b and λ′=a′u′/b′. Hence the steady-state distribution for lineage measurements can be recovered from the generating function Flin(z) given in Equation 31 by taking the derivatives at z=−1, i.e.
where Γ(n,λ)=∫λ∞tn−1e−tdt is the incomplete gamma function. Similarly, the steady-state distribution for population measurements can be recovered from the generating function Fpop(z) given in Equation 32 and is given by
The second case occurs when mRNA is unstable and the gene switches slowly between the two states. In this case, we have a,b≪1 and thus
From Equations 26 and 27, the generating function Fss(t,z) at any time within a cell cycle is given by
Hence it follows from Equation 31 that the lineage distribution is given by
where δ0(n) is the Kronecker delta which takes the value of 1 when n=0 and takes the value of 0 otherwise, and it follows from Equation 32 that the population distribution is given by
Recall that the steady-state means of the mRNA number for lineage and population measurements are given by
If gene replication is not taken into account, i.e. w=1, we have shown that the steady-state mRNA distribution at any time within a cell cycle is given by
Inserting Equation 88 into Equation 77 gives the explicit expressions of the lineage and population means, i.e.
The explicit expression in the general case is too complicated and is omitted here.
We have seen that the mRNA means for the two types of measurements have different expressions. A natural question is how far the lineage statistics deviates from the population one. Figures S2A and S2B show the the ratio of the lineage mean to the population mean, R1=⟨n⟩lin/⟨n⟩pop, as functions of β, η, w, and σ1′/σ1. Clearly, R1 is always greater than 1, which means that the lineage mean is greater than the population mean.^46^ As well, R1 is largest when the mRNA synthesis rate scales with cell volume (β=1), mRNA is unstable (η≫1), and there is no dosage compensation (σ1′=σ1). When these three conditions are satisfied, it follows from Equation 35 that
In this case, the lineage and population means can be obtained exactly as
and it is easy to see that R1 attains its maximum of R1max≈1.1 when w≈0.6 (Figure S2B). In other words, the lineage mean can differ from the population mean by at most 10%.
Moreover, we have also compared the variances of the mRNA number for lineage and population measurements and the ratio R2 of the lineage variance to the population variance is shown in Figures S2C and S2D as functions of β, η, w, and σ1′/σ1. Similarly, R2 is also large when the mRNA synthesis rate scales with cell volume, mRNA is unstable, and there is no dosage compensation.
Here we give the proof of Equation 38. For simplicity, we only focus on population measurements; the proof for lineage measurements is totally the same. Note that the Fano factor of the mRNA number in a population of cells can be decomposed as
When mRNA synthesis is balanced (β=1), it follows from Equations 35 and 37 that the mean mRNA number scales with the birth volume Vb and the second factor moment scales with Vb2, i.e.
This is because the parameters u, u′, λ, λ′, μ, and μ are all proportional to Vb when β=1. On the other hand, the cell volume distribution for population measurements is given by Equation 41. It is easy to see that the mean and variance of cell volume in a population of cells are given by
Hence the Fano factor of cell volume is given by
which also scales with Vb. Inserting Equations 81 and 82 into Equation 80, we immediately obtain Equation 38.
Another important question is whether the conventional telegraph model without a cell cycle description can capture the dynamic properties of the detailed telegraph model with such a description. Note that this is impossible when gene replication is taken into account since the mRNA distribution for the former has at most two modes, while the latter can have more than two modes. Hence, in the following, we only focus on the case when gene replication is not taken into account, i.e. w=1.^49^
Recall that the conventional telegraph model (with no volume-dependent rates) is characterized by the effective reactions
where deff=d+g is the effective decay rate of mRNA. To make a fair comparison of the two models, the mRNA synthesis rate ρ¯ of the conventional model is taken to be the time-average of the detailed model within a cell cycle, i.e.
It is well known^12^ that the steady-state mRNA mean for the conventional model is given by
where λ=au/b and η=d/g. Interestingly, when the mRNA synthesis rate is proportional to cell volume (β=1), we have ⟨n⟩conv=λ/(log2), which agrees with the lineage mean given in Equation 79. When the mRNA synthesis rate is independent of cell volume (β=0), we have ⟨n⟩conv=ηλ/(η+1), which coincides with the population mean given in (79). In addition, it follows from Equations 79 and 83 that the means for the conventional and detailed models are related by
This shows that the mean ratio R′=⟨n⟩conv/⟨n⟩pop for the two models only depends on β (Figure S3A). It attains its minimum Rmin′=1 when β=0 and attains its maximum Rmax′=1/2(log2)2≈1.04 when β=1. As a result, the conventional model can accurately capture the mRNA mean of the detailed model.
While the conventional model can capture the mean of the detailed model, it cannot accurately capture the mRNA distribution. To see this, recall that the steady-state distribution pnconv of the mRNA number for the conventional telegraph model has the generating function^12^
where a¯=σ1/deff, b¯=(σ0+σ1)/deff, and u¯=ρ¯/deff. Figure S3B illustrates the Hellinger distance between the distributions for the two models as a function of β and η. It can be seen that they coincide with each other when the transcription rate is volume-independent (β=0) and mRNA is unstable (η≫1). In other cases, they deviate from each other significantly — the conventional model has a much smaller gene expression noise than the detailed model (Figure S3C). This can be explained as follows. When w=1, the steady-state inactive probability of the gene at birth is poffb=(b−a)/b and the active probability of the gene at birth is ponb=a/b. Hence for unstable mRNAs, the time-dependent generating function can be simplified significantly as
In particular, when β=0 and η≫1, the steady-state mRNA distribution at any time within a cell cycle is independent of time t. Hence the generating functions of the lineage and population distributions are the same and are given by
where a=σ1/d, b=(σ0+σ1)/d, and u=ρ/d. When β=0 and η≫1, we have ρ¯=ρ and d/deff≈1. In this case, the two generating functions given in Equations 84 and 86 are approximately equal and thus the conventional model makes the correct prediction. We emphasize that there have been numerous studies that estimated the rate constants of stochastic gene expression dynamics based on the conventional telegraph model.^9^^,^^19^ Our results suggest that parameters estimated using this approach maybe unreliable.
While the conventional model fails to capture the mRNA distribution of the detailed model, we find that that it is capable of capturing the modality (unimodality or bimodality) of the distribution. To see this, following,^91^^,^^92^ we define the strength of bimodality as
where Hlow is the height of the lower peak, Hhigh is the height of the higher peak, and Hvalley is the height of the valley between them. For unimodal distributions, κ is automatically set to be 0. For bimodal distributions, κ is a quantity between 0 and 1 since Hvalley<Hlow≤Hhigh. In general, to display strong bimodality, the following two conditions are (i) the two peaks should have similar heights and (ii) there should be a deep valley between them. The former ensures that the time periods spent in the low and high expression states are comparable, while the latter guarantees that the two expression levels are distinguishable. Clearly, κ is large if both conditions are satisfied and is small if any one of the two conditions is violated. Hence, κ serves as an effective indicator that characterizes the strength of bimodality.
Figures S3D and S3E illustrate κ as a function of the gene switching rates σ0 and σ1 for the two models. Clearly, both models display unimodality in the regime of fast gene switching and display bimodality in the slow switching regime. Furthermore, we find that the regions in parameter space where the two models show bimodality are very close to each other, except that the detailed model needs a larger gene activation rate σ1 to obtain the same strength of bimodality (shown by the blue and orange dashed lines in Figures S3D and S3E). This indicates the two models in general show the same modality but the heights of the modes may be different.
Here we focus on two special cases for mRNA to display concentration homeostasis. The first special case takes place where gene replication is not taken into account (w=1), the gene activation and inactivation probabilities at birth are given by poffb=(b−a)/b and ponb=a/b. This implies that μ=0 and thus the steady-state mean at birth can be simplified to a large extent as
Inserting this equation into Equation 35 yields
In particular, when β=1, we have ⟨n⟩ss(0)=λ and ⟨n⟩ss(t)=λegt=(λ/Vb)V(t). In this case, the mRNA displays concentration homeostasis, i.e. constant mean concentration across the cell cycle.
Another special case occurs when the mRNA is produced in a bursty manner (σ0≫σ1) and when dosage compensation is perfect (σ1′=σ1/2). In this case the gene is mostly off and thus poffb≈poffr≈1. This further implies that λ′=2βw−1λ and μ=μ′=0. Inserting these equations into Equation 36 shows that the steady-state mean at birth is still given by Equation 87 and the time-dependent mean under cyclo-stationary conditions is still given by Equation 88. In particular, when β=1, we have ⟨n⟩ss(0)=λ and ⟨n⟩ss(t)=λegt=(λ/Vb)V(t). Hence when mRNA synthesis is balanced and bursty and when dosage compensation is perfect, the mRNA also displays concentration homeostasis.
To determine the value of σϵ in real systems, we examined the lineage data collected in E. coli and haploid fission yeast cells using a mother machine.^36^^,^^84^ The E. coli data set contains the lineage measurements of cell length at three different temperatures (25o C, 27o C, and 37o C). The fission yeast data set contains the lineage measurements of cell area under seven different growth conditions (Edinburgh minimal medium (EMM) at 28o C, 30o C, 32o C, and 34o C and yeast extract medium (YE) at 28o C, 30o C, and 34o C).
The inference of σϵ can be divided into the following three steps. First, the mean of the birth volume Vb across all generations gives an estimate of v¯. Next, the slope of the regression line of the division volume Vd on the birth volume Vb gives an estimate of α. Finally, since α and v¯ have been determined, σϵ can be estimated as the sample standard deviation of ϵ=Vd−αVb−(2−α)v¯. The estimated values of σϵ in E. coli and fission yeast under all growth conditions are listed in Table S2.
Given the values of the four variables αb, Nb, Vb, and T, the generating function F at any time t∈[0,T] within a cell cycle is given by Equation 44. To proceed, let Π(k)(i,n,x,τ)=P(k)(αb=i,Nb=n,Vb=x,T=τ) denote the joint distribution of the four variables in generation k. Then the generating function F at any fixed proportion θ∈[0,1] of the cell cycle in that generation is given by
We next focus on the mRNA distribution in generation k+1. Similarly, once the values of the four variables are fixed, the generating functions Fi, i=0,1 at division are given by Equation 14, i.e.
where the initial conditions Fi(0,z), i=0,1 are determined by Equation 45. To proceed, let αd denote the state of daughter copy A at division, let Nd denote the mRNA number at division, and let Vd denote the cell volume at division. From Equation 90, we know the conditional joint distribution of αd and Nd, i.e.
Hence the joint distribution of αd, Nd, Vb, and T in generation k is given by
From this it follows that the joint distribution of αd, Nd, and Vd in generation k is given by
Since we have assumed binomial partitioning of molecules at division, the joint distribution of αb, Nb, and Vb in generation k+1 is given by
Hence the joint distribution of the four variables αb, Nb, Vb, and T in generation k+1 is given by
where the conditional distribution of T given Vb can be computed from Equation 43 as
Applying Equations 89 and 91 repeatedly, we can compute the exact mRNA distribution at any time within a cell cycle in all generations. Finally, the steady-state joint distribution of the four variables αb, Nb, Vb, and T is given by
Inserting this equation into Equation 89 gives the transient mRNA distribution under cyclo-stationary conditions.
Once the joint distribution of the four variables is known, it follows from Equation 89 that the steady-state mRNA distribution for lineage measurements is given by
where ⟨T⟩≈log(2)/g is the mean doubling time. The exact expression for the steady-state population distribution is difficult to obtain. However, when the variability is cell cycle duration is small (σϵ≪1), an approximation of the steady-state population distribution is given by
Note that when the cell volume dynamics is stochastic, the analytical mRNA distributions involve multiple integration, which is difficult to compute using the conventional multi-grid method. An alternative strategy to compute the multiple integration is to use the Monte Carlo method with the joint distribution Πss(i,n,x,τ) being randomly sampled.
Recall that before replication, there is only one gene copy in the cell and the EDM is given by Equation 29. After replication, there are two gene copies in the cell and the EDM is given by Equation 30. At cell birth, there is only one gene copy and thus the ENM approximation of the mRNA distribution is given by
where pEDM(n|V) is the steady-state distribution of the EDM given in Equation 29 and
is the cell volume distribution at birth which is Gaussian.^83^ Similarly, we can construct the ENM approximations for the mRNA distributions at replication and division.
We next construct the ENM approximation for the mRNA distribution of lineage measurements. This is much more complicated because for a given cell of volume V, it is unclear whether it has one or two gene copies. To determine the number of gene copies in the cell, we also need the information of the birth volume Vb and the cell cycle duration T. Once the values of V, Vb, and T are known, the age of the cell is t=(1/g)log(V/Vb). It has only one gene copy when t<wT and has two gene copies when ≥wT . Hence the ENM approximation for the lineage distribution is given by
Here pEDM(n|t,x,τ) is steady-state distribution of the EDM for a cell of age t given that the birth volume is x and the cell cycle duration is τ, and
is the joint distribution of Vb and T. Note that the reaction scheme given in Equation 29 should be used for the EDM when t<wτ and the reaction scheme given in Equation 30 should be used when t≥wτ.
We are grateful to Prof. Qiwen Sun and Prof. Feng Jiao for pointing us to useful references. C.J. acknowledges support from National Natural Science Foundation of China with grant No. U1930402, No. 12271020, and No. 12131005. R.G. acknowledges support from the Leverhulme Trust (RPG-2020-327).
R.G. conceived the original idea. C.J. performed the theoretical derivations and numerical simulations. C.J. and R.G. interpreted the theoretical results and jointly wrote the manuscript.
The authors declare that they have no competing interests.
Published: January 20, 2023