Quantitative Understanding of Probabilistic Behavior of Living Cells Operated by Vibrant Intracellular Networks

For quantitative understanding of probabilistic behaviors of living cells, it is essential to construct a correct mathematical description of intracellular networks interacting with complex cell environments, which has been a formidable task. Here, we present a novel model and stochastic kinetics for an intracellular network interacting with hidden cell environments, employing a complete description of cell state dynamics and its coupling to the system network. Our analysis reveals that various environmental effects on the product number fluctuation of intracellular reaction networks can be collectively characterized by Laplace transform of the time-correlation function of the product creation rate fluctuation with the Laplace variable being the product decay rate. On the basis of the latter result, we propose an efficient method for quantitative analysis of the chemical fluctuation produced by intracellular networks coupled to hidden cell environments. By applying the present approach to the gene expression network, we obtain simple analytic results for the gene expression variability and the environment-induced correlations between the expression levels of mutually noninteracting genes. The theoretical results compose a unified framework for quantitative understanding of various gene expression statistics observed across a number of different systems with a small number of adjustable parameters with clear physical meanings.

Popular Summary

Live cells with the same DNA produce different numbers of protein molecules, which results in variations in the biological behaviors of a clonal population of cells. Quantitatively understanding intracellular chemical fluctuations and their impact on the probabilistic behaviors of living cells is one of the most challenging goals of modern biophysical science. To achieve the goal, it is necessary to construct a correct mathematical model of intracellular networks interacting with cell environment, which has remained a formidable task because cell environments and their coupling to the system network are too complex to be fully represented by conventional models. Here, we present a novel model and stochastic kinetics optimized for intracellular networks; we treat the cell state dynamics and their coupling to the system network in an implicit and exact manner while modeling the system network explicitly.

Taking a new theoretical approach, we discover a general principle governing the chemical fluctuations produced by intracellular networks. On the basis of our results, we propose an efficient method for quantitatively analyzing intracellular chemical fluctuations. By applying our approach to the gene expression network, we obtain simple analytic results for the gene expression level variation among a clonal population of cells. Our results provide excellent quantitative explanations of the various gene expression statistics observed across a number of different gene expression systems in E. coli and S. cerevisiae in a unified manner.

Our results present a new paradigm for a quantitative understanding of intracellular chemical fluctuations.

Figure 1

Concept of vibrant reaction process coupled to complex and hidden cell environment. Schematic representations of (a) an elementary reaction process under a homogeneous environment and (b) a vibrant reaction process with the rate coefficient being a stochastic variable that differs from cell to cell and fluctuates over time (Fig. 1 of Supplemental Material [32]). (c) Coupling of a reaction system to a complex and hidden cell environment makes the rate-coefficient stochastic variable dependent on the cell environmental state, the correct, explicit model of which cannot be easily constructed for lack of available information. A reaction process whose rate coefficient is coupled to hidden environment will be designated by vibrant reaction process and represented by the wavy arrow.

Figure 2

Effects of environment-coupled fluctuations in the rate coefficient on the fluctuation in the product number. (a) Schematic representation of the vibrant reaction and ensuing decay of the product molecule with constant rate $\gamma$. The reaction rate $R$ is factored into the reactant population-dependent rate factor $\theta$ and the rate coefficient $k$, which are dependent on environmental variable $r$. ${\eta }_z^2(\equiv ⟨\delta z^2⟩/⟨z{⟩}^2)$ is given by Eq. (2). (b),(c) Comparison between the prediction of Eq. (2) for Fano factor $F_z(={\eta }_z^2⟨z⟩)$ and stochastic simulation results on two different models of environment-coupled reaction. We focus on the situation where the fluctuation in the reaction rate $R$ is dominated by the fluctuation in the rate coefficient $k$ by setting $\theta$ equal to unity. $F_z-1$ is proportional to $F_k(=⟨\delta k^2⟩/⟨k⟩)$ with the proportionality coefficient given by ${\widehat{\varphi }}_k(\gamma )={\int }_0^{\infty }dt{e}^{-\gamma t}{\varphi }_k(t)$ with ${\varphi }_k(t)$ being the normalized time-correlation function of rate-coefficient fluctuation [Eq. (c1)]. Examples of stochastic time trace of $k$, which is a function of stochastic variable $r(t)$ (left-hand panel), ${\varphi }_k(t)$ (middle panel), and $F_z(\gamma )-1$ (right-hand panel), are displayed for the following models. (b) $k(r)=k_0+\kappa r^2$ with stochastic dynamics of $r$ being the stationary Gaussian Markov process known as Ornstein-Uhlenbeck process [46]. $\lambda$ stands for the relaxation rate of the fluctuation in $r$ defined by $⟨r(t)r(0)⟩=⟨r^2⟩e^{-\lambda t/2}$. $F_z$ is a monotonically decaying function of $\gamma$. (c) $k(r)=\kappa \exp [-\beta r]$ with stochastic dynamics of $r$ being the stationary Gaussian non-Markov process describing the position of the weakly damped harmonic oscillator with mass $m$ under the frictional force $2m\lambda dr(t)/dt$ [47]. For this model, $⟨r(t)r(0)⟩$ is given by $⟨r^2⟩e^{-\lambda t}[\cos (\omega t)+(\lambda /\omega )\sin (\omega t)]$ with $\omega =\sqrt{{\omega }_0^2-{\lambda }^2}$, where ${\omega }_0$ is the natural frequency of the oscillator. When ${\varphi }_k(t)$ is an oscillatory function, $F_z-1$ can be a nonmonotonic function of $\gamma$. See Appendix pp3 for the details of the models.

Figure 3

Analyses of the gene expression variability in the synthetic E. coli with T7 RNAP driven transcription. (a) Our gene expression network model. RNAP bound fraction $\theta$ of promoter is explicitly modeled as a function of the number $N_{Rp}$ of RNAP, distribution of which is under control in our system. For the model, the simple analytic result of the noise ${\eta }_{m(p)}^2$ in mRNA (protein) level can be obtained, which has the same structure as Eq. (2b). (b),(c) The best global fit of Eq. (6) to the experimental data for the dependence of the protein noise on the mean RNAP level and RNAP noise. See Eq. (i1) for the explicit expression of Eq. (6) in terms of the control variables in the experiment. (d),(e) Dependence of protein noise on the mean RNAP occupancy of promoter. The protein noise comprises the intrinsic noise, $(1+\alpha )/⟨p⟩$, and the extrinsic noise, ${\eta }_{p,\theta }^2(={\beta }_{p\theta }{\eta }_{\theta }^2)$ and ${\eta }_{p,k}^2(={\beta }_{p{k}_{TX}}{\eta }_{k_{TX}}^2+{\chi }_{p{k}_{TL}}{\eta }_{k_{TL}}^2)$, propagated from environment-coupled fluctuations in $\theta$ and from that in $k_{TX}$ and $k_{TL}$. The RNAP noise among the sampled cells is controlled to be $0.1\pm 0.002$. The error bar or the standard deviation of each experimental result is obtained over 10 different samples, each of which contains about 1500 cells. The size of the error bar is smaller than or comparable to the symbol size (see Fig. 6 of Supplemental Material [32]). Details of the data analysis are presented in Appendix pp9.

Figure 4

(a),(b) The best global fit of Eq. (7) to the experimental data for the dependence of the mean protein level on the mean RNAP level and RNAP noise. See Eq. (i2) for the explicit expression of Eq. (7) in terms of the control variables in the experiment. (c) Dependence of the mean protein level scaled by its maximum or the mean RNAP occupancy of promoter, $⟨p⟩/{⟨p⟩}_{\max }[=⟨\theta (\zeta )⟩]$, on $\theta (\overline{\zeta })$. Inset: Dependence of $⟨\theta (\zeta )⟩$ on $\overline{\zeta }$ (see Fig. 9 of Supplemental Material [32]). $⟨\theta ⟩=\theta (\overline{\zeta })$ for ${\eta }_{\zeta }^2=0$ (dashed line). ${\eta }_{\zeta }^2=0.44$ and 0.59 for CFP and mCherry (Table 1). (d) Dependence of ${\eta }_{p,\theta }^2$ scaled by its maximum ${\eta }_{p,\theta }^{2\ast }$ on $\theta (\overline{\zeta })$. For both CFP and mCherry, ${\eta }_{p,\theta }^2/{\eta }_{p,\theta }^{2*}\cong {[1-\theta (\overline{\zeta })]}^2$. Data for ${\eta }_{p,\theta }^2$ are obtained from experimental data for ${\eta }_p^2$ by ${\eta }_{p,\theta }^2={\eta }_p^2-[{({\eta }_p^{\mathrm{int}})}^2+{\eta }_{p,k}^2]$, with ${({\eta }_p^{\mathrm{int}})}^2$ given by Eq. (f1) and ${\eta }_{p,k}^2$ given in Table 1. The value of the RNAP noise ${\eta }_{N_{Rp}}^2$ in (c) and (d) is the same as in Figs. 3 and 3.

Figure 5

Unified quantitative explanation of various patterns in the Fano factor data observed across different gene expression systems. Two parameters, $q(\equiv {\eta }_{p,k}^2/{\eta }_{p,\theta }^{2*})$ and ${\eta }_{\zeta }^2$, determine the pattern in the dependence of Fano factor $F_p$ on $⟨\theta ⟩$ [Eq. (g5)]. (a) When $q<1/3$, $F_p(⟨\theta ⟩)$ is a nonmonotonic function with the single maximum at $⟨\theta ⟩={⟨\theta ⟩}^{\ast }(<2/3)$. (b) ${⟨\theta ⟩}^{\ast }$ decreases with ${\eta }_{\zeta }^2$. (c) When $q>1/3$, $F_p(⟨\theta ⟩)$ monotonically increases with $⟨\theta ⟩$. (d),(e) Dependence of the intensity Fano factor $F_I(\propto F_p)$ on $⟨\theta ⟩$ for CFP and mCherry in our synthetic E. coli with T7 RNAP transcription (d) and for yEGFP gene expression in S. cerevisiae under two different transcriptional controls, reported in Ref. [24] (e). In (d), the theory curves are the results of Eqs. (6) and (7) with the predetermined parameters values extracted from the analysis in Fig. 3 (Table 1). Equations (6) and (7) also provide a unified quantitative explanation of all the Fano factor data shown in (e) with a small number of adjustable parameters (Table I of Supplemental Material [32]). The extracted value of $q$ is far smaller than $1/3$ for data in (d) and (e), where the Fano factor has the maximum at about ${⟨\theta ⟩}^{\ast }\cong 3^{-1}(2-\sqrt{1-3q})\cong 1/3$. The small value of parameter $q$ means that the extrinsic noise originating from RNAP-promoter association step makes the dominant contribution to the total extrinsic protein noise. In the opposite case, $q$ becomes very large and the Fano factor becomes a linear function of $⟨\theta ⟩$, which is the case for the abundant proteins in the wild-type E. coli in (f) [see also Figs. 13(b) and 14 in Supplemental Material 32]. See text for the analytic expression of $F_p(⟨\theta ⟩)$.

Figure 6

Quantitative analysis of the environment-induced correlation between the expression levels of dual reporter genes. (a) Our gene expression network model for the dual reporter system expressing gene $A$ (CFP) and gene $B$ (mCherry). (b) The best global fit of Eq. (10) with Eq. (11) to the experimental data for the dependence of the mean scaled correlation on the mean RNAP level and RNAP noise. See Eq. (i3) for the expression of Eq. (10) in terms of the control variables in the experiment. (c)–(e) The environment-induced correlation between the CFP $(p_A)$ and mCherry $(p_B)$ levels measured by $C_p[\equiv ⟨\delta p_A\delta p_B⟩/(⟨p_A⟩⟨p_B⟩)]$, ${\varphi }_p(\equiv ⟨\delta p_A\delta p_B⟩/\sqrt{⟨\delta p_A^2⟩⟨\delta p_B^2⟩})$, and $⟨\delta p_A\delta p_B⟩$. They have two components: one originating from the correlated fluctuations of ${\theta }_A$ and ${\theta }_B$ (red surface) and the other from the correlated fluctuation of $k_{TX(TL)}^{(A)}$ and $k_{TX(TL)}^{(B)}$ (blue surface). (f)–(h) Dependence of $C_p$, ${\varphi }_p$, and $⟨\delta p_A\delta p_B⟩$ on $\theta ({\overline{\zeta }}_A)$. In our system, $\theta ({\overline{\zeta }}_B)$ is dependent on $\theta ({\overline{\zeta }}_A)$, as shown in the inset of (f). They are not the same because the value of ${\overline{K}}_A$ is different from that of ${\overline{K}}_B$ (Table 1). The red (blue) line represents the contribution of the red (blue) surface in (c)–(e). In (f), the theory curve represents the best fit of Eq. (10) to experimental data. In (g) and (h), the theory curves represent the prediction of Eqs. (6), (7), and (10) with predetermined parameter values in Table 1. (i) $C_p$ linearly increases with the probability that both promoters of CFP and mCherry are unoccupied by RNAP. (j) The joint distribution of CFP and mCherry levels represented by heat map (experiment) and the contour plot (joint gamma distribution). See Fig. 16 of Supplemental Material [32] and Appendix pp11 for more details. Near the most probable protein level, $(p_A^{\ast },p_B^{\ast })$, the joint distribution becomes the two-dimensional Gaussian (Note 8 in and Fig. 17 in Supplemental Material [32]) with standard deviation ${\sigma }_{+}$ (${\sigma }_{-}$) along the positive (negative) diagonal direction. (k) Dependence of ${\sigma }_{+}/{\sigma }_{-}$, ${\eta }_{\text{dual}}^{\text{ext}}/{\eta }_{\text{dual}}^{\mathrm{int}}$, and ${\overline{\eta }}_p^{\text{ext}}/{\overline{\eta }}_p^{\text{int}}$ on $\theta ({\overline{\zeta }}_A)$. ${\sigma }_{+}/{\sigma }_{-}>{\overline{\eta }}_p^{\text{ext}}/{\overline{\eta }}_p^{\text{int}}\ge {\eta }_{\text{dual}}^{\text{ext}}/{\eta }_{\text{dual}}^{\mathrm{int}}$. The value of the RNAP noise ${\eta }_{N_{Rp}}^2$ is the same as in Figs. 3 and 3.