A mathematical model exhibiting the effect of DNA methylation on the stability boundary in cell-fate networks

Figures & data

Figure 1. The DNA methylation cycle as currently understood [31,32]. (a) The TET-mediated cytosine DNA demethylation cycle starts with transferring a methyl group onto the C5 position of the cytosine to form 5-methylcytosine (5mC) by the DNA (cytosine-5)-methyltransferase (mainly DNMT3A/B for de novo methylation). The 5mC is then further iteratively oxidized by TET enzymes to form 5hmC, 5fC, and 5caC. The 5fC and 5caC are then recognized by TDG that yields an abasic site that is repaired by BER and returns to the unmodified cytosine state. (b) A simplified de-methylation cycle for Oct4 in the PSCC. Due to the fact that 5mC and 5hmC are two leading abundant oxidized forms of cytosine, we can effectively reduce the full cycle to a three state model

Figure 2. A detailed diagram of the PSCC network showing N, T, O interactions. TET, Oct4, Nanog denote the unbound genes. *TET and *Oct4 denote NT binding to the respective promoter sites, TET* and Oct4* denote Oct4 bound to the respective promoter sites, and *TET^{*, *}Oct4* denote both NT and Oct4 binding to the respective promoter sites. The black dashed line indicates binding locations of the Nanog-guided TET protein complex NT, while the blue dashed line indicates binding locations of Oct4. The coloured shaded area emphasizes the DNA demethylation cycle

Figure 3. Graphical illustration of how a shift in the boundary separating two regions of attraction can result in phenotype switches under small perturbations. In the left panel, a perturbation of the silenced state leads to a transient configuration that belongs to the basin of attraction of the active steady state, resulting in a phenotype switch to the active steady state (i.e, the somatic state). In the right panel, the boundary has shifted, reflecting that fact that the relative sizes of the two basins of attraction have exchanged dominance. Now the same perturbation keeps the state in the domain of attraction of the silenced state, and the perturbed state eventually returns to that attractor. The shift in boundary could be due, for instance, to an increase in the methylation rate parameter $\gamma$

Figure 4. The two-pathway activation scheme for TF-gene binding. In the full model, such an activation scheme describes the promoter binding activation for both Oct4 and TET genes. Considering Oct4 as an example, there are two pathways for the activation of its promoter. Each of NT and Oct4 can independently bind to the promoter of Oct4. When both NT and Oct4 bind to the promoter region of the Oct4 gene, it will be activated

Figure 5. The diagram on the left shows how the triangle topology can be further simplified into a two-species CRN. Since the 5hmC methylation state has a slower dynamics compared to the 5mC methylation state, one may employ a quasi-steady states approximation to reduce the triangle topology to the two-state model depicted on the right

Table 1. The CRN model of the PSCC. Reactions have been grouped into eight modules and labelled from R1 to R18

Oct4 promoter(R1-R4)	$\begin{array}{c}{D}_{00}^O+[NT]\underset{a_{nt}}{\overset{K_{nt}\ast a_{nt}}{\rightleftharpoons }}{D}_{10}^O\\ D_{01}^O+[NT]\underset{a_{nt}}{\overset{K_{nt}\ast a_{nt}}{\rightleftharpoons }}{D}_{11}^O\\ D_{00}^O+O\underset{a_O}{\overset{K_O\ast a_O]}{\rightleftharpoons }}{D}_{01}^O\\ D_{10}^O+O\underset{a_O}{\overset{K_O\ast a_O}{\rightleftharpoons }}{D}_{11}^O\end{array}$	TET promoter (R5-R8)	$\begin{array}{c}{D}_{00}^T+[NT]\underset{a_{nt}}{\overset{K_{nt}\ast a_{nt}}{\rightleftharpoons }}{D}_{10}^T\\ D_{01}^T+[NT]\underset{a_{nt}}{\overset{K_{nt}\ast a_{nt}}{\rightleftharpoons }}{D}_{11}^T\\ D_{00}^T+O\underset{a_O}{\overset{K_O\ast a_O}{\rightleftharpoons }}{D}_{01}^T\\ D_{10}^T+O\underset{a_O}{\overset{K_O\ast a_O}{\rightleftharpoons }}{D}_{11}^T\end{array}$
Nanog promoter (R9)	$D_0^N+O\underset{a_O}{\overset{K_O\ast a_O}{\rightleftharpoons }}{D}_1^N$	Nanog-TET Dimer (R10)	$N+T\underset{a}{\overset{K_d\ast a}{\rightleftharpoons }}[NT]$
Protein decay and production (R11-R13)	$\begin{array}{c}N\stackrel{\delta }{⟶}\mathrm{\varnothing };D_1^N\stackrel{{\alpha }_N}{⟶}{D}_1^N+N\\ T\stackrel{\delta }{⟶}\mathrm{\varnothing };D_{11}^T\stackrel{{\alpha }_T}{⟶}{D}_{11}^T+T\\ O\stackrel{\delta }{⟶}\mathrm{\varnothing };D_{11}^O\stackrel{{\alpha }_O}{⟶}{D}_{11}^O+O\end{array}$	Basal Rate (R14-R16)	$\begin{array}{rl}& \mathrm{\varnothing }\stackrel{m_1}{⟶}N\\ & \mathrm{\varnothing }\stackrel{m_2}{⟶}T\\ & \mathrm{\varnothing }\stackrel{m_3}{⟶}O\end{array}$
Methylation (R17)	$D_{00}^O\stackrel{\gamma }{⟶}{D}_m$	Demethylation (R18)	$D_m+NT\stackrel{\theta }{⟶}{D}_{00}^O+NT$

Figure 6. The relative stability (measured by BoAp) of the silenced steady state increases with respect to the DNAm rate ($\gamma$). We show simulations for a representative set of different rates for the full model, where extensive simulations are shown in SI-$\mathcal{\S }$1.1. Each curve in the figure correspond to a particular set of parameters. Each curve ends when bistability is lost and BoAp jumps to 100%. This behaviour is illustrated in SI-$\mathcal{\S }$ 2.1. The simulations also show that, at small $\gamma$ region, the BoAp exhibits a sudden slowing down at around $\gamma ={10}^{-1.5}$

Figure 7. The relative stability (measured by BoAp) of the silenced state decreases with respect to the DNA demethylation rate ($\theta$). Shown are a representative set of simulations for the full model, where extensive simulations are shown in SI-$\mathcal{\S }$1.1. Each curve in the figure corresponds to a particular set of parameters

Figure 8. The BoAp of the silenced state shows a biphasic behaviour when the methylation cycle transitions from slow to fast, where R defined as the methylation association rate R = $\frac{\gamma }{\theta }$. The plot shows the BoAp while $\gamma$ is increased with a fixed methylation ratio, i.e., $\theta$ is increased so as to keep the ratio R constant. The system loses bistability when R exceeds 30. The kinetic rate parameters used are: $K_O$ = 0.3, $K_{nt}$ = 0.2, $K_d$ = 0.1, $a$= 1, $a_{nt}$ = 100, $a_O$ = 100, ${\alpha }_T$ = 1, ${\alpha }_O$ = 1, ${\alpha }_N$ = 1, $\delta$ = 1, $R\in \{0.01,0.1,1,10,20,30\}$, $\gamma \in [{10}^{-2},{10}^2]$

Figure 9. Phase plane analysis of the two-dimensional model (11). (a) The directions of the vector field and the BoAs in the case of bistability with stable steady states $s_0,s_1$, where $s_0$ denotes the silenced steady state and $s_1$ denotes the active steady state. The arrows in each region denote the direction of the vector field. (b) The area of the $\epsilon$-independent BoA of the silenced steady state $s_0$ increases when the methylation ratio increases. Recall that $\epsilon \to \mathrm{\infty }$ denotes the case of a slow methylation kinetics (relative to protein kinetics), while $\epsilon \to 0$ denotes the case of fast methylation kinetics

Kasper Dindler KH. Rasmussen, “Role of TET enzymes in DNA methylation, development, and cancer,”. Genes Dev. 2016;30(7):733–750.

 Google Scholar

Wu X, Zhang Y. TET-mediated active DNA demethylation: mechanism, function and beyond. Nat Rev Genet. 2017;18(9):517.

 Google Scholar