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

ABSTRACT

Cell-fate networks are traditionally studied within the framework of gene regulatory networks. This paradigm considers only interactions of genes through expressed transcription factors and does not incorporate chromatin modification processes. This paper introduces a mathematical model that seamlessly combines gene regulatory networks and DNA methylation (DNAm), with the goal of quantitatively characterizing the contribution of epigenetic regulation to gene silencing. The ‘Basin of Attraction percentage’ is introduced as a metric to quantify gene silencing abilities. As a case study, a computational and theoretical analysis is carried out for a model of the pluripotent stem cell circuit as well as a simplified self-activating gene model. The results confirm that the methodology quantitatively captures the key role that DNAm plays in enhancing the stability of the silenced gene state.

KEYWORDS:

• DNA methylation
• stability boundary
• cell-fate
• transcriptional regulation
• epigenetics
• gene silencing
• basin of attraction
• bistability

Introduction

Cell-fate determination in developmental biology is a multi-step biological process in which cellular functions are modified or specialized, ultimately resulting in a differentiated or proliferative state. A central role in this process is played by the dynamics of a core network of genes, often referred to as a cell-fate network (CFN). The expression levels of these genes can trigger a cascade of events that determine, in principle irreversibly, the fate of a given cell along a specific lineage. The metaphor of marbles rolling down a hill whose shape is shaped by a CFN (the Waddington landscape [1,2]) provides a way to visualize this process. Well-known examples of CFNs are the PU.1/GATA.1 gene regulatory network in haematopoietic progenitors [3,4], the pluripotent stem cell network [4,5], and a circuit including SNAIL, miR-34, ZEB1, and miR-200 that regulates epithelial-mesenchymal transitions in tumour metastasis [6].

Phenotypes associated to CFNs had classically been seen as irreversible. However, in their pioneering work [7], Takahashi and Yamanaka artificially induced a pluripotent state in mouse somatic cells through a process of overexpression of Oct3/4, Sox2, c-Myc, and Klf4. This success notwithstanding, simply overexpressing genes has been experimentally found to be grossly inefficient [8]. Hence, there has been a great interest in obtaining a quantitative theoretical understanding of CFNs in order to guide the process of reprogramming and increase its efficiency [9,10]. Such a quantitative understanding could have a huge impact in the field of regenerative medicine and stem cell therapy [11,12].

Traditionally, CFNs have been modelled within the wider theoretical framework of Gene Regulatory Networks (GRNs) [13,14]. A GRN is defined as a set of genes, each expressing a protein. The expressed proteins can act as Transcription Factors (TFs) by binding to the various promoters in the network to inhibit or enhance expression of the corresponding genes [15,16]. When used to model a CFN, a GRN must be able to display multistability. This means that the long-term expression levels of genes can settle on multiple distinct steady state attractors, each corresponding to a distinct cell lineage. A well-known multistable GRN is the toggle switch [17] whose architecture consists of two mutually inhibiting genes. Similar architectures can also give rise to multistability and have been subject of great interest [2].

Despite their versatility and variety of applications, GRN models of CFNs do not typically account for epigenetic effects, such as DNA methylation, histone modifications, or chromatin remodelling. There have been both theoretical [18–20] and experimental [21] works aiming at the understanding of how epigenetic regulation and gene regulation interact. It is well known that CFNs employ epigenetic regulation as a mechanism to ensure the irreversibility of the cell lineage [22]. DNA methylation (DNAm), in particular, has been well studied in the context of developmental CFNs [23]. DNAm at the transcriptional start sites (TSSs) is a highly heritable, hard to reverse, and robust silencer of genes [24–26]. In comparison to its well-known silencing function at the regulatory regions, less is known about DNAm’s role at intragenic regions [27].

In this work, we develop mathematical models of GRNs that incorporate DNAm at TSSs, and thus can quantitatively explain its effect on gene silencing. As a metric to quantify the gene silencing abilities of DNAm, we consider the shift in the stability boundary of the basin of attraction (BoA) of the silenced steady state. More precisely, we define the BoA percentage (BoAp) of a steady state as the volume fraction of a predetermined region of the state space that is occupied by the BoA. Equipped with this metric, one can then computationally and theoretically study the sensitivity of the BoAp to parameters such as DNAm rates or the time scale of DNAm.

We carry out this program computationally, illustrating it with a three-gene model of the pluripotent stem cell network, and quantifying how DNAm effectively increases the BoAp compared to a standard GRN. In addition, and in order to gain understanding, we also consider an ideal single self-activating gene, for which a more theoretical study is possible. We view our work as a first step towards integrating epigenetic mechanisms into the standard GRN paradigm studied in systems biology. Furthermore, we envision our model (and its derivatives) as computational guides uncovering the possible behaviours of a network and evaluating the feasibility of the interventions that can be used for regulation. For instance, as we will discuss later, the DNAm rate in our model is proportional to the abundance of DNMT, and hence our mathematical framework can quantitatively predict the effect of DNMT on the relative stability of the silenced (methylated) and active steady states. This would allow our model to function as a computational aid to design the reprogramming protocol in the context of the PSCC. We will elaborate further on the biological implications of our model in the Discussion section.

Basic biological concepts and definitions

DNA methylation as epigenetic regulation

DNAm is one of the main epigenetic regulation mechanisms studied in modern molecular cell biology. Methylation plays a crucial role [28] in understanding the dynamics of gene silencing. It is associated with methylation of Cytosine-Phosphate-Guanine (CpG) islands, which are regions of DNA with a high G +C content and a high frequency of CpG dinucleotides relative to the bulk genome [29]. Methylation of promoter regions directly affects TF-promoter binding kinetics, and is often associated with transcriptional repression [24–26,30].

The current understanding of the DNAm cycle is shown in Figure 1. As depicted, cytosines at DNA promoter regions are initially methylated by DNA methyltransferase (DNMT) to become 5-methylcytosine (5mC). The methylated gene promoter is further oxidized to 5-Hydroxy Methylcytosine (5hmC). The oxidized form of the cytosine can further be oxidized by the TET protein into the 5fC and 5caC forms. The cycle closes by a TET/TDG (thymine DNA glycosylase)/BER (base excision repair)-dependent pathway that restores the unmodified cytosine state [31,32].

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

Experiments have shown that hydroxymethylation and formylation are relatively stable cytosine modifications in genomic DNA of both dividing and nondividing cells [33,34]. The 5hmC and 5fC forms in this enzymatic oxidation process are observed to be comparatively longer lived, though still transient, states. Moreover, it also has been shown [31] that the pool size of the transient oxidized cytosine 5hmC is significantly larger than the other forms of oxidized cytosines (5fC, 5caC). Therefore, in this work, we further simplify the DNA demethylation cycle to the cycle in Figure 1(b) which consists of just three major steps, which we will refer to as the ‘triangle topology.’ The cycle starts with an initial de novo methylation on the unbounded gene promoter, which is followed by a TET oxidization enzymatic process oxidizing 5mC to 5hmC. The last stage corresponds to a replication-dependent dilution bringing 5hmC back to the unbound promoter state.

The underlying DNAm dynamics is crucial for understanding how epigenetics regulates CFNs. DNAm plays a central role not only in early embryogenesis, but also in maintaining the correct pattern of methylation along DNA in somatic cells. As a model system of great interest, we review next how the pluripotent stem cell circuit (PSCC) is affected by the DNAm cycle.

The pluripotent stem cell circuit (PSCC)

Some of the key Transcription Factors (TFs) that are identified as being crucial to the PSCC include Nanog, Oct4, Sox2, Klf4, and the TET protein family [7,31]. Many of such TFs are ‘pioneering’ transcription factors that are able to bind directly to the condensed chromatin. According to the induction experiments in [7], the overexpression of the aforementioned TFs can lead to the induction of the pluripotent state. For the PSCC, the role of each TF in the GRN has been extensively studied [35]. Much research has been devoted to elucidating the interplay between the Nanog, Oct4 and TeT1 proteins [36,37]. Oct4 is central to the pluripotency machinery, has been very well studied [38], and it can by itself reprogram a differentiated cell into pluripotency [39]. Other experiments regarding the induction of the pluripotent state in PSCCs have highly heterogeneous results [8,40]. More recent research [41,42] has shown that, at a mechanistic level, Nanog is able to guide TET protein to bind to a particular region of the methylated region on the chromatin via the formation of a Nang-TET compound that is actively involved in the DNA demethylation process.

Notably, DNAm of the Oct4 gene plays a key role in the regulation of the PSCC. In Yamanaka’s lab early works [7,43], it was observed that the promoters of Oct4 (and other pluripotency-related genes) are unmethylated in embryonic stems (ES) cells, while they are methylated in mouse embryonic fibroblasts. Furthermore, it was observed that Oct4 becomes demethylated just prior to reprogramming [44], and that the differentiation of ES cells is associated with the silencing of Oct4 and de novo methylation of its regulatory region [45,46]. Hence, activation of the PSCC requires demethylation of the corresponding promoters [47], which means that understanding DNAm is a prerequisite of understanding the PSCC. However, very few mathematical models have incorporated DNAm in their GRNs. In this work, we model the interplay between epigenetic regulation via the Nanog-guided TET-mediated active demethylation cycle and the role of pioneering transcription factors such as Oct4 and Nanog. A mechanistic view of this interplay helps one understand how emergent cell states can arise phenotypically from the underlying genotypic level interaction. Although the full modelling of epigenetic regulation at the molecular level would involve many additional epigenetic regulation mechanisms beyond methylation such as histone modifications [31], DNAm plays a major role in the coarse grained picture of the process of the induction of the pluripotent state which we are mostly interested in.

Basin of attraction percentage (BoAp)

Using high-dimensional attractors, and specifically stable steady states, to represent cell fates in the epigenetic landscape is a widely used paradigm for studying cell fate [14,48]. In a biological context, a steady state attractor corresponds to an observable phenotype, e.g., a somatic state with a silenced Oct4 promoter. In the dynamical cellular environment, perturbations can temporarily change the expression levels away from steady-state. Intuitively, the GRN should be able to return to its attractor (phenotype) for a sufficiently small perturbation, while there is a risk of jumping to another attractor (phenotype) for a large perturbation. By ‘perturbation’ we mean an abrupt change in the levels of the dynamical variables of interest such as TF abundances and DNAm levels. An attractor that can withstand a large perturbation is more ’stable’. Hence, to characterize the stability boundary of such attractors we use the concept of Basin of Attraction (BoA). The BoA of a given attractor (phenotype) is the set of all perturbed levels that will return back eventually to the steady state attractor. Figure 3 illustrates these ideas.

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$

The concept of BoA has been frequently used by researchers and it is standard in dynamical systems when analysing multi-stable systems; moreover, it has recently been systematically used by the physics community in that context [49]. In order to quantify the BoA of a given attractor, we propose to use the concept of BoA Percentage (BoAp), which is defined as the ratio of the ‘size’ of the BoA of interest to the total size of the space of interest. A larger BoAp for a given attractor means that its BoA is larger, and hence it is relatively more stable, while a smaller BoAp means that it is relatively less stable.

In our simulation analysis, we adopt this concept to quantify the stability boundary of the somatic cell state. In a bistable one-dimensional model, the BoAp of the smaller steady state $s_0$ can be defined as the distance between $s_0$ and the unstable state relative to the length of a segment of interest. However, this is not an accurate measure for high-dimensional systems. Therefore, we adopt a volumetric definition for the BoA. We start by fixing a region of interest in the state space. With respect to such a region, BoAp can be defined as follows:

(1) $\mathrm{B}\mathrm{o}\mathrm{A}\mathrm{p}:=\frac{\mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{B}\mathrm{o}\mathrm{A}\mathrm{o}\mathrm{f}{s}_0}{\mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}.$(1)

Equipped with this definition, we can study the effect of DNAm and demethylation constants on BoAp. The silenced steady state corresponds to a highly methylated promoter. Hence, as the methylation rate $\gamma$ increases, the BoA of the silenced steady state expands. On the other hand, as the demethylation rate $\theta$ increases, the BoA of the silenced steady state shrinks. Our subsequent results will show that the BoAp of the silenced steady state increases as the methylation rate increases, and decreases as demethylation rate increases. Note that we will use the term ‘steady state’ to mean a ‘steady state attractor’ from now on.

Results

The PSCC with epigenetic regulation

In our full model description of the PSCC, the network can be viewed as made up of two parts. The first part is the gene-TF regulatory network [50] which is a core underlying network defining how TFs bind and unbind to DNA as shown in Figure 2 without the shaded area. The main regulatory sub-network for the activation of the Oct4 and Nanog genes is depicted within the two-pathway activation scheme in Figure 4. Both the Oct4 and the TET gene promoter binding mechanisms will follow such scheme, which has two pathways that lead to the activation of the gene. Oct4 protein self-binds to the first site while the Nanog-guided TET protein complex binds to the second site. The promoter is active only when both sites are occupied. In this GRN, the Nanog gene is activated by the Oct4 gene. Furthermore, Nanog and TET form a heterodimer which models TET guiding Nanog to the target promoters.

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

The GRN with DNA demethylation is summarized in Figure 2. The DNA demethylation cycle with the triangle topology can be further simplified into two states, which are shown in Figure 5. The parameter $\gamma$ is thought of as the effective methylation rate and the parameter $\theta$ as the effective demethylation rate. The CRN for the PSCC with a two-state methylation cycle can be written as in Table 1. (For simplicity, we omit explicit variables for mRNA and other intermediate steps such as protein maturation and post-translation modifications.) We model this network deterministically through a system of ODEs. The variables in this system can be thought of as mean numbers for the various species in a population. The full model is described in the Methods section.

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$

Model reduction via quasi-steady state approximation

The full model CRN is computationally intensive to simulate, and it is conceptually difficult to interpret its dynamical behaviour. Hence, we perform a quasi-steady state approximation by distinguishing between two time-scales: the ‘fast’ time-scale of TF-promoter kinetics (binding/unbinding), and the ‘slow’ time-scale of protein expression/decay and DNA methylation kinetics. Note that the methylation cycle in our model is a phenomenological reduction of the cycle depicted in Figure 1, which includes a dilution step as well as the coordinated action of many TFs and enzymes. Hence, the methylation cycle is not faster than protein kinetics, which are known to be much slower than TF-promoter kinetics [15]. This is also supported by experimental evidence, as the re-establishment of methylation patterns can take up to 20 minutes [51], and de novo methylation kinetics take days which are even slower than protein kinetics as has been observed in various de novo methylation studies [52,53]. It can be noted also that inherited methylation patterns have very slow kinetics in general [54,55].

Therefore, we assume that the reaction rates of the promoter-TF binding and unbinding dynamics are faster than the other reaction rates in the CRN Table 1. By applying time scale separation or quasi-steady state approximation, we reduce the full dynamical system model from a 17-dimensional state space to a 4-dimensional state space. The full model ODEs are given in the Methods section. The reduced dynamical system’s ODEs can be written as below:

$\frac{dN}{dt}=m_1-\delta N+\frac{{\alpha }_{N}O}{K_O+O},$

$\frac{dT}{dt}=m_2-\delta T+\frac{{\alpha }_{T}NTO}{\left(NT+K_{nt}{K}_d\right)\left(K_O+O\right)},$

$\frac{dO}{dt}=\frac{-K_{nt}{K}_d\delta O^2-\delta NT{O}^2+K_O{K}_{nt}{K}_d{m}_3+{\alpha }_{O}NTO+K_{nt}{K}_d{m}_{3}O}{\left(NT+K_{nt}{K}_d\right)\left(K_O+O\right)}$

$+\frac{K_O{m}_{3}NT+m_{3}NTO-K_O\delta NTO-K_O{K}_{nt}{K}_d\delta O-{\alpha }_O{D}_{m}NTO}{\left(NT+K_{nt}{K}_d\right)\left(K_O+O\right)}$

$\frac{d{D}_m}{dt}=\frac{-K_O{K}_{nt}{K}_d^2\gamma D_m+K_O{K}_{nt}{K}_d^2\gamma -\theta K_O{N}^2{T}^2{D}_m-\theta N^2{T}^{2}O{D}_m}{K_d\left(NT+K_{nt}{K}_d\right)\left(K_O+O\right)}$

$-\frac{\theta K_O{K}_{nt}{K}_{d}NT{D}_m+\theta K_{nt}{K}_{d}NTO{D}_m}{K_d\left(NT+K_{nt}{K}_d\right)\left(K_O+O\right)}$

Effect of methylation rate

For the full model, Figure 6 displays how the BoAp changes as one changes the epigenetic reaction rate $\gamma$. Each plot in the figure corresponds to one set of parameters, sampled from a wide range. The parameter set for the full model contains the following rates: $K_O$, $K_{nt}$, $K_d$, $a$, $a_{nt}$, $a_O$, ${\alpha }_T$, ${\alpha }_O$, ${\alpha }_N$, $\delta$, $\gamma$, $\theta$, $m_1$, $m_2$ and $m_3$. Biologically speaking, changing the reaction rate $\gamma$ is equivalent to changing the effective methylation rate. Our subsequent study of the single gene model will display the same trend of BoAp’s dependence on $\gamma$.

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}$

The phenotypic somatic state is characterized [36] by low expression of Oct4, Nanog, TET, and methylated promoters. With simulations performed over 6200 parameter sets as shown in Supplementary Figure 1(a), the BoAp vs. $\gamma$ plot shows that the BoAp is monotonically increasing with respect to $\gamma$. In Figure 6, we show a representative sample of these curves. If $\gamma$ is increased beyond a certain threshold, then bistability is lost and the BoAp of the silenced steady state becomes 100%. Note the different behaviour for intermediate values of $\gamma$ compared to low or larger values of $\gamma$. Another interesting observation is that the BoAp is bounded below by $10\mathrm{\%}$ for all parameter sets when $\gamma$ is greater than ${10}^{-1.5}$.

On the other hand, changing the reaction rate $\theta$ means changing the demethylation rate. As we reduced the triangle topology for the DNA demethylation cycle to a two-state cycle in the full model, the effective demethylation rate can be thought of as representing the combination effect of 5mC to 5hmC oxidation rate affected and the rate of oxidized 5hmC reverting to the unbound promoter state via either TDG-BER pathway or dilution via cell replication. In Supplementary Figure 1(b), the simulation shows that the BoAp is monotonically decreasing with respect to $\theta$. Over 6200 parameter sets have been simulated and we plot BoAp vs. $\theta$ over 4 orders of magnitude to show the trend. In Figure 7, we show a representative sample of these curves. We choose the value of the range for $\theta$ to be the same as the range of $\gamma$ in the BoAp vs. $\gamma$ plot, which both are from ${10}^{-2}$ to ${10}^2$.

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

These results confirm that our model captures the key role that methylation rate plays in enhancing the stability of the somatic steady state. The higher the effective methylation rate, the larger the BoAp. This means that it is more likely for the underlying dynamical system to stay at the somatic (silenced) steady state, which implies that it is relatively more stable. Furthermore, the methylation rate can be increased to the point of eliminating the existence of the active steady state (the pluripotent steady state) and making the silenced steady state immune to any abrupt perturbations. Note that the methylation rate is proportional to the abundance of DNMT, and hence all the conclusions above can be interpreted as applying to the DNMT abundances. We will elaborate on the practical implications of our model further in the Discussion section.

Effect of slow/fast methylation kinetics

To understand how the relative methylation rate affects the BoAp, we define the methylation association rate as R = $\frac{\gamma }{\theta }$.

The simulations in Figure 8 show that BoAp vs. $\gamma$ across 4 orders of magnitude in $\gamma$ while $\theta$ is increased so as to keep the ratio $R$ constant with different ratios of $R$ that span 3 orders of magnitude. Assuming a normalized protein decay rate ($\delta =1$), Figure 8 allows us to disentangle two effects: the methylation ratio ($R$) and the time-scale of the methylation cycle. The plots display a biphasic trend, i.e., the BoAp monotonically increases and then monotonically decreases. Therefore, there is a maximum value of the BoAp which occurs in the slow range $[{10}^{-1.5},{10}^{-0.5}]$ for various methylation ratios $R$, while the BoAp decreases as the methylation cycle becomes faster (which is consistent with the analytical picture to be shown in Figure 9). These results mean that the stability of the silenced state (or the efficiency of methylation) is maximized with a higher methylation ratio and a relatively slower methylation kinetics (roughly an order of magnitude slower than the protein decay rate).

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

It should be noted that the biphasic nature of the trends is not a result of high-dimensional dynamics. Rather, the same behaviour can be observed in the single gene model (see SI-$\mathcal{\S }$ 1.2). Furthermore, the biphasic trend for the fixed ratio $R$ can roughly be inferred from the PSCC simulations by varying each $\gamma$ and $\theta$ independently. As per the trend we have shown in Figure 6, the BoAp increasing rate exhibits a sudden slowing down at $\gamma$ around ${10}^{-1.5}$, while there is no corresponding change for the demethylation rate $\theta$ shown in Figure 7. Therefore, it is not surprising that the methylation association ratio $R$ has a biphasic behaviour.

Model of a single self-activating gene

In order to complement the numerical investigations in the previous sections, we consider next a simplified model which includes Oct4 only. Its interaction with the rest of the network is modelled as a self-activating loop. This model is more amenable to an analytical study as will be discussed below.

Our simplified model consists of a single self-activating gene. A CRN model for it can be written as follows:

(2) $2X+D_i\underset{{\beta }_{-}}{\overset{\beta /2}{\rightleftharpoons }}{D}_a$(2)

(3) $D_a\stackrel{k}{⟶}{D}_a+X$(3)

(4) $D_i\underset{\gamma }{\overset{\theta }{\rightleftharpoons }}{D}_m$(4)

(5) $X\stackrel{\delta }{⟶}\mathrm{\varnothing },$(5)

where $X$ is the protein expressed from the active promoter state $D_a$ and $D_i$ is the inactive promoter state. (As earlier, we omit explicit variables for mRNA and other intermediate steps such as protein maturation and post-translation modifications.) To model self-activation, we assume that an $X$ homodimer can bind to the inactive promoter site to activate the gene. The inactive promoter region can be methylated to a state that we denote by $D_m$. Reversely, $D_m$ can be epigenetically modified back to the $D_i$ state. We assume that protein $X$ has degradation rate $\delta$. Our aim is to study the effect of the methylation ratio $\gamma /\theta$ on the dynamics.

Since the total concentration of promoters is conserved, the ODE model can be reduced by writing $D_i=1-D_a-D_m$, where we scale the total number of promoters to 1. The system of ODEs that results is given as follows:

(6) $\frac{dx}{dt}=k{D}_a-\delta x-\beta x^2(1-D_a-D_m)+{\beta }_{-}{D}_a,$(6)

(7) $\frac{d{D}_a}{dt}=\beta x^2(1-D_a-D_m)-{\beta }_{-}{D}_a,$(7)

(8) $\frac{d{D}_m}{dt}=\gamma (1-D_a-D_m)-\theta D_m.$(8)

The steady states can be evaluated analytically to yield a cubic equation. There always exists a steady state in which $x$$=$$0$ (and also $D_a$$=$$0$). We call such a state the silenced steady state. It is easy to verify by linearization that this state is locally stable, for any parameter values. A bifurcation occurs and three steady states exist if and only if the following inequality holds:

(9) $\frac{k^2}{{\delta }^2}>\left(\frac{\gamma }{\theta }+1\right)\frac{4}{B},$(9)

where $B:={\beta }_{-}/\beta$ is (double) the dissociation ratio, and $k/\delta$ and $\gamma /\theta$ are the production and methylation ratios, respectively. One of two additional steady states is always locally asymptotically stable and the other one is unstable (a saddle). We call the second stable steady state the active steady state. It follows from condition (9) that a sufficiently high methylation ratio renders the active steady state non-existent. Hence, a point attractor will coincide with the silenced steady state. This recapitulates our observation in Figures 6, 7 that increasing the methylation ratio beyond a certain threshold will eliminate bistability, and hence the active (pluripotent) steady state disappears.

The model can be reduced further by utilizing the fact that transcription factor binding/unbinding is fast relative to protein expression and epigenetic modifications. We thus carry out a quasi-steady state (QSS) approximation, in which we set the state $D_a$ to:

(10) $D_{a,qss}=\frac{(1-D_m)x^2}{B+x^2}.$(10)

Substituting this expression into the ODE system, we obtain a two-dimensional system whose two variables are the promoter methylated state $y(t):=D_m(t)$ and the protein $x$, as follows:

(11) $\dot{x}=k(1-y)\frac{x^2}{B+x^2}-\delta x$(11)

$\epsilon \dot{y}=\gamma (1-y)\frac{B}{B+x^2}-\theta y,$

where $\epsilon$ is a parameter that represents the relative time-scale separation between protein and methylation kinetics. This system must be studied for $(x,y)\in [0,\mathrm{\infty })\times [0,1]$. We next analyse this reduced system.

Global behaviour

Observe that all solutions of the reduced system are bounded. Indeed, $y$ is bounded by definition, and $\dot{x}<0$ whenever $x>k/\delta$, so the rectangle $[0,k/\delta ]\times [0,1]$ is an attractor.

The steady states of the reduced system are in a one-to-one correspondence with the states of the full system, in the sense that for each steady state $(x,y)$ there is a unique steady state $(x,D_a,D_m)$ of the full system, where $D_m=y$ and $D_a$ is given by the expression in Eq. (10)(10) $D_{a,qss}=\frac{(1-D_m)x^2}{B+x^2}.$(10) . Thus, there are one or more steady states depending on condition (9). We may think of the system (11) as a ‘toggle switch’ in which $Y$ (the methylated promoter) represses $X$ (protein), and conversely $X$ represses $Y$. If silencing is stronger, we expect the protein state to remain at zero, while if activation is stronger, we expect protein values to converge to a higher steady state.

In order to analyse the global behaviour of system (11), we first evaluate its Jacobian matrix:

$\left[\begin{array}{cc}\frac{2kBx(1-y)}{{(B+x^2)}^2}-\delta & -k\frac{x^2}{B+x^2}\\ \frac{-2\gamma Bx(1-y)}{\epsilon {(B+x^2)}^2}& -\frac{1}{\epsilon }\left(\frac{\gamma B}{B+x^2}+\theta \right)\end{array}\right].$

Both off-diagonal entries are non-positive, which makes this system a monotone system, and specifically, a system that is cooperative with respect to the cone defined by the $(-,-)$ orthant, as discussed in [56] (see also [57] for exposition and biological applications). In [58], M.W. Hirsch proved (Theorem 2.2) that every bounded solution of a planar cooperative system converges to an equilibrium; see also Theorem 2.2 in [56]. Coupled from boundedness of solutions for our system, we conclude that every trajectory converges to an equilibrium (which depends on the initial conditions of the trajectory).

Phase plane analysis

Our two-dimensional reduced model can be analysed via nullclines. A summary of the analysis is depicted in Figure 9 for a bistable model. The nullclines and the steady states are independent of $\epsilon$. On the other hand, the separatrix, which is the boundary between the BoAs, is $\epsilon$-dependent. It is shown that each BoA can be partitioned into $\epsilon$-dependent and $\epsilon$-independent regions (Figure 9(b)) shows that the area of the $\epsilon$-independent BoA of $s_0$ increases as the methylation ratio is increased. Hence, it can be argued that the BoAp of the silenced steady state increases as the methylation ratio increases, which recapitulates our observations for the full model. Furthermore, we have performed numerical simulations for the single-gene model and they show the same trends. Refer to SI-$\mathcal{\S }$ 1.2.

Slow methylation kinetics

We have argued previously that DNAm can be even slower than protein kinetics; hence, we focus on case $\epsilon \to \mathrm{\infty }$ in this subsection. Figure 9 shows that, in this limit, the BoA of $s_0$ consists only of the $\epsilon$-independent region. Next, we analytically investigate the aforementioned limit.

Computing the asymptotic separatrix

Setting $\dot{x}=0$, we see that the $x$-nullcline consists of two components, which we call the silenced and active arcs, because the dynamics along each arc converge to the silenced and active steady states, respectively. The arcs can be written as follows:

1. Silenced/Somatic: The arc is $x=0$, which gives the slow dynamics of

$\dot{y}=\gamma B(1-y)-\theta y.$

This is a linear system, and it has an asymptotically stable steady state at $y=\frac{B\gamma /\theta }{B\gamma /\theta +1}$. Notice that the location of this steady state is an increasing function of the methylation ratio $\gamma /\theta$.

2. Active/Pluripotent: The arc is given by the implicit equation: $\frac{k}{\delta }(1-y)x=B+x^2$. Solving for $x$ in terms of $y$:

(12) $x=\frac{k}{\delta }(1-y)\pm \sqrt{\frac{k^2}{4{\delta }^2}{(1-y)}^2-B}.$(12)

The slow dynamics are given as (note ${lim sup}_{t\to \mathrm{\infty }}y(t)\ne 1$ for all $t$ and any initial condition with $y(0)\ne 1$):

$\dot{y}=\frac{\gamma \delta B}{k}-\theta y.$

Hence, $y=\frac{\gamma \delta B}{\theta k}$ is the asymptotic steady state, which is proportional to the methylation ratio $\gamma /\theta$.

Given an initial condition $(x(0),y(0))$, we want to determine which arc will be representative of the dynamics. A necessary condition for the existence of the second arc (i.e., the chance to converge to the active steady state) is that the quantity under the square root in (12) is real. Hence, we have:

(13) $y(0)\le 1-\frac{2\delta \sqrt{B}}{k},\mathrm{a}\mathrm{n}\mathrm{d}1-\frac{2\delta \sqrt{B}}{k}>0$(13)

Hence, $y(0)>1-\frac{2\delta \sqrt{B}}{k}$ automatically implies that the silenced arc represents the dynamics regardless of the initial $x(0)$. Furthermore, if the dissociation ratio is high compared to the protein production ratio (i.e., $k/\delta <\frac{1}{2}\sqrt{B}$) then the dynamics are also represented solely by the silenced arc. Assuming the necessary condition is satisfied, we need also that

(14) $x(0)>\frac{k}{2\delta }(1-y)-\sqrt{\frac{k^2}{4{\delta }^2}{(1-y)}^2-B}.$(14)

To summarize, if both (13), (14) are satisfied, then the second arc is representative of the dynamics; otherwise, the first arc represents the dynamics. This conclusion is summarized in the dark grey region in Figure 9(a).

A remark concerning reprogramming by over-expression of $X$. Our global analysis allows us to intuitively understand the role of $X$ in over-expression experiments (such as the ones carried out in Yamanaka’s lab [7]). In this scenario, we assume that $(x(0),y(0))=(0,\frac{B\gamma /\theta }{B\gamma /\theta +1})$, i.e., the initial condition is a silenced steady state. Over-expression means that a sudden dose of $x$ will be injected to the system such that $x(0_{+})=u$, where $u$ is a constant positive number. Our aim is to drive the system to the BoA of the active steady state.

It follows from the analysis above that

$\frac{B\gamma /\theta }{B\gamma /\theta +1}\le 1-\frac{2\delta \sqrt{B}}{k}$

is necessary for reprogrammability. Hence, the expression ratio $\frac{k}{\delta }$ and the methylation ratio $\frac{\gamma }{\alpha }$ must satisfy:

(15) $\frac{k}{\delta }>2\sqrt{B},\frac{\gamma }{\theta }\le \frac{1}{B}\left(\frac{k}{2\delta \sqrt{B}}-1\right).$(15)

If both of these conditions are satisfied, then overexpression will steer the trajectory to the active arc if the following holds:

(16) $u>\frac{k}{2\delta }\left(\frac{1}{B\gamma /\theta +1}\right)-\sqrt{\frac{k^2}{4{\delta }^2}{\left(\frac{1}{B\gamma /\theta +1}\right)}^2-B}.$(16)

In a nutshell, eqns. (15)(15) $\frac{k}{\delta }>2\sqrt{B},\frac{\gamma }{\theta }\le \frac{1}{B}\left(\frac{k}{2\delta \sqrt{B}}-1\right).$(15) , (16(16) $u>\frac{k}{2\delta }\left(\frac{1}{B\gamma /\theta +1}\right)-\sqrt{\frac{k^2}{4{\delta }^2}{\left(\frac{1}{B\gamma /\theta +1}\right)}^2-B}.$(16) ) imply that we need a sufficiently high Oct4 expression ratio, a sufficiently high overexpression value ($u$), and a sufficiently low methylation ratio. This again underscores the importance of the methylation ratio.

Comparison of the two models. The results for the single gene model (presented analytically above, and numerically in SI-$\mathcal{\S }$ 1.2 parallel those for the four-dimensional model (which was, in turn, obtained through QSSA from the full 17-dimensional CRN)). For example, we see a qualitatively similar behaviour regarding the BoAp’s dependence on the (de)methylation rate in both models. In addition, the results in the reduced model with 4D state space again agree with those for the single gene model with 2D state space.

The full model has the advantage of incorporating all the regulation loops and displaying the effect of the methylation cycle on the steady-state values of the Nanog and TET proteins. However, the single gene model can be reduced to a 2D model which enables the analytical computation of the steady states and the application of phase plane methods. Basic trends in the BoAp can be intuitively read from the shift in the separatrix within the predefined cube. Unfortunately, performing similar analysis for the 4D model is intractable, but the numerical simulations have shown the same basic behaviour predicted by the single gene model.

The single gene model is not equipped with the additional regulation loops mediated by the NT heterodimer. Nevertheless, this model demonstrates that the trends displayed by the full model are intrinsic effects of the proposed methylation/demethylation cycle. In both cases, a slow epigenetic process – in our case, a methylation cycle – shifts the stability boundary of the underlying system consistently in the same direction.

Discussion

We have proposed a mathematical modelling paradigm to elucidate the role of DNAm in the activation and silencing of genes. This has been achieved by integrating DNAm into the standard paradigm of GRNs, which rely on TF-promoter binding/unbinding. Our model can be considered as an effective model which considers essential chemical reactions that are relevant to the core cell Oct4 reprogramming network. We have provided a quantitative model of the effect of methylation rates on gene silencing and activation, particularly in the context of ES cell reprogramming. We have shown that increasing the methylation rate increases the BoAp of the silenced state until the active state disappears and the silenced state becomes the only stable steady state. Since the methylation rate is directly proportional to the abundance of DNMT, our results predict that reducing DNMT abundances will lead to gene activation via the emergence of the new active steady state and the loss of stability of the silenced state by BoA shrinkage. Indeed, experimental evidence shows that the application of DNMT inhibitors leads to the reactivation of the associated genes in various contexts [59–62]. In the context of the PSCC, repressing DNMT3B helps in generating pluripotent cells and underscores its role as a regulator of the circuit [63]. Furthermore, downregulation of DNMT3a/b via miRNAs has been suggested as a mechanism for regulating pluripotency [64,65]. On the other hand, DNMTs 1/3a/3b genes can be eliminated without any negative impact on the pluripotent state [66]. Similar results were reproduced in human ES cells with DNMT 3a/b knockout [67]. These results are consistent with our prediction that a highly reduced methylation rate will consolidate the stability of the active steady state. As discussed in the introduction, a smaller basin of attraction of the silenced state means that natural perturbations in the cell environment will have a greater chance of leading to a state transition to the active state.

Since DNAm and demethylation rates can be changed dynamically during the developmental process, the model can be applied at different stages. Also, it can be used to predict the intervention needed (for instance in terms of DNMT up- or down-regulation) to activate or silence a gene of interest. The relative stability of a certain phenotype (measured by the BoAp) can be also estimated which correlates with the average time (e.g., number of days) needed to transition between different phenotypes. Hence, our model can be used as a quantitative guide to design the reprogramming protocol in the context of the PSCC. Furthermore, the DNAm rate can be measured directly. Experimental tools have been developed to trace the dynamic changes of DNAm at single base resolution [68]. In principle, a time series of the DNAm density at the TSS over the course of the developmental process can be constructed from such data, and the methylation rate can be determined accordingly. Such time-series data have been used in a recent experiment [53], which studied the reprogramming of B cells to induced pluripotent cells. That paper showed that the last stage of the successful pluripotency induction displays successive waves of hydroxymethylation associated with a decrease in DNA methylation and active demethylation. Although such techniques are in their infancy, we predict that our model can be used, in conjunction with more extensive data collection methods, to predict the likelihood of the activation of a gene based on the measurable methylation and demethylation rates.

Our model can also be used in conjunction with more direct experimental manipulation methods. In addition to DNMT inhibitors mentioned above, a recent intriguing technique in synthetic epigenetics is the newly developed synthetic read-write modules [69] implemented in human cells via m6A DNA modification. Such synthetic epigenetic regulatory system can allow us, in principle, to regulate the DNAm patterns and influence the resulting phenotypes by controlling the (de)methylation rate. Using our model, the relative stability of a phenotype is reflected in the duration of the time that a cell takes to successfully reprogram to the pluripotent state due to the stochastic nature of the underlying processes. From an engineering perspective, with the help of synthetic epigenetic read-write modules, one can aim at re-engineering the DNAm patterns to explore how changes of (de)methylation rates affect the reprogramming time during the developmental process. In addition to directly manipulating the DNAm or demethylation rates, one can also think of indirect ways for interacting with DNAm. Evidence suggest that histone modifications like H3K4me1/3, H3K27me3, H3K27ac and the polycomb family of proteins are intricately interacting with DNAm and are playing key roles in recruiting some co-factors and interacting with enhancer regions [70]. Therefore, one can experimentally engineer specific pathways for perturbing the DNAm rate through histone level interventions.

Along with the discovery of the role of TET protein in regulation of the DNA (de)methylation cycle, several theoretical models have been proposed to understand the underlying mechanistic picture of how TFs and DNA interact with each other, and additionally how epigenetic factors such as DNAm and histone modifications affect gene regulatory networks in various biological contexts. In the following, we discuss the relation between these models and ours. Many related core gene regulatory networks have been proposed at various levels of complexity. An early attempt using a probabilistic Boolean network that describes the interplay between gene expression, chromatin modifications, and DNAm has been proposed in [19]. In addition to lacking mechanistic interpretations, Boolean models do not allow one to quantify the effect of DNAm on the stability boundary between the multiple phenotypes. As the role of the TET protein family in gene regulation became clearer [31,42,71,72], various CRN models and ODE models that have partially considered TET-mediated demethylation cycle [36,73] and core gene-TF regulatory mechanisms have been proposed. Notably in [73], an ODE model for a core gene-TF regulatory network of the PSCC has been studied. The system involves Oct4, Sox2 and Nanog proteins and exhibits multistability under certain conditions. Regarding epigenetic regulation, that work modelled the DNAm cycle on top of Oct4 auto-activation. The model did not consider the interaction between Nanog and TET, nor how NT heterodimer regulates the DNA demethylation cycle. Despite the fact that the DNAm model in [73] did not include a two-pathway activation-like mechanism for gene activation, it represented an early attempt to mechanistically link gene regulation and chromatin modifications. In [36], a mechanism was proposed for the core gene-TF network without epigenetics, and then a two-state model without NT regulating the DNA demethylation process was investigated on top of this core network. However, the underlying chemical reaction model in [36] was not spelled-out clearly. In addition to this, without the NT complex involved in the DNA demethylation process, the model in [36] couples the core gene-TF network and DNA demethylation cycle, which does not reflect the experimental evidence of the role of NT in the regulation of DNA demethylation. Specifically, the model lacks the reaction in which NT complex binds to the 5hmC. In the recent experimental work [74], four chromatin regulators, including DNMT3B which causes de novo methylation of cytosine-guanine dinucleotides (CpGs), have been studied. Though a descriptive three-state model was proposed in [74] for explaining how each chromatin regulator affects the silencing and reactivation of gene expression, a detailed mechanism at a molecular level was not provided. In comparison with all these previous attempts, we model the underlying CFN with a detailed CRN model at a molecular level, and we clearly show how epigenetic regulation, specifically the DNA demethylation cycle, regulates the CFN by the Nanog-guided TET protein complex.

In [75], a purely phenomenological model of how epigenetic feedback affects gene expression dynamics is proposed. The main idea is to represent the kinetic rate parameters by epigenetic variables that have their own dynamic equations. More recently, a similar idea has been employed in the context of epithelial-mesenchymal transition (EMT) and cancer metastasis [18]. Such a view of epigenetic feedback regulation of kinetic rates is a valuable complementary perspective.

So far, we have explored a cell fate network combining gene-TF interactions and a TF-guided DNA demethylation process. For the CRN we proposed, we have considered how the TET-mediated demethylation cycle potentially affects the gene regulatory network, and subsequently affects the cell fate decision associated with the transcription factors. We have simplified the demethylation cycle into a two-state cycle. Our mechanistic modelling approach using CRN agrees with the intuition on how DNAm affects the cell fate network. The faster the methylation rate the more stable the silenced (somatic) state will be. However, in the detailed molecular picture of the cell fate network that involves Oct4, Nanog, TET and potentially many other transcription factors, there might exist all sorts of self or non-self interactions in the network. For example, if one adds Nanog protein self activation, one would expect such a model to potentially exhibit tristability. Although our model explored the simplest case in which the system gives rise to bistability, the basin of attraction boundary of other models with tristability can still be studied by the relative BoAp within a predefined hypercube.

DNAm’s role in gene silencing at transcriptional start sites has been well established where genomic maps have shown that regulatory genes are indeed active when they are unmethylated [76]. Furthermore, its repressing roles have been documented in many specific cases such as that of the PSCC [46], X-chromosome inactivation [77], and in ageing and cancer [78,79]. However, it is worth mentioning that a full and universal understanding of its role across different genes and its effect on phenotypic variability is not completely achieved yet. This is since different epigenetic mechanisms at different scales are intricately related to DNAm pattern and dynamics. Additionally, its role can vary in different localized biological contexts [76]. For example, although DNAm generally represses TF binding at transcriptional start sites, methylation in the gene body can have a different role [25]. For instance, it has been shown recently that DNAm can indirectly facilitate TF binding in some special cases [80]. Another issue for debate is whether DNAm is a master regulator or a by-product of the underlying GRN [80,81]. Therefore, the analysis provided in this paper is more predictive for the majority of cases in which DNAm plays a repressing role in a GRN. A future direction will be to develop a more encompassing framework for the various roles of DNAm.

As for future directions, even though the two-pathway activation is a mechanism-based model, it does not capture the full complexity of the whole picture of gene activation. The implicit assumption for this activation process is that genes are only activated by TFs binding to promoters. Hence, we are studying the extension of our methods to the analysis of more sophisticated transcriptional control mechanisms, such as the formation of super-enhancers [82]. We are also interested in combining machine learning approaches with our CRN model so as to discover new mechanisms when the relevant data is available. Finally, when generalizing the BoAp to a stochastic setting, the effective landscape of the dynamical system and the mean first passage time between equilibria are the concepts which are most relevant to the analysis of the stochastic stability boundary. Furthermore, stochastic analysis via effective potentials can be used for defining the BoAp and exploring gene-TF interactions with epigenetic regulation.

Methods

Basic mathematical concepts and definitions

Mathematically, our model is initially described as a Chemical Reaction Network (CRN). CRNs provide a natural formalism for representing biological interactions, and in particular biochemical processes, and CRN descriptions map into systems of Ordinary Differential Equations (ODEs) using standard procedures.

Chemical reaction network (CRN) framework

It is essential for our study to keep track of promoter occupancy, since DNAm at the level of promoters is commonly understood as a slow process compared to TF binding/unbinding, and it can be even slower than protein kinetics. However, our main model assumes that DNAm and protein kinetics are in the same time-scale.

Following [83], each gene is associated with promoter states and one or more a protein states.¹ After the model is set up, we later reduce it by the substitution of quasi steady-state approximations of fast variables, based on appropriate relative time-scales.

We briefly review the general formalism, see, e.g., [84,85]. A CRN is specified by a set of species $S=\{Z_1,..,Z_n\}$ and a set of reactions $R=\{{\mathrm{R}}_1,\dots ,{\mathrm{R}}_{\nu }\}$. A reaction ${\mathrm{R}}_j$ can be written as: ${\sum }_{i=1}^n{\alpha }_{ij}{Z}_i\to {\sum }_{i=1}^n{\beta }_{ij}{Z}_j$. The associated stoichiometry matrix $\mathrm{\Gamma }\in {\mathbb{R}}^{n\times \nu }$ is defined elementwise as $[\mathrm{\Gamma }{]}_{ij}={\beta }_{ij}-{\alpha }_{ij}$. Each reaction ${\mathrm{R}}_j$ can occur with a rate function $R_j:{\mathbb{R}}_{\ge 0}^n\to {\mathbb{R}}_{\ge 0}^{\nu }$. We assume that $R_j$ takes the form of Mass-Action kinetics: $R_j(z)={\prod }_{i=1}^n{k}_j{z}_j^{{\alpha }_{ij}}$, where $k_j$ is a kinetic constant. Letting $z(t)\in {\mathbb{R}}_{\ge 0}^n$ be the vector of species concentrations at time $t$, the associated ODE can be written as: $\dot{z}=\mathrm{\Gamma }R(z)$, where $R:=[R_1,..,R_{\nu }{]}^T$.

The goal of our models is to grasp the essence and the most relevant dynamics that give rise to the CFN with bistability. Unfortunately, it is difficult to decide on appropriate ranges of kinetic parameters for such coarse grained level models. One restriction that we have imposed on parameters is to require bistability, meaning (for both the single gene and the full model) that there should be two stable attractors (and a saddle node) for the underlying dynamics. One stable attractor corresponds to the somatic (silenced) state, and the other one corresponds to the active (pluripotent) state.

Full model description

We grouped our CRN for the full model into eight modules with the reactions numbered from R1 to R18 as shown in Table 1.

1. Oct4 promoter module (R1-R4): This module contains four reactions. In our modelling context, we assume that the Oct4 gene has two binding sites. The first binding site is for Oct4 protein itself [86,87], and the second binding site is for Nanog-TET heterodimer complex [42] as shown in the two-pathway activation scheme in Figure 4. Although we have shown NT binding to the first site and Oct4 binding to the second site in Figure 4, the cases of Oct4 binding to the first site and NT binding to the second site are the same for our modelling. Both cases are considered as a single-site occupancy and are denoted by $D_{10}^O$ in Table 1. As long as two sites are occupied, the corresponding gene will get activated. The parameter ${\alpha }_{NT}$ is the forward rate for all reactions that involve NT binding to a promoter, and $K_{NT}$ is the dissociation rate for all reversible reactions associated with NT binding to a promoter. The parameter ${\alpha }_O$ is the forward rate for all reactions that involve Oct4 binding to a promoter, and $K_O$ is the dissociation rate for all reversible reactions associated with Oct4 binding to a promoter.

2. TET promoter module (R5-R8): This module contains 4 reactions. Similar to the Oct4 promoter, the TET promoter region is also able to bind with the NT dimer and the Oct4 protein [31,36,42]. For the reaction rate, we adopt the same rate parameter for TET as for Oct4.

3. Nanog promoter module (R9): This module describes the reversible reaction of Oct4 protein binding to the Nanog promoter [87].

4. Nanog-TET module (R10): This module describes the fact that the TET-mediated DNA demethylation cycle is guided by the Nanog protein in the form of an NT complex [53]. The parameters $a_1,K_d$ are the association and dissociation rates of the heterodimerization of Nanong and TET.

5. Protein decay and production module (R11-R13): In this module, we assume Nanog, TET, and Oct4 all have the same degradation rate $\delta$. When each gene promoter has been fully occupied by TFs, the corresponding gene will be activated. The protein production rates for each activated promoter $D_1^N$, $D_{11}^T$, $D_{11}^O$ are given by the reaction rates ${\alpha }_N$, ${\alpha }_T$ and ${\alpha }_O$.

6. Protein basal production module (R14-R16): In this module, we assume that TET protein has a basal production rate $m_2$. The rates $m_1$ (basal rate for Nanog) and $m_3$ (basal rate for Oct4) are 0. We used the values of such basal protein production rates from [36].

7. Oct4 promoter Methylation module (R17): This module describes the de novo methylation of the Oct4 gene promoter region with a methylation rate $\gamma$. We assume that the Oct4 promoter can be methylated only if its promoter is unbound. This is supported by experimental observations that show that Oct4 protects its promoter from methylation when it is bound to it [88].

8. Nanog-guided TET demethylation module (R18): In this module, we assume that the triangle topology of DNA demethylation cycle has a slow dynamics on the 5hmC oxidized promoter state [31]. Using a quasi-steady state approximation, we reduced the demethylation cycle from the triangle topology to a two-state CRN with the effective demethylation rate given by $\theta$. (See Figure 5)

In our model, the DNA methylation rate refers to the effective rate combining de novo methylation rate, mainly contributed by DNMT3a/b, and maintenance methylation (upon replication) contributed by both DNMT1 and DNMT3a/b [89,90]. DNMT3 is a family of DNA methyltransferases that can methylate hemimethylated and unmethylated CpGs at the same rate, while DNMT1 has the tendency to methylate hemimethylated DNA after replication [90]. In our model of the DNA demethylation cycle, the 5mC state is considered as a fully methylated state, and our model is not affected by the gene methylation state being hemi-methylated state for a transient short period.

In order to understand the basic effects of DNAm on silencing, and to maintain the intelligibility and tractability of the model, we assume that the methylation status of a promoter is binary: it is either methylated or not. As we have illustrated in Figure 1, this is a simplified representation of the real picture where methylation is a complex process. Recent results suggest that the methylation of individual cytosines can be important in certain cases [91]. Furthermore, the methylation of Oct4 starts initially from two sites: the proximal enhancer and the distal promoter, spreading later to neighbouring regions [45]. In addition, the average promoter methylation does not usually reach to either 0% or 100% [92]. However, these intricate details are not relevant to our reduced model. For instance, if 80% methylation of Oct4 promoter is sufficient for silencing, and the cell is somatic at phenotypic level, we assume that an average DNAm level of 80% to be functionally equivalent to a fully methylated promoter.

ODE system for the full CRN model

The corresponding ODE model for the PSCC (Table 1) is given below:

$\frac{dN(t)}{dt}=-aNT+a{k}_d[NT]+{\alpha }_N{D}_1^N-\delta N+m_1$

$\frac{dT(t)}{dt}=-aNT+a{k}_d[NT]+{\alpha }_T{D}_{11}^T-\delta T+m_2$

$\frac{d[NT](t)}{dt}=aNT-a{k}_d[NT]-a_{nt}[NT]D_{00}^O+K_{nt}{a}_{nt}{D}_{10}^O-a_{nt}[NT]D_{01}^O+K_{nt}{a}_{nt}{D}_{11}^O$

$-a_{nt}[NT]D_{00}^T+K_{nt}{a}_{nt}{D}_{10}^T-a_{nt}[NT]D_{01}^T+K_{nt}{a}_{nt}{D}_{11}^T$

$\frac{d{D}_{00}^O(t)}{dt}=-a_{nt}[NT]D_{00}^O+K_{nt}{a}_{nt}{D}_{10}^O-a_{O}O{D}_{00}^O+K_O{a}_O{D}_{01}^O-\gamma D_{00}^O+\theta D_m[NT]$

$\frac{d{D}_{10}^O(t)}{dt}=a_{nt}[NT]D_{00}^O-K_{nt}{a}_{nt}{D}_{10}^O-a_{O}O{D}_{10}^O+K_O{a}_O{D}_{11}^O$

$\frac{d{D}_{01}^O(t)}{dt}=-a_{nt}[NT]D_{01}^O+K_{nt}{a}_{nt}{D}_{11}^O+a_{O}O{D}_{00}^O-K_O{a}_O{D}_{01}^O$

$\frac{d{D}_{11}^O(t)}{dt}=a_{nt}[NT]D_{01}^O-K_{nt}{a}_{nt}{D}_{11}^O+a_{O}O{D}_{10}^O-K_O{a}_O{D}_{11}^O$

$\frac{dO(t)}{dt}=-a_{O}O{D}_{00}^O+K_O{a}_O{D}_{01}^O-a_{O}O{D}_{10}^O+K_O{a}_O{D}_{11}^O-a_{O}O{D}_{00}^T+K_O{a}_O{D}_{01}^T$

$-a_{O}O{D}_{10}^T+K_O{a}_O{D}_{11}^T-a_{O}O{D}_0^N+K_O{a}_O{D}_1^N+{\alpha }_O{D}_{11}^O-\delta O+m_3$

$\frac{d{D}_{00}^T(t)}{dt}=-a_{nt}[NT]D_{00}^T+K_{nt}{a}_{nt}{D}_{10}^T-a_{O}O{D}_{00}^T+K_O{a}_O{D}_{01}^T$

$\frac{d{D}_{10}^T(t)}{dt}=a_{nt}[NT]D_{00}^T-K_{nt}{a}_{nt}{D}_{10}^T-a_{O}O{D}_{10}^T+K_O{a}_O{D}_{11}^T$

$\frac{d{D}_{01}^T(t)}{dt}=-a_{nt}[NT]D_{01}^T+K_{nt}{a}_{nt}{D}_{11}^T+a_{O}O{D}_{00}^T-K_O{a}_O{D}_{01}^T$

$\frac{d{D}_{11}^T(t)}{dt}=a_{nt}[NT]D_{01}^T-K_{nt}{a}_{nt}{D}_{11}^T+a_{O}O{D}_{10}^T-K_O{a}_O{D}_{11}^T$

$\frac{d{D}_0^N(t)}{dt}=-a_{O}O{D}_0^N+K_O{a}_O{D}_1^N$

$\frac{d{D}_1^N(t)}{dt}=a_{O}O{D}_0^N-K_O{a}_O{D}_1^N$

$\frac{d{D}_m(t)}{dt}=\gamma D_{00}^O-\theta D_m[NT]$

Steady state calculation

For the models defined above (the single gene model is (2)-(5), and the full model given by the CRN in Table 1), we first calculate the number of steady states that the model admits given a certain set of parameters. We use a Homotopy Continuation (HC) method to find all possible steady state solutions for our dynamical system models. The Global HC method is currently the best candidate to find all steady states solution for the system [93].

Sampling method

For the single gene model

We calculate the BoAp with various parameter sets that give rise to three steady states. Then, we calculate the BoAp against $\gamma$ and $\theta$ individually. Though the model contains four state variables $O$, $D_i$, $D_a$, $D_m$, it is essentially a three-dimensional model due to the conservation of the total concentration of promoters. For simplicity, we normalize total concentration of promoters to an arbitrary constant, say 5.

In order to calculate the BoAp, we need to specify a predefined volume of interest. We have defined the hyper-cube $[0,o]\times [0,d_1]\times [0,d_2]$ where $d_1,d_2,o$ are 1.5 times the corresponding high steady state. Because the conservation law puts a constraint on the total concentration of promoters, for each $D_i$, $D_a$ and $D_m$, their concentration should be all less than the total (in the simulation we choose the concentration total = 5). We sample the state space (with a given grid size) and count the number of points in the defined grided space that converges to the low steady state. We then plot BoAp with respect to $\gamma$ and $\theta$ for various numbers of parameter sets and show that they display a consistent trend.

For the full model

Given the CRN from Table 1, we first calculate steady state solutions that show bistability. Then, we use quasi-steady steady approximation to reduce the full model from 17 dimensions to 4 dimensions and simulate the BoAp in 4D space with a predefined hypercube similar to the previous paragraph.

Parameters

The ‘free’ parameters below are constrained to the set which can give rise to bistability. The following ‘free’ parameters can be thought of as normalized constants that can be re-scaled to their experimental values when applied to a biological context.

In the single gene model, the kinetic rate parameters are: $\beta$, ${\beta }_{-}$, $k$, $\delta$, $\gamma$, $\theta$. In the simulation of BoAp vs. $\gamma$ plot, we choose the parameter ($\beta$, ${\beta }_{-}$, $k$, $\delta$, $\theta$) values ranging from 0 to 5 with a step size of 0.5. For the simulation of BoAp vs $\theta$ plot, the same range is used for $\beta$, ${\beta }_{-}$, $k$, $\delta$, $\gamma$.

The range can be chosen to be $[0,a]$, for any $a>0$. This is since the ODE (2)-(5) follows Mass-action kinetics and all the rates can be scaled to have values from 0 to $a$ for the purpose of BOAp investigation.

In the full model, the kinetic rate parameters that define the system are: $K_O$, $K_{nt}$, $K_d$, $a$, $a_{nt}$, $a_O$, ${\alpha }_T$, ${\alpha }_O$, ${\alpha }_N$, $\delta$, $\gamma$, $\theta$, $m_1$, $m_2$ and $m_3$ as shown in Table 1. Although the values for the three dissociation parameters $K_O$, $K_{nt}$, $K_d$ were given in [36], in our full model simulation results, we are varying $K_O$, $K_{nt}$, $K_d$ from 0 to 0.5 with 0.1 step size, and $\delta$ from 1 to 5 with 0.5 step size, and the parameter $\theta$ from 1 to 20 with a step size of 2. We set protein production rates ${\alpha }_T$, ${\alpha }_O$ and ${\alpha }_N$ for Nanog, TET and Oct4 to be 1, the parameter $a=1$, and finally the parameters $a_{nt}$ and $a_O$ are set to be 100. For those parameters that are set to be constants, we did extensive simulations that have shown that the BoAp is insensitive to large changes in their values. Hence, our presented simulations only contain parameter sets that have noticeable influence on the BoAp.

Methylation association rate ($R$)

In the single gene model, we choose the methylation association rate ($R$) to be from 0.2 to 1.3 so that it covers both the “$R<1$“ regime and the “$R>1$“ regime of the dynamics with a particular parameter set. In the full model, we choose the methylation association rate ($R$) to span across several orders of magnitude, as possible, ranging from 0.01 to 30. The values of $R$ for each model depend on the parameter set we choose, but the qualitative trend for all the parameter sets we tested is the same and shows a biphasic trend. For the full model in Figure 8 we have chosen the following parameters: $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\}$, and with $\gamma$ spanning the range $[{10}^{-2},{10}^2]$. For those $R$ values outside the range, the system does not exhibit bistability for this particular parameter set.

Code availability

The programming code that was used to analyze the raw data that supports the findings of this study is available from the corresponding author upon request.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Science Foundation [1849588]; National Science Foundation [1716623].

Notes

1. We assume that transcription is fast enough that no explicit consideration of mRNA abundances is required, but it would be straightforward to add mRNA intermediates to the model.