Stochastic epigenetic dynamics of gene switching

Epigenetic modifications of histones crucially affect eukaryotic gene activity, while the epigenetic histone state is largely determined by the binding of specific factors such as the transcription factors (TFs) to DNA. Here, the way in which the TFs and the histone state are dynamically correlated is not obvious when the TF synthesis is regulated by the histone state. This type of feedback regulatory relation is ubiquitous in gene networks to determine cell fate in differentiation and other cell transformations. To gain insights into such dynamical feedback regulations, we theoretically analyze a model of epigenetic gene switching by extending the Doi-Peliti operator formalism of reaction kinetics to the problem of coupled molecular processes. Spin-1 and spin-1/2 coherent-state representations are introduced to describe stochastic reactions of histones and binding or unbinding of TFs in a unified way, which provides a concise view of the effects of timescale difference among these molecular processes; even in the case that binding or unbinding of TFs to or from DNA is adiabatically fast, the slow nonadiabatic histone dynamics gives rise to a distinct circular flow of the probability flux around basins in the landscape of the gene state distribution, which leads to hysteresis in gene switching. In contrast to the general belief that the change in the amount of TF precedes the histone state change, flux drives histones to be modified prior to the change in the amount of TF in self-regulating circuits. Flux-landscape analyses shed light on the nonlinear nonadiabatic mechanism of epigenetic cell fate decision making.

Figure 1

Schematic of a self-regulating gene. Close-up views of chromatin at around the promoter region are illustrated at the bottom, where histones (brown) are actively marked (red) or repressively marked (green). Histones in the chromatin are collectively marked through their mutual interactions. When the activating mark is dominant, chromatin has an open structure, which enhances the gene transcription activity ($s=1$), and when the repressive mark is dominant, condensed chromatin suppresses the gene activity ($s=-1$). When neither mark is dominant, the histone state is neutral ($s=0$). The transcription factor (TF) is a dimer of the product protein (blue oval). The bound TF recruits histone modifiers (black), which trigger the collective histone modification as was observed in engineered cells [40, 41, 44] by perturbing the transition rates $r_{s{s}^{\prime }}$ among the histone states. The protein production rate $g_{\sigma s}$ depends both on the TF binding status, $\sigma$, and on the histone state, $s$. See the text for the definition of rates denoted by the arrows.

Figure 2

Relaxation of a self-activating gene to the repressive state. Equation (9) was numerically calculated by applying the adiabatic approximation of TF binding or unbinding. The protein concentration $p\left(t\right)$ (black line) and the histone state $\zeta \left(t\right)$ (red line) were derived by averaging 1000 simulated trajectories. See text for the experimental estimation of $\zeta \left(t\right)$ (green dots) from the data of [40]. $k$ was set to $k=2{\mathrm{d}}^{-1}$, and the other rates were defined in units of $k$; for $t<0, r_{10}^0=r_{0-1}^0=0.7, r_{01}^0=r_{-10}^0=0.1$, and $\mathrm{\Delta }{r}_{10}=\mathrm{\Delta }{r}_{0-1}=1$; for $t\ge 0, r_{10}^0=r_{0-1}^0=r_{01}^0=0.01, r_{-10}^0=0.7$, and $\mathrm{\Delta }{r}_{10}=\mathrm{\Delta }{r}_{0-1}=0.01$. For both $t<0$ and $t\ge 0, \mathrm{\Delta }{r}_{01}=\mathrm{\Delta }{r}_{-10}=0, h_0/f=200, g_{11}=1, g_{-11}=g_{10}=0.2, g_{-10}=g_{1-1}=g_{-1-1}=0, \mathrm{\Omega }=100$, and $c=20$.

Figure 3

Landscape $U(p,\zeta )$ and probability flux $\mathbf{J}(p,\zeta )$ calculated on the 2D plane of the protein concentration $p$ and the histone state $\zeta$ in the adiabatic approximation of TF binding or unbinding. $U$ is shown with contours and $\mathbf{J}$ is shown with yellow arrows. (a–c) Self-activating and (d) self-repressing genes. (a) $h_0/f=2$, (b) $h_0/f=20$, (c) $h_0/f=200$, and (d) $h_0/f=2$. The rate parameters are scaled in units of $k$ by setting $k=1$; $g_{11}=1, g_{-11}=g_{10}=0.2, g_{-10}=g_{1-1}=g_{-1-1}=0, \mathrm{\Delta }{r}_{10}=\mathrm{\Delta }{r}_{0-1}=1, \mathrm{\Delta }{r}_{01}=\mathrm{\Delta }{r}_{-10}=0$ for (a)–(c) and $g_{-11}=1, g_{11}=g_{-10}=0.2, g_{10}=g_{1-1}=g_{-1-1}=0, \mathrm{\Delta }{r}_{10}=\mathrm{\Delta }{r}_{0-1}=0$, and $\mathrm{\Delta }{r}_{01}=\mathrm{\Delta }{r}_{-10}=1$ for (d). For (a)–(d), $\mathrm{\Omega }=100, c=20$, and $r_{10}^0=r_{0-1}^0=r_{01}^0=r_{-10}^0=1$.

Figure 4

Temporal correlation and optimal paths in self-regulating genes. (a, c) Self-activating and (b, d) self-repressing genes. (a, b) The normalized difference between two-time cross correlations, $A\left(t\right)$, is plotted as a function of $t$ in units of $1/k$. (c, d) The optimal paths for the on-to-off (green) and off-to-on (purple) directions are superposed on the flux $\mathbf{J}$ on the 2D plane of the protein concentration $p$ and the histone state $\zeta$. Calculations were performed using the adiabatic approximation of TF binding or unbinding. The rates of the histone state change are scaled by the parameter $u$ as $r_{10}^0=r_{0-1}^0=r_{01}^0=r_{-10}^0=u$ with $\mathrm{\Delta }{r}_{10}=\mathrm{\Delta }{r}_{0-1}=u$ in (a, c) or $\mathrm{\Delta }{r}_{01}=\mathrm{\Delta }{r}_{-10}=u$ in (b, d). In (a) and (b), $A\left(t\right)$ is plotted with $u=1.5$ (green), $u=1$ (red), and $u=0.5$ (black). In (c) and (d), $u=1$. The other parameters of (a) and (c) are the same as in Fig. 3 and those of (b) and (d) are the same as in Fig. 3.

Figure 5

Landscape, probability flux, and optimal paths of a self-activating gene calculated in the 3D space of the protein concentration $p$, the TF binding status $\xi$, and the histone state $\zeta$, in the nonadiabatic kinetics. (a) Flux $\mathbf{J}$ is superposed on the landscape $U$ with $1<U\le 2$ (light blue), $-0.5<U\le 1$ (red), and $U\le -0.5$ (purple). We find a gene-off state at $p\approx 0, \xi \approx -1$, and $\zeta \approx 0\sim -1$ and a gene-on state at $p\approx \xi \approx \zeta \approx 1$. (b) Flux and optimal paths of on-to-off (green) and off-to-on (purple) directions. The rate parameters are scaled in units of $k$ with $h_0=2$ and $f=0.1$. The other parameters are the same as in Fig. 3.

Figure 6

Landscape $U(p,\zeta )$ and probability flux $\mathbf{J}(p,\zeta )$ with the exogenous TF contribution $p_{\text{ex}}$. (a) $p_{\text{ex}}=0.1$ and (b) $p_{\text{ex}}=0.6$. Calculations were performed on the 2D plane of the protein concentration $p$ and the histone state $\zeta$ in the adiabatic approximation of TF binding or unbinding. The parameters are the same as in Fig. 3.