Shaping protein distributions in stochastic self-regulated gene expression networks

In this work, we study connections between dynamic behavior and network parameters, for self-regulatory networks. To that aim, a method to compute the regions in the space of parameters that sustain bimodal or binary protein distributions has been developed. Such regions are indicative of stochastic dynamics manifested either as transitions between absence and presence of protein or between two positive protein levels. The method is based on the continuous approximation of the chemical master equation, unlike other approaches that make use of a deterministic description, which as will be shown can be misleading. We find that bimodal behavior is a ubiquitous phenomenon in cooperative gene expression networks under positive feedback. It appears for any range of transcription and translation rate constants whenever leakage remains below a critical threshold. Above such a threshold, the region in the parameters space which sustains bimodality persists, although restricted to low transcription and high translation rate constants. Remarkably, such a threshold is independent of the transcription or translation rates or the proportion of an active or inactive promoter and depends only on the level of cooperativity. The proposed method can be employed to identify bimodal or binary distributions leading to stochastic dynamics with specific switching properties, by searching inside the parameter regions that sustain such behavior.

Figure 1

Schematic representation of the transcription-translation mechanism under study. The promoter associated to the gene of interest is assumed to switch between active $\left({\mathrm{DNA}}_{\mathrm{on}}\right)$ and inactive $\left({\mathrm{DNA}}_{\mathrm{off}}\right)$ states, with rate constants $k_{\mathrm{on}}$ and $k_{\mathrm{off}}$ per unit time, respectively. In this study, the transition is assumed to be controlled by a feedback mechanism induced by the binding and unbinding of a given number of $X$-protein molecules, which makes the network self-regulated. Transcription of messenger RNA $\left(\text{mRNA}\right)$ from the active $\text{DNA}$ form and translation into protein $X$ are assumed to occur at rates (per unit time) $k_1$ and $k_2$, respectively. $k_{\varepsilon }$ is the rate constant associated to transcriptional leakage. Both $\text{mRNA}$ and $X$-protein degradation are assumed to occur by first-order processes with rate constants ${\gamma }_1$ and ${\gamma }_2$, respectively.

Figure 2

Transient distributions for positive feedback $\left(H=-4\right)$ with leakage $\varepsilon =0.2$. In order to compare accuracy, distributions were computed from CME (green lines), Eq. (4) (black lines), and Gillespie SSA realizations (on the order of $10000$) to reconstruct the distributions (blue lines). Parameters employed to simulate the self-regulation mechanism are presented in Table 1. The initial condition is given by a delta function centered at $x=0$ to represent the absence of protein at $\tau =0$. The snapshot at the right of the plot corresponds to a stationary distribution.

Figure 3

The qualitative shape of all possible distributions for $H=-1$ depending on the intersections between $\rho \left(x\right)$ and $r\left(x\right)$ in Eq. (8) (black and green lines, respectively, on the first column plots). The different $r\left(x\right)$ are marked with a letter to identify which shape corresponds to each line. For $r\left(0\right)\ge 1$ there is at most one positive intersection that corresponds with graded responses (rows a and b). For $0<r\left(0\right)<1$ there can be two positive intersections leading to binary distributions (row c). For $r\left(0\right)\le 0$ no intersection occurs (row d).

Figure 4

The qualitative shape of all possible distributions for $H<-1$ depending on the intersections between $\rho \left(x\right)$ and $r\left(x\right)$ in Eq. (8) (black and green lines, respectively, on the first file plots). The different $r\left(x\right)$ are marked with a letter to identify which shape corresponds to each line. For $r\left(0\right)>1$ (row a), the number of intersection points ranges from one to three as the slope decreases, which corresponds with graded and bimodal distributions. For $0<r\left(0\right)\le 1$ there can be at most two positive intersections that correspond with binary distributions (row b and c). For $r\left(0\right)\le 0$ no intersection occurs (row d).

Figure 5

The upper plot depicts function $\rho \left(x\right)$ as given in (2) for $H<-1$ and two possible functions $r_1\left(x\right)$ and $r_2\left(x\right)$ tangent to $\rho \left(x\right)$ at points $({\overline{x}}_1$ and ${\overline{x}}_2)$, respectively. The lower plot represents the first derivative ${\rho }^{\prime }\left(x\right)$ and a possible value of $r^{\prime }\left(x\right)$ represented by the horizontal dashed line. Values $x\in \mathcal{S}$ which happen to be within the interval $({\overline{x}}_1,{\overline{x}}_2)$ will satisfy condition (12) for a minimum since ${\rho }^{\prime }\left(x\right)<r^{\prime }\left(x\right)$. On the other hand, maxima will lie in intervals $(0,{\overline{x}}_1)$ and $({\overline{x}}_2,\infty )$.

Figure 6

Possible shapes of the binary/bimodal regions for $K=70$ and $H=-4$ [according to expression (20), ${\varepsilon }^{*}=0.36$ for these parameter values]. (a) Striplike region corresponding to $\varepsilon <{\varepsilon }^{*}$. (b) Hornlike region corresponding to $\varepsilon >{\varepsilon }^{*}$. Contours in plot (a) provide an indication of the relative error between the analytical distribution and the distribution computed from the CME. For most of the region the error remains quite low, rising to about $4\%$ only near the lower bound of the region. Plots (c) and (d) represent the stationary distributions associated to the points marked in the corresponding regions (a and b, respectively). For comparison purposes, distributions are computed from the analytical solution, Eq. (6), as well as from the CME. Plots (e) and (f) depict the corresponding stochastic dynamics (SSA). Plot (e) reflects a scenario with rare transitions between high and low protein levels, while plot (f) represents a noisy behavior with frequent transitions between high and low protein levels.

Figure 7

Intersections of $\rho \left(x\right)$ (continuous line) with $r\left(x\right)$ [Eq. (9)] associated to stochastic analysis (dashed line), and $r_d\left(x\right)$ [Eq. (24)] associated to deterministic ODEs (dotted line), for a self-regulatory network with parameters $H=-4,K=70,\varepsilon =0.2,a=9,b=18$. $r_d\left(x\right)$ is shifted with respect to $r\left(x\right)$, which leads to one intersection point instead of three. Consequently the deterministic system shows one single (monostable) equilibrium state, whereas the stochastic counterpart has a bimodal distribution compatible with a transitional dynamics, switching between two positive states.

Figure 8

A self-regulatory network with parameters $H=-4,K=70,\varepsilon =0.2,a=9,b=18$. (a) Stationary distribution computed from the CME (dashed line) and the analytical distribution (continuous line). The two most frequent states in the distribution correspond with the intersections between $\rho \left(x\right)$ and $r\left(x\right)$ shown in Fig. 7. (b) The corresponding stochastic dynamics (SSA) with regular transitions taking place between a low and a high protein level.

Figure 9

Trajectories exhibited by the corresponding deterministic system. (a) Trajectories associated to a network with parameters $H=-4,K=70,\varepsilon =0.2,a=9,b=18$ (see also Fig. 8) evolving to a stable equilibrium point. (b) Trajectories associated to a network with parameters $H=-4,K=70,\varepsilon =0.1,a=40,b=3.7$. Red circle and black dots correspond with the unstable and stable equilibrium points, respectively.

Figure 10

Regions sustaining bimodal and binary distributions in the parameters space. Regions in the $a-b$ parameter space reduce to the strip shaded area. Areas in gray (shaded region 2) denote regions with bimodal distributions. Dashed lines represent the bistable region boundaries of the deterministic model. (a) Effect of $H$ on the region area for $K=70$, and $\varepsilon =0.2$, below the critical value ${\varepsilon }^{*}=0.25$ (for $H=-3$), ${\varepsilon }^{*}=0.36$ (for $H=-4$), and ${\varepsilon }^{*}=0.6694$ (for $H=-10$). (b) Effect of $K$ on the region area, with $K=40,K=70$, and $K=100$, for $\varepsilon =0.2$ and $H=-4$. (c) Effect of $\varepsilon$ on the region area, with $\varepsilon =0.1,\varepsilon =0.2$, and $\varepsilon =0.3$, for $K=70$ and $H=-4$.