Mathematical analysis of the limiting behaviors of a chromatin modification circuit

Abstract

In the last decade, the interactions among histone modifications and DNA methylation and their effect on the DNA structure, i.e., chromatin state, have been identified as key mediators for the maintenance of cell identity, defined as epigenetic cell memory. In this paper, we determine how the positive feedback loops generated by the auto- and cross-catalysis among repressive modifications affect the temporal duration of the cell identity. To this end, we conduct a stochastic analysis of a recently published chromatin modification circuit considering two limiting behaviors: fast erasure rate of repressive histone modifications or fast erasure rate of DNA methylation. In order to perform this mathematical analysis, we first show that the deterministic model of the system is a singular singularly perturbed (SSP) system and use a model reduction approach for SSP systems to obtain a reduced one-dimensional model. We thus analytically evaluate the reduced system’s stationary probability distribution and the mean switching time between active and repressed chromatin states. We then add a computational study of the original reaction model to validate and extend the analytical findings. Our results show that the absence of DNA methylation reduces the bias of the system’s stationary probability distribution toward the repressed chromatin state and the temporal duration of this state’s memory. In the absence of repressive histone modifications, we also observe that the time needed to reactivate a repressed gene with an activating input is less stochastic, suggesting that repressive histone modifications specifically contribute to the highly variable latency of state reactivation.

1 Introduction

Multicellular organisms are composed of cells with different phenotypes, even if all cells share the same genetic sequence, and this phenotypic distinction is maintained despite the unavoidable presence of noise. This property is known as epigenetic cell memory (ECM). For instance, ECM allows human differentiated cells to have different identities, even if they share the same genetic sequence, and to maintain these identities across cell division. In the past years, several studies have shown how the structure of the DNA, defined as chromatin state and determined by histone modifications and DNA methylation, affects gene expression and then has a critical role in ECM [1, 2].

This is the reason why several models describing the chromatin dynamics have been developed and analyzed. While some of them include either histone modifications or DNA methylation but not both [3, 4], and others are suitable only for computational analysis [5,6,7], a chemical reaction model including both DNA methylation and histone modifications has only recently appeared [8]. The circuit comprises positive feedback loops generated by the cooperation and competition among chromatin modifications.

In this paper, we focus on determining the specific contributions of histone modifications and DNA methylation to the features of the stationary probability distribution and temporal duration of cell memory. To this end, we perform a mathematical analysis of this model using the theory of singular singularly perturbed systems [9]. More precisely, we exploit this theory to reduce the chromatin modification model to a one-dimensional system, which we use to create a one-dimensional Markov chain suitable for analytical study. Then, we consider two limiting cases. In the first case, we consider the limit in which the erasure rate of DNA methylation becomes much larger than the erasure rate of the other chromatin modifications, obtaining a reduced system in which DNA methylation is absent. In the second case, we consider the limit in which the erasure rate of the repressive histone modifications becomes much larger than the erasure rate of the other chromatin modifications, thereby obtaining a reduced system with no repressive histone modifications. In all cases considered, the expression for the stationary distribution is obtained by exploiting detailed balance [10], while first step analysis [11] is applied to analytically evaluate the temporal duration of memory. Finally, we validate and extend the analytical results with computational simulations of the original reaction model using Gillespie’s stochastic simulation algorithm (SSA) [12].

Fig. 1

Reactions and diagram of the chromatin modification circuit. a List of the reactions associated with the full chromatin modification circuit. The numbers associated with the reactions are described in the main text. The boxes enclose reactions associated with activating histone marks (blue), repressive histone marks (pink), and DNA methylation (brown). Dark shades are associated with the establishment and light shades are associated with erasure. b Full chromatin modification circuit diagram. c Simplified circuit diagram in which DNA methylation is absent. d Simplified circuit diagram in which repressive histone modification is absent. In b–d, each arrow corresponds to reactions in a associated with the same number. More precisely, solid arrows represent the establishment and erasure reactions of chromatin modifications, while dashed arrows represent the increase of the establishment and erasure reaction rate due to the presence of another species (Color figure online)

2 Models

The chromatin modification circuit model analyzed in this work is developed in [8]. It includes H3K9 methylation (H3K9me3) and DNA methylation (CpGme), associated with repressed chromatin state [13], H3K4 methylation/acetylation (H3K4me3/ac), associated with active chromatin state ([1], Chapter 3 and [14]), and their known interactions.

The basic unit of the model is D, that is, the nucleosome with DNA wrapped around it. Then, we have $\mathrm{D^A}$, that is, nucleosome with H3K4me3/ac, $\mathrm{D^R_2}$, that is, nucleosome with H3K9me3, $\mathrm{D^R_1}$, that is, nucleosome with CpGme, and $\mathrm{D^R_{12}}$, that is, nucleosome with both H3K9me3 and CpGme. In terms of key molecular interactions considered in the model, all the modifications can be de novo established (reactions , , and ). Then, the read-write mechanism, in which histone modifications recruit marks of the same kind to nearby nucleosomes, generates auto-catalytic loops (reactions , ). Similarly, the cooperation between DNA methylation and repressive histone modification, through which each mark enhances the creation of the other, generates cross-catalytic loops (reactions ). Finally, basal erasure and recruited erasure, wherein activating marks recruit repressive mark’s eraser enzymes and vice versa, are represented by reactions , , and , , , respectively. In this reaction model, it is introduced the assumption that the rate of the establishment, auto- and cross-catalysis, and erasure of H3K9me3 (DNA methylation) does not change if the other repressive chromatin modification is present on the same nucleosome. All the reactions are listed in Fig. 1a, and the diagrams of the chromatin modification models are represented in Fig. 1b–d. More precisely, Fig. 1b shows the full chromatin modification circuit, Fig. 1c shows the simplified chromatin modification circuit including only activating and repressive modifications, and Fig. 1d shows the simplified chromatin modification circuit that only includes activating histone modifications and DNA methylation.

Now, let us introduce the ordinary differential equation (ODE) model associated with the full chromatin modification circuit. More precisely, by letting $n_{D^A}$, $n_{D^R_1}$, $n_{D^R_2}$, $n_{D^R_{12}}$ and $n_{D}$ represent the number of $\mathrm{D^A, D^R_1, D^R_2, D^R_{12}, D}$, we introduce the ODE model in terms of the fractions ${\bar{D}}^A=n_{D^A}/{\text {D}}_{\textrm{tot}}$, ${\bar{D}}^R_1=n_{D^R_1}/{\text {D}}_{\textrm{tot}}$, ${\bar{D}}^R_2=n_{D^R_2}/{\text {D}}_{\textrm{tot}}$, ${\bar{D}}^R_{12}=n_{D^R_{12}}/{\text {D}}_{\textrm{tot}}$ and ${\bar{D}}=n_{D}/{\text {D}}_{{tot}}$, with ${\text {D}}_{\textrm{tot}}$ the total number of modifiable units, that is the total number of nucleosomes within the gene of interest. This can be done by assuming ${\text {D}}_{\textrm{tot}}$ sufficiently large, such that $n_{D^A}$, $n_{D^R_1}$, $n_{D^R_2}$, $n_{D^R_{12}}$ and $n_{D}$ can be considered real-valued. Now, let us introduce $D_{{tot}}={\text {D}}_{\textrm{tot}}/\Omega$, with $\Omega$ the reaction volume, and let us define the normalized inputs as ${\bar{u}}^R_1=u^R_{10}+u^R_1$, ${\bar{u}}^R_2=u^R_{20}+u^R_2$ and ${\bar{u}}^A=u^A_0+u^A$, with

$$\begin{aligned} u^R_{i0}=\frac{k^i_{W0}}{k^A_MD_{{tot}}}, u^R_{i}=\frac{k^i_{W}}{k^A_MD_{{tot}}},\hbox {for\;}i=1,2,\;\;\hbox {and\;} u^A_{0}=\frac{k^A_{W0}}{k^A_MD_{{tot}}}, u^A=\frac{k^A_{W}}{k^A_MD_{{tot}}}.\nonumber \\ \end{aligned}$$

(1)

We consider ${\bar{u}}^R_1,{\bar{u}}^R_2$ and ${\bar{u}}^A$ as inputs of our dynamical system because they can be modulated by transcription factors external to the chromatin modification circuit [8]. Now, let us define the parameters $\alpha = k_M/k^A_M$, ${\bar{\alpha }}= {\bar{k}}_M/k^A_M$ and $\alpha '= k^{'}_M/k^A_M$: the first one represents the dimensionless rate constant of the auto-catalytic loops, while the last two represent the dimensionless rate constants of the cross-catalytic loops. In our analysis, without loss of generality, we introduce the simplifying assumption that these three parameters have the same order. Finally, let us also define

$$\begin{aligned} \varepsilon =\frac{\delta +\bar{k}^A_E}{k^A_M D_{{tot}}},\;\;\varepsilon '=\frac{k^A_E}{k^A_M},\;\;\mu =\frac{k^R_E}{k^A_E},\;\;\mu '=\frac{k^{'*}_T}{k^A_E}, \end{aligned}$$

(2)

with $b=O(1)$ such that $(\delta +{\bar{k}}^R_E)/({\delta + {\bar{k}}^A_E}) = b\mu$ and $\beta =O(1)$ such that $({\delta ^{'}+k^{'}_T})/({\delta + {\bar{k}}^A_E}) = \beta \mu '$. More precisely, $\mu$ represents the ratio between the erasure rates of repressive histone modifications and activating histone modifications and $\mu '$ represents the ratio between the erasure rates of DNA methylation and activating histone modifications. Furthermore, based on the previous definitions, we have that $(\delta +{\bar{k}}_E^R)/(k^A_M D_{{tot}})=b\varepsilon \mu$ and $(\delta '+k_T')/(k^A_M D_{{tot}})=\beta \varepsilon \mu '$. This implies that the dimensionless parameter $\varepsilon$ scales the ratio between the rate of the basal erasure and the one of the auto-/cross-catalytic loop of each mark. Finally, given that $k_E^R/k_M^A=\mu \varepsilon '$ and $k_T'^*/k_M^A=\mu ' \varepsilon '$, the dimensionless parameter $\varepsilon '$ scales the ratio between the rate of the recruited erasure and the one of the auto-/cross-catalytic loop of each mark. We collect the definitions and interpretations of these parameters in Table 1. Now defining the normalized time $\tau =tk^A_MD_{{tot}}$, the ODEs associated with the full chromatin modification circuit are

Table 1 Definitions and interpretations of $\varepsilon$, $\varepsilon '$, $\mu$, and $\mu '$

$$\begin{aligned} \frac{\hbox {d}\bar{D}^R_1}{\hbox {d} \tau }&=({\bar{u}}^R_1+\alpha '(\bar{D}^R_2+\bar{D}^R_{12})){\bar{D}}+\mu (b \varepsilon +\varepsilon '\bar{D}^A)\bar{D}^R_{12}\nonumber \\&\quad -(u^R_{20}+\alpha (\bar{D}^R_2+\bar{D}^R_{12})+\bar{\alpha }(\bar{D}^R_1+\bar{D}^R_{12})+\mu '(\beta \varepsilon +\varepsilon '\bar{D}^A))\bar{D}^R_1\nonumber \\ \frac{\hbox {d}\bar{D}^R_2}{\hbox {d} \tau }&=({\bar{u}}^R_2+\alpha (\bar{D}^R_2+\bar{D}^R_{12})+\bar{\alpha }(\bar{D}^R_1+\bar{D}^R_{12})){\bar{D}}+\mu '(\beta \varepsilon +\varepsilon '\bar{D}^A)\bar{D}^R_{12}\nonumber \\&\quad -(u^R_{10}+\alpha '(\bar{D}^R_2+\bar{D}^R_{12})+\mu (b \varepsilon +\varepsilon '\bar{D}^A))\bar{D}^R_2\nonumber \\ \frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d} \tau }&=(u^R_{10}+\alpha '(\bar{D}^R_2+\bar{D}^R_{12}))\bar{D}^R_2+(u^R_{20}+\alpha (\bar{D}^R_2+\bar{D}^R_{12})+\bar{\alpha }(\bar{D}^R_1+\bar{D}^R_{12}))\bar{D}^R_1\nonumber \\&\quad -(\mu {'}(\beta \varepsilon +\varepsilon '\bar{D}^A)+\mu (b\varepsilon +\varepsilon '\bar{D}^A))\bar{D}^R_{12}\nonumber \\ \frac{\hbox {d}\bar{D}}{\hbox {d}\tau }&=\mu '(\beta \varepsilon + \varepsilon ' \bar{D^A}){\bar{D}}^R_1+\mu (b \varepsilon + \varepsilon '\bar{D^A}){\bar{D}}^R_2\nonumber \\&\quad +(\varepsilon +\varepsilon '({\bar{D}}^R_1+{\bar{D}}^R_{12})+\varepsilon '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})\nonumber \\&\quad +{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+{\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A){\bar{D}}\nonumber \\ \frac{\hbox {d} \bar{D}^A}{\hbox {d} \tau }&= ({\bar{u}}^A+\bar{D}^A){\bar{D}}-(\varepsilon +\varepsilon '(\bar{D}^R_2+\bar{D}^R_{12})+\varepsilon '(\bar{D}^R_1+\bar{D}^R_{12}))\bar{D}^A, \end{aligned}$$

(3)

with ${\bar{D}}+{\bar{D}}^A+{\bar{D}}^R_1 +{\bar{D}}^R_2 + {\bar{D}}^R_{12} = 1$ as constraint.

We next write the system in singular perturbation form by assuming that the rates associated with the auto- and cross-catalysis are much faster than the rate of the erasure processes. This assumption is in agreement with empirical findings suggesting that the natural erasure of chromatin marks is a slow process [15]. We thus let $\varepsilon =c\varepsilon '$, with $c=O(1)$. Then, introducing the new time variable ${\bar{\tau }} =\tau \varepsilon '$ in the ODEs (3), the system of equations can be rewritten as follows:

$$\begin{aligned} \varepsilon '\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}&= ({\bar{u}}^A+{\bar{D}}^A){\bar{D}}-\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\ \nonumber \varepsilon '\frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}&=(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^R_2+(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}^R_1\nonumber \\&\quad -\varepsilon '(\mu (b c+ \bar{D^A})+\mu '(\beta c+ \bar{D^A})){\bar{D}}^R_{12}\nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}^R_1}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\nonumber \\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}^R_1-(\varepsilon '\mu '(\beta c+ \bar{D^A})){\bar{D}}^R_1 \nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\nonumber \\&\quad -(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+\varepsilon '\mu (b c+ \bar{D^A})){\bar{D}}^R_2\nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}}{\hbox {d}\bar{\tau }}&=\varepsilon '(\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_1+\mu (b c+ \bar{D^A}){\bar{D}}^R_2)\nonumber \\&\quad +\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A-({\bar{u}}^A+{\bar{D}}^A){\bar{D}}\nonumber \\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+{\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}. \end{aligned}$$

(4)

Here, $\varepsilon '$ is the small parameter that we exploit in the model reductions performed in Sect. 4. Furthermore, when we investigate the limiting behavior of the system for large $\mu '$, we define the parameter $U':=\varepsilon '\mu '$ and assume that it is $\varepsilon '$-independent, so that the product $\varepsilon '\mu '$ does not vanish as $\varepsilon '$ approaches zero. Similarly, in order to study the behavior of the system for large $\mu$, we will define $U:=\varepsilon '\mu$ and assume that it is $\varepsilon '$-independent, so that $\varepsilon '\mu$ does not vanish as $\varepsilon '$ approaches zero. Each case will lead to a different reduced model, as shown in Sect. 4.

3 Singular singularly perturbed system and model reduction approach

In this section, we introduce the definition of singular singularly perturbed system and the model reduction approach developed in [9], which we will apply to the full chromatin modification circuit model (3). Let us first give the definition of integral manifold S provided in [16, 17]:

Definition 3.1

(Integral manifold) Given a general dynamical system $\frac{\hbox {d}x}{\hbox {d}t} = f(x,y,t)$ with $x\in \mathbb {R}^n$, $t\in \mathbb {R}$, let us define a smooth surface S in $\mathbb {R}^n \times \mathbb {R}$ as an integral manifold of the system if any trajectory (x(t), t) of the system, that has at least one point in common with S, lies entirely in S.

Definition 3.2

(Singular singularly perturbed system)

Assuming $x\in \mathbb {R}^m$ and $y_2\in \mathbb {R}^n$, let us introduce the system

$$\begin{aligned} \begin{aligned} \varepsilon '\dot{x} = f_1(x,y_2,t,\varepsilon ')\;\;\;\;\; \varepsilon '\dot{y}_2 = f_2(x,y_2,t,\varepsilon '), \end{aligned} \end{aligned}$$

(5)

with functions $f_1$ and $f_2$ sufficiently smooth, and let us define the matrix

$$\begin{aligned} A(x,y_2,t,\varepsilon ')= \begin{pmatrix} \frac{\partial f_1}{\partial x}&{}\frac{\partial f_1}{\partial y_2}\\ \frac{\partial f_2}{\partial x}&{}\frac{\partial f_2}{\partial y_2}\\ \end{pmatrix}= \begin{pmatrix} f_{1x}&{}f_{1y_2}\\ f_{2x}&{}f_{2y_2}\\ \end{pmatrix}. \end{aligned}$$

(6)

If $A(x,y_2,t,0)$ is singular on some subspace of $\mathbb {R}^m\times \mathbb {R}^n \times \mathbb {R}$, then we define the system in (5) as a singular singularly perturbed system [9].

Now, let us consider the following conditions [9]:

• C1: $f_2(x,y_2,t,0)=0$ has a smooth isolated root $y_2=\phi (x,t)$ with $t\in \mathbb {R}$ and $x\in \mathbb {R}^m$;

• C2: the matrix A, defined in (6), with $y_2$ = $\phi (x,t)$ and $\varepsilon '=0$ has a kernel of dimension m and m corresponding linearly independent eigenvectors, and the matrix

  $$\begin{aligned} B(x,\phi (x,t),t,0)=\frac{\partial f_2(x, \phi (x,t), t, 0)}{\partial y_2} \end{aligned}$$

  (7)

  has n eigenvalues $\lambda _i(x,t):Re(\lambda _i)\le -2\theta$, with $\theta >0$;

• C3: defining the domain $\mathcal {X}$ as $\mathcal {X}=\{(x,y_2,t,\varepsilon ')|x \in$ $\mathbb {R}^m,||y_2 - \phi (x,t)||\le \rho ,t\in \mathbb {R}, 0\le \varepsilon '\le \varepsilon '_0\}$, the functions $f_1$, $f_2$ and the matrix A are continuously differentiable $(k+2)$ times in $\mathcal {X}$, with $k\ge 0$ for some positive $\varepsilon '_0$ and $\rho$.

Now, let us introduce the new variables $y_2=y_1+\phi (x,t)$ in system (5), that can then be rewritten as follows:

$$\begin{aligned} \begin{aligned} \varepsilon '\dot{x}= & {} C(x,t)y_1+F_1(x,y_1,t)+\varepsilon 'X(x,y_1,t,\varepsilon ')\\ \varepsilon '\dot{y}_1= & {} B(x,t)y_1+F_2(x,y_1,t)+\varepsilon 'Y(x,y_1,t,\varepsilon '), \end{aligned} \end{aligned}$$

(8)

with

$$\begin{aligned} C(x,t)&=f_{1y_2}(x,\phi (x,t),t,0),\nonumber \\ B(x,t)&=f_{2y_2}(x,\phi (x,t),t,0),\nonumber \\ F_1(x,y_1,t)&=f_1(x,y_1+\phi (x,t),t,0)-C(x,t)y_1,\nonumber \\ F_2(x,y_2,t)&=f_2(x,y_1+\phi (x,t),t,0)-B(x,t)y_1,\nonumber \\ \varepsilon 'X(x,y_1,t,\varepsilon ')&=f_1(x,y_1+\phi (x,t),t,\varepsilon ')-f_1(x,y_1+\phi (x,t),t,0), \nonumber \\ \varepsilon 'Y(x,y_1,t,\varepsilon ')&=f_2(x,y_1+\phi (x,t),t,\varepsilon ')-f_2(x,y_1+\phi (x,t),t,0), \end{aligned}$$

(9)

in which $F_1$ and $F_2$ are such that $||F_1(x,y_1,t)||=O(||y_1||^2)$, $||F_2(x,y_1,t)||=O(||y_1||^2)$, and $\varepsilon ^{'-1}F_1(x,\varepsilon 'y,t)$ and $\varepsilon ^{'-1}F_2(x,\varepsilon 'y,t)$, with $y=y_1/\varepsilon$, are continuous in $\mathcal {X}$, with $\mathcal {X}$ defined in C3 [9]. Let us then consider the following theorem and remark:

Theorem 3.1

(Theorem 7.1 from [9]) If conditions C1–C3 are verified, then there exists an $\varepsilon '_1$, $0<\varepsilon '_1<\varepsilon '_0$, such that, for any $\varepsilon '\in (0,\varepsilon '_1)$, system (8) has a unique integral manifold, $y_1=\varepsilon 'h(x,t,\varepsilon ')$, that is exponentially attractive. The motion along this integral manifold is described by the following equation:

$$\begin{aligned} \begin{aligned} \dot{{\bar{x}}} = X_1({\bar{x}},t,\varepsilon '),\\ \end{aligned} \end{aligned}$$

(10)

in which $X_1({\bar{x}},t,\varepsilon ')=C({\bar{x}},t)h({\bar{x}},t,\varepsilon ')+X({\bar{x}},\varepsilon 'h,t,\varepsilon ')+\varepsilon ^{'-1}F_1({\bar{x}},\varepsilon 'h,t)$, with the function $h(x,t,\varepsilon ')$ continuously differentiable $k-$times with respect to x and t.

Remark 3.1

Given the fact that the integral manifold is exponentially attractive for a sufficiently small $\varepsilon '$, then, for any solution $(x(t),y_1(t))$ of (8) with initial conditions $(x(t_0),y_1(t_0)=(x^0,y^0_1)$ such that $|y^0_1-\varepsilon 'h(x^0,t_0,\varepsilon ')|$ is sufficiently small, we have a solution of the reduced system (10) such that

$$\begin{aligned} \begin{aligned} x(t)=\bar{x}(t)+\zeta _1(t),\;\;\;\;y_1(t)=\varepsilon 'h({\bar{x}}(t),t,\varepsilon ')+\zeta _2(t), \end{aligned} \end{aligned}$$

with $\zeta _i(t)=O(e^{-(\theta /\varepsilon ')(t-t_{0})})$, $i=1,2$, and $t\ge t_0$ ([9, 18, 19] Chapter 6). This implies that the behavior of the original system’s trajectories near the integral manifold can be determined by studying the reduced system’s trajectories (10).

Furthermore, as described in [9, 17], we can obtain $h(x,t,\varepsilon ')$ by introducing the change of variable $y=y_1/\varepsilon '$ in (8) and rewriting it in standard singular perturbation form:

$$\begin{aligned} \begin{aligned} \dot{x} = {\tilde{X}}(x,y,t,\varepsilon ')\;\;\;\;\; \varepsilon '\dot{y} = {\tilde{Y}}(x,y,t,\varepsilon '), \end{aligned} \end{aligned}$$

(11)

with $x\in \mathbb {R}^m,y\in \mathbb {R}^n,t\in \mathbb {R}$, ${\tilde{X}}(x,y,t,\varepsilon ') = C(x,t)y + \varepsilon ^{'-1}F_1(x,\varepsilon 'y,t) + X(x,\varepsilon 'y,t,\varepsilon ')$, ${\tilde{Y}}(x,y,t,\varepsilon ') = B(x,t)y + \varepsilon ^{'-1}F_2(x,\varepsilon 'y,t)+Y(x,\varepsilon 'y,t,\varepsilon ')$. Given that $F_i$, with $i=1,2$, are such that $||F_i(x,y_1,t)||=O(||y_1||^2)$ in $\mathcal {X}$, then $\varepsilon ^{'-1}F_i(x,\varepsilon 'y,t)$ are well defined as $\varepsilon '$ approaches zero [9]. Then, defining the smooth isolated root of ${\tilde{Y}}(x,y,t,0)=0$ as $y=h_0(x,t)$, it is possible to show that since conditions C1–C3 are verified, the eigenvalues $\lambda _i$ of the matrix $(\partial {\tilde{Y}}/\partial y)(x,h_0(x,t),t,0)$ satisfy the inequality $Re(\lambda _i)\le -2\theta$, with $\theta >0$. Then, the integral manifold $y=y_1/\varepsilon '=h(x,t,\varepsilon ')$ can be calculated as an asymptotic expansion in integer powers of $\varepsilon '$, $h(x,t,\varepsilon ')=h_0(x,t)+\varepsilon ' h_1(x,t)+\cdots +\varepsilon ^{'k} h_k(x,t)+\cdots$, whose coefficients are smooth function with bounded norm [17] and they can be found by substituting the expansion in the second equation of (11), obtaining [9]:

$$\begin{aligned} \begin{aligned} \varepsilon '\frac{\partial h}{\partial t}+\varepsilon '\frac{\partial h}{\partial x}{\tilde{X}}(x,h,t,\varepsilon ')={\tilde{Y}}(x,h,t,\varepsilon '). \end{aligned} \end{aligned}$$

(12)

4 Results

In this section, we obtain reduced versions of the full chromatin modification system (4) in the limit where $\varepsilon '$ approaches zero. We first show that, considering $\varepsilon '$ as small parameter, the ODE model (4) is a singular singularly perturbed system, and then we apply the approach introduced in Sect. 3 to obtain a reduced model. This allows us to obtain a one-dimensional reduced model with one-dimensional associated Markov chain, whose stochastic properties can be analytically determined, allowing a mechanistic understanding of how the parameters defined in Table 1 affect system behavior. Then, in order to validate the trends analytically determined, that rely on a deterministic quasi-steady state approximation [20], we conduct a computational study of the original reaction system and show that the analytically derived trends are mirrored by the original system.

In order to single out the contributions of repressive histone modification and DNA methylation to the system stochastic properties, we conduct the model reduction for two limiting cases. In particular, in order to single out the contribution of repressive histone modifications, we first introduce the $\varepsilon '$-independent parameter $U':=\varepsilon '\mu '$, so that $\varepsilon '\mu '$ does not vanish as $\varepsilon '$ approaches zero in the model reduction, and then we consider the limiting behavior as $U'\rightarrow \infty$. Similarly, in order to single out the contribution of DNA methylation, we repeat the model reduction introducing the $\varepsilon '$-independent parameter $U:=\varepsilon '\mu$, so that $\varepsilon '\mu$ does not vanish as $\varepsilon '$ approaches zero, and then we consider the limiting behavior as $U\rightarrow \infty$.

4.1 Behavior of the full chromatin modification circuit with $\varepsilon '$ as small parameter

Let us summarize the results of the model reduction in the following proposition:

Proposition 4.1

Let $\varepsilon =c\varepsilon '$, with $c=O(1)$, and let us consider the following system:

$$\begin{aligned} \frac{{\text {d}}\bar{D}^A}{{\text {d}}\tau }&=\left( \frac{(\mu (b \varepsilon + \varepsilon '{\bar{D}}^A)\mu '(\beta \varepsilon + \varepsilon '\bar{D^A})){\bar{K}}({\bar{u}}^A+{\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12}}\right) {\bar{D}}^R_{12}\nonumber \\&\quad -\left( \frac{(\varepsilon +2\varepsilon '{\bar{D}}^R_{12})({\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12})}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12}}\right) {\bar{D}}^A\nonumber \\ \frac{{\text {d}}\bar{D}^R_{12}}{{\text {d}}\tau }&=\left( \frac{(\varepsilon +2\varepsilon '{\bar{D}}^R_{12})({\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12})}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12}}\right) {\bar{D}}^A\\&\quad -\left( \frac{(\mu (b \varepsilon + \varepsilon '{\bar{D}}^A)\mu '(\beta \varepsilon + \varepsilon '\bar{D^A})){\bar{K}}({\bar{u}}^A+{\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha '){\bar{D}}^R_{12}}\right) {\bar{D}}^R_{12},\nonumber \end{aligned}$$

(13)

with ${\bar{D}}^A + {\bar{D}}^R_{12}=1$. Furthermore, let $({\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))=M({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau ))$ represent the system’s unique integral manifold. Then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau )$, ${\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))$ of (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_1(0),{\bar{D}}^R_2(0))-M({\bar{D}}^A(0),{\bar{D}}^R_{12}(0))|$ is sufficiently small, we have that, for $\varepsilon '$ sufficiently small, the solution of (13), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_{12}(\tau ))$, is such that

$$\begin{aligned} \begin{aligned} ({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau ))&=({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_{12}(\tau ))+\zeta _1(\tau ),\\ ({\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))&=M({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_{12}(\tau ))+\zeta _2(\tau ), \end{aligned} \end{aligned}$$

(14)

in which $\zeta _i(\tau )=O(e^{-(\theta /\varepsilon ')\tau })$, with $i=1,2$, $\theta >0$, and $\tau \ge 0$.

Proof

Let us consider the ODE model (4), and let us define x, $y_2$, $f_1$ and $f_2$ as follows:

$$\begin{aligned} x&= \begin{pmatrix} {\bar{D}}^A\\ {\bar{D}}^R_{12}\\ \end{pmatrix},y_2=\begin{pmatrix} {\bar{D}}^R_1\\ {\bar{D}}^R_2\\ {\bar{D}}\\ \end{pmatrix},f_1=\begin{pmatrix} f_{11}\\ f_{12}\\ \end{pmatrix},f_2=\begin{pmatrix} f_{21}\\ f_{22}\\ f_{23}\\ \end{pmatrix},\\ f_{11}&=({\bar{u}}^A+{\bar{D}}^A){\bar{D}}-\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A,\\ f_{12}&=(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^R_2+(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}^R_1\\&\quad -\varepsilon '(\mu (b c+ \bar{D^A})+\mu '(\beta c+ \bar{D^A})){\bar{D}}^R_{12},\\ f_{21}&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+\varepsilon '\mu '(\beta c+ \bar{D^A})){\bar{D}}^R_1,\\ f_{22}&=({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\\&\quad -((u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12}))+\varepsilon '\mu (b c+ \bar{D^A})){\bar{D}}^R_2,\\ f_{23}&=\varepsilon '(\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_1+\mu (b c+ \bar{D^A}){\bar{D}}^R_2\\&\quad +(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A) -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1\\&\quad +{\bar{D}}^R_{12})+{\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A){\bar{D}}. \end{aligned}$$

Now, it is possible to show that $\phi$, defined in C1, is equal to $\phi (x)=(0,0,0)$, and that matrix A, defined in (6), with ${\bar{D}} = {\bar{D}}^R_1= {\bar{D}}^R_2=0$ and $\varepsilon '=0$ can be written as follows:

$$\begin{aligned} A(x,y_2=\phi (x),t,0)= \left( \begin{array}{cc} 0_{2,2}&{}{\bar{A}}_{2,3}\\ 0_{3,2}&{}{\bar{A}}_{3,3}\\ \end{array}\right) \end{aligned}$$

(15)

with

$$\begin{aligned} {\bar{A}}_{2,3}&= \left( {\begin{matrix} 0&{}\qquad 0&{}\qquad ({\bar{u}}^A+{\bar{D}}^A)\\ (u^R_{20}+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12})&{}\qquad (u^R_{10}+\alpha '{\bar{D}}^R_{12})&{}\qquad 0\\ \end{matrix}}\right) ,\\ {\bar{A}}_{3,3}&= \left( {\begin{matrix} {\hat{A}}_{2,2}&{}\qquad {\hat{A}}_{2,1}\\ 0_{1,2}&{}\qquad {\hat{A}}_{1,1}\\ \end{matrix}}\right) ,\;\;\;\; {\hat{A}}_{2,1}= \left( {\begin{matrix} ({\bar{u}}^R_1+\alpha '{\bar{D}}^R_{12})\\ ({\bar{u}}^R_2+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12})\\ \end{matrix}}\right) ,\\ {\hat{A}}_{2,2}&= \left( {\begin{matrix} -(u^R_{20}+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12})&{}\qquad 0\\ 0&{}\qquad -(u^R_{10}+\alpha '{\bar{D}}^R_{12})\\ \end{matrix}}\right) ,\\ {\hat{A}}_{1,1}&= \left( {\begin{matrix}-({\bar{u}}^A+{\bar{D}}^A)-({\bar{u}}^R_1+{\bar{u}}^R_2+(\alpha +{\bar{\alpha }}+\alpha '){\bar{D}}^R_{12}) \end{matrix}}\right) . \end{aligned}$$

The matrix A in (15) is singular, and this implies that system (4) is singular singularly perturbed (Def. 3.2). More precisely, A has a twofold zero eigenvalue, with two associated linearly independent eigenvectors, and matrix $B={\bar{A}}_{3,3}$, with the definition of B given in (7). When no external inputs are applied ($u^A=u^R_1=u^R_2=0$ and then ${\bar{u}}^A=u^A_0$, ${\bar{u}}^R_1=u^R_{10}$, and ${\bar{u}}^R_2=u^R_{20}$), matrix B has three eigenvalues with negative real part if $u^R_{10}, u^R_{20}, u^A_0\ge l$ with $l>0$. This implies that we can apply Theorem 3.1 to reduce our system. To this end, let us first introduce the new variables ${\tilde{D}}={\bar{D}}/ \varepsilon '$, ${\tilde{D}}^R_1={\bar{D}}^R_1/ \varepsilon '$ and ${\tilde{D}}^R_2={\bar{D}}^R_2/ \varepsilon '$ in (4):

$$\begin{aligned} \varepsilon '\frac{\hbox {d}\tilde{D}^R_1}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_1+\alpha '(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})){\tilde{D}}+\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\nonumber \\&\quad -(u^R_{20}+\alpha (\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})+\varepsilon '\mu '(\beta c+ \bar{D^A})){\tilde{D}}^R_1\nonumber \\ \varepsilon '\frac{\hbox {d}\tilde{D}^R_2}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_2+\alpha (\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})){\tilde{D}}+\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\nonumber \\&\quad -((u^R_{10}+\alpha '(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12}))+\varepsilon '\mu (b c+ \bar{D^A})){\tilde{D}}^R_2\nonumber \\ \varepsilon '\frac{\hbox {d}\tilde{D}}{\hbox {d}\bar{\tau }}&=\mu '(\beta c+ \bar{D^A})\varepsilon ^{' }{\tilde{D}}^R_1+\mu (b c+ \bar{D^A})\varepsilon ^{' }{\tilde{D}}^R_2\nonumber \\&\quad +(c+(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})+(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\&\quad -({\bar{u}}^R_2+\alpha (\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})){\tilde{D}}\nonumber \\&\quad -({\bar{u}}^R_1+\alpha '(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A){\tilde{D}}\nonumber \\ \frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}&=(u^R_{10}+\alpha '(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})){\tilde{D}}^R_2\nonumber \\&\quad +(u^R_{20}+\alpha (\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})){\tilde{D}}^R_1\nonumber \\&\quad -(\mu (b c+ \bar{D^A})+\mu '(\beta c+ \bar{D^A})){\bar{D}}^R_{12}\nonumber \\ \frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}&= ({\bar{u}}^A+{\bar{D}}^A){\tilde{D}}-(c+(\varepsilon ^{' }{\tilde{D}}^R_1+{\bar{D}}^R_{12})+(\varepsilon ^{' }{\tilde{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A. \end{aligned}$$

(16)

Now, to determine the integral manifold $M({\bar{D}}^A,{\bar{D}}^R_{12})=({\bar{D}},{\bar{D}}^R_1,{\bar{D}}^R_2)$, let us find the expression for the asymptotic expansion of ${\tilde{D}}$, ${\tilde{D}}^R_1$ and ${\tilde{D}}^R_2$:

$$\begin{aligned} {\tilde{D}}&= h_0({\bar{D}}^A, {\bar{D}}^R_{12}, \varepsilon ')= h_{00}({\bar{D}}^A, {\bar{D}}^R_{12}) + \varepsilon 'h_{01}({\bar{D}}^A, {\bar{D}}^R_{12})+O(\varepsilon ^{'^2}),\nonumber \\ {\tilde{D}}^R_1&= h_1({\bar{D}}^A, {\bar{D}}^R_{12}, \varepsilon ')= h_{10}({\bar{D}}^A, {\bar{D}}^R_{12}) + \varepsilon 'h_{11}({\bar{D}}^A, {\bar{D}}^R_{12})+O(\varepsilon ^{'^2}), \nonumber \\ {\tilde{D}}^R_2&= h_2({\bar{D}}^A, {\bar{D}}^R_{12}, \varepsilon ')= h_{20}({\bar{D}}^A, {\bar{D}}^R_{12}) + \varepsilon 'h_{21}({\bar{D}}^A, {\bar{D}}^R_{12})+O(\varepsilon ^{'^2}). \end{aligned}$$

(17)

To this end, let us substitute (17) in the first three equations of (16) to obtain

$$\begin{aligned} \varepsilon '\frac{\hbox {d}h_1}{\hbox {d}\bar{\tau }}&=\varepsilon '\left( \frac{\partial h_1}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_1}{\partial {\bar{D}}^R_{12}}\frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}\right) \nonumber \\&=({\bar{u}}^R_1+\alpha '(\varepsilon ^{' }h_2+{\bar{D}}^R_{12}))h_0+\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\nonumber \\&\qquad -(u^R_{20}+\alpha (\varepsilon ^{' }h_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+{\bar{D}}^R_{12})+\varepsilon '\mu '(\beta c+ \bar{D^A}))h_1\nonumber \\ \varepsilon '\frac{\hbox {d}h_2}{\hbox {d}\bar{\tau }}&=\varepsilon '\left( \frac{\partial h_2}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_2}{\partial {\bar{D}}^R_{12}}\frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}\right) \nonumber \\&=({\bar{u}}^R_2+\alpha (\varepsilon ^{' }h_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+{\bar{D}}^R_{12}))h_0+\mu '(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\nonumber \\&\qquad -(u^R_{10}+\alpha '(\varepsilon ^{' }h_2+{\bar{D}}^R_{12})+\varepsilon '\mu (b c+ \bar{D^A}))h_2\nonumber \\ \varepsilon '\frac{\hbox {d}h_0}{\hbox {d}\bar{\tau }}&=\varepsilon '\left( \frac{\partial h_0}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_0}{\partial {\bar{D}}^R_{12}}\frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}\right) \nonumber \\&=\mu '(\beta c+ \bar{D^A})\varepsilon ^{' }h_1\nonumber \\&\qquad +\mu (b c+ \bar{D^A})\varepsilon ^{' }h_2+(c+(\varepsilon ^{' }h_1+{\bar{D}}^R_{12})+(\varepsilon ^{' }h_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\&\qquad -({\bar{u}}^R_2+\alpha (\varepsilon ^{' }h_2+{\bar{D}}^R_{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+{\bar{D}}^R_{12})+{\bar{u}}^R_1\nonumber \\&\qquad +\alpha '(\varepsilon ^{' }h_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A)h_0. \end{aligned}$$

(18)

Now, we can obtain $h_{i0}$ and $h_{i1}$, with $i=0,1,2$, by equating the terms of the right and left hand side of the equations above multiplied by the same power of $\varepsilon '$. More precisely, given that $\frac{\partial h_{i0}}{\partial {\bar{D}}^R_{12}}$ and $\frac{\partial h_{i0}}{\partial {\bar{D}}^A}$ are bounded $\forall i=0,1,2$, and then $\varepsilon '\frac{\partial h_{i0}}{\partial {\bar{D}}^R_{12}},\varepsilon '\frac{\partial h_{i0}}{\partial {\bar{D}}^A}\ll 1$ for sufficiently small values of $\varepsilon '$, we can write $h_{i0}$ and $h_{i1}$, with $i=0,1,2$, as follows:

$$\begin{aligned} h_{00}&=\frac{(c+2{\bar{D}}^R_{12}){\bar{D}}^A}{{\bar{u}}^R_2+\alpha {\bar{D}}^R_{12}+{\bar{\alpha }}{\bar{D}}^R_{12}+{\bar{u}}^R_1+\alpha '{\bar{D}}^R_{12}+{\bar{u}}^A+{\bar{D}}^A},\nonumber \\ h_{10}&=\frac{({\bar{u}}^R_1+\alpha '{\bar{D}}^R_{12})h_{00}+\mu \left( b c+ {\bar{D}}^A\right) {\bar{D}}^R_{12}}{u^R_{20}+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12}},\nonumber \\ h_{20}&=\frac{({\bar{u}}^R_2+\alpha {\bar{D}}^R_{12}+{\bar{\alpha }}{\bar{D}}^R_{12})h_{00}+\mu '\left( \beta c+ \bar{D^A}\right) {\bar{D}}^R_{12}}{u^R_{10}+\alpha '{\bar{D}}^R_{12}},\\ h_{01}&=\frac{(\mu (c b + {\bar{D}}^A)h_{20}+\mu '(\beta c + {\bar{D}}^A)h_{10})}{{\bar{u}}^R_2+\alpha {\bar{D}}^R_{12}+{\bar{\alpha }}{\bar{D}}^R_{12}+{\bar{u}}^R_1+\alpha '{\bar{D}}^R_{12}+{\bar{u}}^A+{\bar{D}}^A},\nonumber \\ h_{11}&=\frac{({\bar{u}}^R_1+\alpha '{\bar{D}}^R_{12})h_{01}-(\alpha h_{20}+{\bar{\alpha }} h_{10}+\mu '(\beta c + {\bar{D}}^A))h_{10}}{u^R_{20}+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12}},\nonumber \\ h_{21}&=\frac{({\bar{u}}^R_2+(\alpha +{\bar{\alpha }}){\bar{D}}^R_{12})h_{01}-(\alpha ' h^2_{20}+\mu (b c + {\bar{D}}^A))h_{20}}{u^R_{10}+\alpha '{\bar{D}}^R_{12}}.\nonumber \end{aligned}$$

(19)

Then, by introducing in the last two equations of (16) the asymptotic expansion of ${\tilde{D}}$, ${\tilde{D}}^R_1$ and ${\tilde{D}}^R_2$ (17) with the expressions for $h_{i0}$ and $h_{i1}$ provided in (19), we obtain the reduced system in (13), in which we have re-introduced the original time variable $\tau ={\bar{\tau }} / \varepsilon '$. The sum of the equations in (13) is equal to zero, implying that $\bar{D}^A+\bar{D}^R_{12}=$ constant. Since $\bar{D}^A+\bar{D}^R_{12}+\bar{D}+\bar{D}^R_1+\bar{D}^R_2=1$ and, for sufficiently small $\varepsilon '$, $\bar{D}=\varepsilon ' {\tilde{D}} \approx 0$, $\bar{D}^R_1=\varepsilon ' \tilde{D}^R_1 \approx 0$, and $\bar{D}^R_2=\varepsilon ' \tilde{D}^R_2 \approx 0$, then ${\bar{D}}^A + {\bar{D}}^R_{12}$ can be approximately set equal to 1 for sufficiently small $\varepsilon '$. Furthermore, given that the integral manifold obtained, $M({\bar{D}}^A,{\bar{D}}^R_{12})=(\bar{D},{\bar{D}}^R_1,{\bar{D}}^R_2)$, is exponentially attractive for a sufficiently small $\varepsilon '$ (see Theorem 3.1), then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau ),{\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))$ of the original system (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_1(0),{\bar{D}}^R_2(0))-M({\bar{D}}^A(0),{\bar{D}}^R_{12}(0))|$ is sufficiently small, we have a solution of the reduced system (13), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_{12}(\tau ))$, that satisfies (14) (see Remark 3.1). $\square$

Now, if we multiply both sides of the ODEs in (13) by $D_\mathrm{{tot}}(k^A_ED_{{tot}})$ and introduce ${\bar{k}}^A_W=k^A_{W0}+k^A_W$, ${\bar{k}}^1_W=k^1_{W0}+k^1_W$, and ${\bar{k}}^2_W=k^2_{W0}+k^2_W$, we can rewrite system (13) as follows:

$$\begin{aligned} \dot{D}^A&=\left( \frac{(\delta +{\bar{k}}^R_E+k^R_ED^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A){\bar{K}}_{dim}({\bar{k}}^A_W+k^A_MD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}\right) D^R_{12}\nonumber \\&\quad -\left( \frac{(\delta +{\bar{k}}^A_E+2k^A_E D^R_{12})({\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12})}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}\right) D^A\nonumber \\ \dot{D}^R_{12}&=\left( \frac{(\delta +{\bar{k}}^A_E+2k^A_E D^R_{12})({\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12})}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}\right) D^A \nonumber \\ \quad&\quad -\left( \frac{(\delta +{\bar{k}}^R_E+k^R_ED^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A){\bar{K}}_{dim}({\bar{k}}^A_W+k^A_MD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}\right) D^R_{12},\nonumber \\ \end{aligned}$$

(20)

with ${\bar{K}}_{dim}=\frac{1}{k^1_{W0}+k^{'}_MD^R_{12}}+\frac{1}{k^2_{W0}+(k_M+{\bar{k}}_M)D^R_{12}}$. This reduced system can be represented by the following chemical reactions:

$$\begin{aligned} {\mathrm{D^{A}} \mathop \rightarrow \limits ^{k_{AR}}\limits \mathrm{D^R_{12}} },\;\;\;{\mathrm{D^R_{12}}} \mathop \rightarrow \limits ^{k_{AR}} \mathrm{D^A} \end{aligned}$$

(21)

with reaction rate coefficients $k_{AR}$ and $k_{RA}$ given by

$$\begin{aligned} k_{AR}&=\frac{(\delta +{\bar{k}}^A_E+2k^A_E D^R_{12})({\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12})}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}, \nonumber \\ k_{RA}&=\frac{(\delta +{\bar{k}}^R_E+k^R_ED^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A){\bar{K}}_{dim}({\bar{k}}^A_W+k^A_MD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_W+{\bar{k}}^1_W+(k_M+{\bar{k}}_M+k_M^{'})D^R_{12}}. \end{aligned}$$

(22)

4.1.1 Mathematical analysis of the stochastic properties

Since the reduced chemical reaction system (21) is characterized by the conservation law $D^R_{12}+D^A\approx D_{{tot}}$, its stochastic behavior can be approximately represented by a one-dimensional Markov chain with state $x=n_{D^R_{12}}\in [0,{\text {D}}_{\textrm{tot}}]$. For any state x, the rate associated with the transition to the next higher state ($x\rightarrow x+1$), $\lambda _x$, and the rate associated with the transition to the next lower state ($x\rightarrow x-1$), $\gamma _x$, can be written as follows:

$$\begin{aligned} \lambda _x&=\left( \frac{(\varepsilon +2\varepsilon '\frac{x}{{\textrm{D}}_{\textrm{tot}}})({\bar{k}}^2_W+{\bar{k}}^1_W+\frac{(k_M+{\bar{k}}_M+k_M^{'})}{\Omega }x)}{{\bar{u}}^A+\frac{({\text {D}}_{\textrm{tot}}-x)}{{\textrm{D}}_{\textrm{tot}}}+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha ')\frac{x}{{\text {D}}_{\textrm{tot}}}}\right) ({{\textrm{D}_\mathrm{{tot}}}}-x),\nonumber \\ \gamma _x&=\left( \frac{\mu (b\varepsilon +\varepsilon ' \frac{({\text {D}}_{\textrm{tot}}-x)}{{\text {D}}_{\textrm{tot}}})\mu '(\beta \varepsilon +\varepsilon ' \frac{({\text {D}}_{\textrm{tot}}-x)}{{\text {D}}_{\textrm{tot}}}){\bar{K}}_x({\bar{k}}^A_W+\frac{k^A_M}{\Omega }({\text {D}}_{\textrm{tot}}-x))}{{\bar{u}}^A+\frac{({\text {D}}_{\textrm{tot}}-x)}{{\text {D}}_{\textrm{tot}}}+{\bar{u}}^R_2+{\bar{u}}^R_1+(\alpha + {\bar{\alpha }} + \alpha ')\frac{x}{{\text {D}}_{\textrm{tot}}}}\right) x. \end{aligned}$$

(23)

Here, we mathematically compute the stationary probability distribution and the time to memory loss of active and repressed chromatin states for this one-dimensional Markov chain.

Proposition 4.2

Let $\lambda _x$ and $\gamma _x$ represent the rate associated with the transition $x\rightarrow x+1$and the rate associated with the transition $x\rightarrow x-1$, respectively. Then, the stationary distribution associated with the one-dimensional Markov chain with rates $\lambda _x$ and $\gamma _x$ can be written as

$$\begin{aligned} \pi (x)= \prod _{i=1}^{x} \frac{\lambda _{i-1}}{\gamma _i}\pi (0)=\frac{\prod _{i=1}^{x} \frac{\lambda _{i-1}}{\gamma _i}}{\left( 1+\sum _{j=1}^{{\textrm{D}}_{\textrm{tot}}}\left( \prod _{i=1}^{j} \frac{\lambda _{i-1}}{\gamma _i}\right) \right) } \end{aligned}$$

(24)

in which $\sum _{x=0}^{{\textrm{D}}_{\textrm{tot}}}\pi (x)=1$.

Proof

By using detailed balance, we can write $\lambda _{x-1}\pi (x-1)=\gamma _x\pi (x)$, for any $x\in [1,{\textrm{D}}_{\textrm{tot}}]$. Then, these equalities can be combined to write $\pi (x)= \prod _{i=1}^{x} (\lambda _{i-1}/\gamma _i)\pi (0)$. Finally, exploiting $\sum _{x=0}^{{\textrm{D}}_{\textrm{tot}}}\pi (x)=1$, we obtain the formula in (24). $\square$

Proposition 4.3

When $\varepsilon \ll 1$, the stationary distribution $\pi (x)$ associated with the one-dimensional Markov chain with rates $\lambda _x$ and $\gamma _x$ as defined in (23) can be written as

$$\begin{aligned} \pi _{\varepsilon \ll 1}(x) \approx {\left\{ \begin{array}{ll} \frac{1}{1+P} &{} {\text {if }}\; x=0 \\ 0 &{} {\text {if }}\; x \ne 0, {\textrm{D}}_{\textrm{tot}}\\ \frac{P}{1+P} &{} {\text {if }}\; x = {\textrm{D}}_{\textrm{tot}} \end{array}\right. } \end{aligned}$$

(25)

with P given by

$$\begin{aligned} \begin{aligned} P&=\frac{({\bar{u}}^A+{\bar{u}}^R_1 + u^R_2+(\alpha + {\bar{\alpha }} + \alpha '))}{({\bar{u}}^A+{\bar{u}}^R_1 + u^R_2+1)}\cdot \frac{{\bar{u}}^R_1 + u^R_2}{\mu \mu 'b\beta \varepsilon {\bar{K}}_{{\textrm{D}}_{\textrm{tot}}}{\bar{u}}^A}\\&\qquad \cdot \prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\left( \frac{2({\bar{u}}^R_1 + u^R_2+(\alpha + {\bar{\alpha }} + \alpha ')\frac{i}{{{\textrm{D}_\mathrm{{tot}}}}})}{\mu \mu '\varepsilon ' \frac{({{\textrm{D}_\mathrm{{tot}}}}-i)}{{{\textrm{D}_\mathrm{{tot}}}}}{\bar{K}}_i({\bar{u}}^A+\frac{({{\textrm{D}_\mathrm{{tot}}}}-i)}{{{\textrm{D}_\mathrm{{tot}}}}})}\right) , \end{aligned} \end{aligned}$$

and ${\bar{K}}_{{{\textrm{D}_\mathrm{{tot}}}}}=\frac{1}{u^R_{10}+\alpha '}+\frac{1}{u^R_{20}+(\alpha +{\bar{\alpha }})}$.

Proof

Assuming that $\varepsilon '\ne 0$ and $\varepsilon \ll 1$, for any $j\in [1,{{\textrm{D}_\mathrm{{tot}}}}-1]$ we have that $\prod _{i=1}^{j} \frac{\lambda _{i-1}}{\gamma _i} \ll \prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}} \frac{\lambda _{i-1}}{\gamma _i}$, with $\lambda _x$ and $\gamma _x$ provided in (23). This implies that when $\varepsilon \ll 1$, the sum $\sum _{j=0}^{{{\textrm{D}_\mathrm{{tot}}}}}\pi (j)=1$ can be approximated as

$$\begin{aligned} 1=\sum _{j=0}^{{{\textrm{D}_\mathrm{{tot}}}}}\pi (j)&=\left[ \sum _{j=1}^{{{\textrm{D}_\mathrm{{tot}}}}}\left( \prod _{i=1}^{j} \frac{\lambda _{i-1}}{\gamma _i}\right) \right] \pi (0)+\pi (0)\\&\approx \prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}} \frac{\lambda _{i-1}}{\gamma _i}\pi (0)+\pi (0)=\pi ({{\textrm{D}_\mathrm{{tot}}}})+\pi (0) \end{aligned}$$

from which, introducing the notation $P=\prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}}(\lambda _{i-1})/(\gamma _i)$ and writing explicitly $\lambda _x$ and $\gamma _x$, the stationary distribution formula (24) can be rewritten as done in (25). $\square$

By studying the expression for $\pi _{\varepsilon \ll 1}(x)$ in (25), it is possible to notice that if $\varepsilon \ll 1$, the only states in which $\pi (x)$ does not vanish are the fully active state $x=0$ and the fully repressed state $x={{\textrm{D}_\mathrm{{tot}}}}$. More precisely, the stationary distribution for $\varepsilon \ll 1$ is bimodal, with two modes in correspondence to $x=0$ and $x={{\textrm{D}_\mathrm{{tot}}}}$, and the probability of having the system in one of the intermediate states is approximately zero. Furthermore, when $\varepsilon$ decreases, P increases and accordingly $\pi _{\varepsilon \ll 1}({{\textrm{D}_\mathrm{{tot}}}})$ increases to the detriment of $\pi _{\varepsilon \ll 1}(0)$. This result is in agreement with the structural asymmetry toward a repressed chromatin state characterizing the chromatin modification circuit because of the cooperation between H3K9me3 and DNA methylation (Fig. 1b). This result (Proposition 4.3) implies that $\varepsilon$ plays crucial role in the duration of memory of the active and repressed states and that, when it is small, the duration of memory increases.

In order to make mathematically precise this qualitative statement, we determine an expression for the temporal duration of the memory of the fully repressed and fully active chromatin states and study how $\varepsilon$ affects it. First, let us provide the definition of time to memory loss and then let us introduce the expression for the time to memory loss of the active and repressed states:

Definition 4.1

(Time to memory loss) Let $t_i^j$ represent the hitting time of $x=j$ starting from $x=i$, that is, $t_i^j :=$ inf$\{t \ge 0:x(t) = j$ with $x(0) = i\}$ with $i,j\in [0,{{\textrm{D}_\mathrm{{tot}}}}]$, where x(t) is the Markov chain described above. The time to memory loss of the fully repressed chromatin state is defined as $\tau _{{{\textrm{D}_\mathrm{{tot}}}}}^0=\mathbb {E}(t_{{{\textrm{D}_\mathrm{{tot}}}}}^0)$. Similarly, the time to memory loss of the active state is defined as $\tau _{0}^{{{\textrm{D}_\mathrm{{tot}}}}}=\mathbb {E}(t_{0}^{{{\textrm{D}_\mathrm{{tot}}}}})$.

Proposition 4.4

The time to memory loss of the repressed chromatin state is given by

$$\begin{aligned} \begin{aligned} \tau _{{{\textrm{D}_\mathrm{{tot}}}}}^0&=\frac{s_{{{\textrm{D}_\mathrm{{tot}}}}-1}}{\gamma _{{{\textrm{D}_\mathrm{{tot}}}}}}\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{1}{s_x}\right) +\frac{1}{\gamma _1}+\sum _{x=2}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\left[ \frac{s_{x-1}}{\gamma _x}\left( 1+\sum _{j=1}^{x-1}\frac{1}{s_j}\right) \right] , \end{aligned} \end{aligned}$$

(26)

in which $s_x=\frac{\lambda _1 \lambda _2\ldots \lambda _x}{\gamma _1 \gamma _2 \ldots \gamma _x}$ and $\lambda _x$ and $\gamma _x$ are defined in (23). The time to memory loss of the active chromatin state is given by

$$\begin{aligned} \begin{aligned} \tau _{0}^{{\textrm{D}_\mathrm{{tot}}}}&=\frac{{\tilde{s}}_{{{\textrm{D}_\mathrm{{tot}}}}-1}}{\lambda _{0}}\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{1}{{\tilde{s}}_x}\right) +\frac{1}{\lambda _{{{\textrm{D}_\mathrm{{tot}}}}-1}}+\sum _{x=2}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\left[ \frac{{\tilde{s}}_{x-1}}{\lambda _{{{\textrm{D}_\mathrm{{tot}}}}-x}}\left( 1+\sum _{j=1}^{x-1}\frac{1}{{\tilde{s}}_j}\right) \right] , \end{aligned}\nonumber \\ \end{aligned}$$

(27)

in which ${\tilde{s}}_x=\frac{\gamma _{{{\textrm{D}_\mathrm{{tot}}}}-1}\gamma _{{{\textrm{D}_\mathrm{{tot}}}}-2}\ldots \gamma _{{{\textrm{D}_\mathrm{{tot}}}}-x}}{\lambda _{{{\textrm{D}_\mathrm{{tot}}}}-1}\lambda _{{{\textrm{D}_\mathrm{{tot}}}}-2}\ldots \lambda _{{{\textrm{D}_\mathrm{{tot}}}}-x}}$.

Proof

By using first step analysis [11], we can write the following equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} \tau _i^0=0 &{} \text{ if } i=0 \\ (\lambda _i+\gamma _i)\tau _i^0-\lambda _i\tau _{i+1}^0-\gamma _i\tau _{i-1}^0=1 &{} \text{ if } i\in [1,{{\textrm{D}_\mathrm{{tot}}}}-1]\\ \gamma _{i}\tau _i^0-\gamma _{i}\tau _{i-1}^0=1 &{} \text{ if } i= {{\textrm{D}_\mathrm{{tot}}}}\\ \end{array}\right. } \end{aligned}$$

(28)

Then, by solving system (28), we obtain the formula for $\tau _{{{\textrm{D}_\mathrm{{tot}}}}}^0$ as in (26). A similar approach can be used to obtain the formula for $\tau _{0}^{{\textrm{D}_\mathrm{{tot}}}}$ as in (27). $\square$

Proposition 4.5

Assuming $\varepsilon '\ne 0$ and normalizing $\tau ^0_{{{D_\mathrm{{tot}}}}}$ and $\tau ^{{{\textrm{D}_\mathrm{{tot}}}}}_0$ with respect to $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ (${\bar{\tau }}^0_{{{\textrm{D}_\mathrm{{tot}}}}} = \tau ^0_{{{\textrm{D}_\mathrm{{tot}}}}}\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ and ${\bar{\tau }}^{{{\textrm{D}_\mathrm{{tot}}}}}_0 = \tau ^{{{\textrm{D}_\mathrm{{tot}}}}}_0\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$), the times to memory loss (26) and (27) in the regime $\varepsilon \ll 1$ can be approximated with the following expressions:

$$\begin{aligned} \begin{aligned} {\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0\approx \frac{G_R}{\mu \mu '\varepsilon ^{2}}\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{G^x_R}{g_1^x(\mu \mu ')}\right) , \end{aligned} \end{aligned}$$

(29)

$$\begin{aligned} \begin{aligned} {\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}\approx \frac{G_A}{\varepsilon }\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{g_2^x(\mu \mu ')}{G^x_A}\right) , \end{aligned} \end{aligned}$$

(30)

in which $g_1^x(\mu \mu ')$ and $g_2^x(\mu \mu ')$ are increasing functions of $\mu \mu '$ with $g_1^x(0)=0$ and $g_2^x(0)=0$, respectively, and $G^x_R$, $G_R$, $G_A$ and $G^x_A$ are functions that do not depend on $\varepsilon$, $\mu '$ and $\mu$.

Proof

By multiplying by $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ the expressions for times to memory loss given in Prop. 4.4, with $\lambda _x$ and $\gamma _x$ defined in (23), we obtain the normalized expressions for times to memory loss. Then, by approximating them with their dominant term (the term $O(1/\varepsilon ^2))$ for ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and the term $O(1/\varepsilon )$ for ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$, respectively), we obtain the expressions (29) and (30). $\square$

By studying the expressions for ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ in (29) and (30), respectively, it is possible to notice that decreasing $\varepsilon$ increases both ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$, implying that lower $\varepsilon$ extends the duration of memory of both the active and repressed chromatin states. However, given the cooperation of the repressive marks and the consequent structural asymmetry of the chromatin modification circuit, ${\bar{\tau }}_{{\textrm{D}_\mathrm{{tot}}}}^0=O(1/\varepsilon ^2)$, while ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}=O(1/\varepsilon )$. This implies that decreasing $\varepsilon$ extends more the repressed state memory than the active state memory.

Now, let us also determine the effect of the asymmetry between the erasure rates of repressive and activating chromatin modifications, encapsulated by the non-dimensional parameters $\mu$ and $\mu '$. From the expression for the stationary distribution in (25) it is possible to notice that, by reducing $\mu '$ or $\mu$ (i.e., reducing the erasure rates of the repressive marks compared to the erasure rate of the active marks), $\pi _{\varepsilon \ll 1}({{\textrm{D}_\mathrm{{tot}}}})$ increases to the detriment of $\pi _{\varepsilon \ll 1}(0)$, i.e., the stationary distribution shifts toward the repressed state. In agreement with these results, when $\mu '$ or $\mu$ decreases, ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ increases, while ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ decreases, that is, the temporal duration of the memory of the repressed state increases, while the duration of memory of the active state decreases.

Fig. 2

Computational analysis of the full chromatin modification circuit, shown in Fig. 1b, using SSA. a The stationary probability distribution, $\pi$, for the chromatin modification circuit represented in Fig. 1b, whose reactions are listed in Fig. 1a. The parameter values considered to generate the plots are in Table S1 of S1 File. In particular, in the left-side plots $\varepsilon =0.28,0.14$, $\mu '=0.675,0.35$ and $\varepsilon '=1$ and in the right-side plots $\varepsilon =0.28,0.14$, $\mu '=0.625,0.35$ and $\varepsilon '=0.4$. In all plots $n_{D^A}$ and $n_{D^R}=n_{D^R_1}+n_{D^R_2}+n_{D^R_{12}}$ represent the number of nucleosomes with activating and repressive modifications. b The stationary distribution for the chromatin modification circuit for different values of $\varepsilon '$. The parameter values considered are listed in Table S1 of S1 File. In particular, $\varepsilon =0.28$ and $\varepsilon '=1,0.01$. c Time trajectories of $n_{D^A}$ and $n_{D^R}$ starting from the fully active state $n_{D^A}=50$, $n_{D^R}=0$ (left) and repressed state $n_{D^A}=0$, $n_{D^R}=n_{D^R_{12}}=50$ (right) for $\varepsilon '=1$ and different values of $\varepsilon$. d Time trajectories of $n_{D^A}$ and $n_{D^R}$, as described in c, but with $\varepsilon '=0.4$. e Time trajectories of the system starting from $n_{D^R}=n_{D^R_{12}}=50,n_{D^A}=0$ and with an input $u^A$ that, at steady state, leads to a unimodal distribution near the active state $n_{D^A}\approx 50$. Each trajectory is represented with a different color. In particular, we set $\varepsilon =0.28$, $\varepsilon '=1$, $\mu =1$ and $\mu '=0.675,0.35$. In c–e, the time is normalized ($\tau =t\frac{k^A_M}{\Omega }{{\textrm{D}_\mathrm{{tot}}}}$, with $\Omega$ the reaction volume) and the parameter values are listed in Table S1 of S1 File

4.1.2 Computational analysis

In this section, we validate the trends determined by the analytical study in the previous section, which exploits a deterministic quasi-steady state approximation [20], and we demonstrate the validity of these results for a broader parameter regime than $\varepsilon ' \ll 1$ and $\varepsilon =c\varepsilon '$, with $c=O(1)$. To this end, we employ the stochastic simulation algorithm (SSA) [12] to study via simulation the original chemical reaction system represented in Fig. 1b, whose reactions are listed in Fig. 1a.

The trend with which $\varepsilon$ and $\mu '$ affect the stationary distribution of the original system $\pi (x)$ is in agreement with the results obtained from the analytical study in Sect. 4.1.1. The parameter $\varepsilon '$ does not significantly vary the way in which $\varepsilon$, $\mu '$, $\mu$ affect the stationary distribution. However, decreasing $\varepsilon '$ compared to $\varepsilon$ leads to less concentrated peaks in the bimodal stationary distribution and, by further decreasing $\varepsilon '$, the distribution can become unimodal (Fig.2b). Now, let us consider a parameter regime in which the system displays a bimodal distribution and let us study how the switching time of the system temporal trajectories depends on $\varepsilon$ (Fig. 2c, d). In particular, in agreement with our analytical findings, it is possible to notice that lowering $\varepsilon$ increases both the time that the system spends at the active state before switching to the repressed state and the time that the system spends at the repressed state before switching to the active state, but the latter one is a stronger effect.

We next determine via simulation how the parameter $\mu '$ affects the reactivation time; that is, the time needed to re-activate an initially repressed chromatin state after a sufficiently large activating input stimulus $u^A$ has been applied (Fig. 2e). It is possible to notice that the time trajectories show a switch-like behavior and that the time needed to see the trajectory switching to the active state after an activating input is applied becomes more variable for small $\mu '$. This implies that for low values of $\mu '$ the reactivation of the gene becomes a very stochastic process.

4.2 Limiting behavior for large $\mu '$

In this section, we determine the stochastic behavior of the chromatin modification circuit when $\mu '$ is large (i.e., when the erasure of DNA methylation is much faster than the erasure of histone modifications). Comparing the results here with those of the previous section allows us to determine how the absence of DNA methylation affects the system’s stochastic features. Since we are interested in the behavior of the system for large $\mu '$, we conduct again the model reduction by assuming that the product $\varepsilon '\mu '$ does not vanish as $\varepsilon '$ approaches zero. This leads to a different reduced model compared to the one obtained in Sect. 4.1. To this end, we define the parameter $U':=\varepsilon '\mu '$ and assume that it is $\varepsilon '$-independent. Then, we introduce $U'$ in the ODE model (4), perform the model reduction with $\varepsilon '$ as a small parameter and consider the limiting case $U^{'}\rightarrow \infty$. Then, we conduct an analytical study of the stochastic behavior of the reduced system, and validate and extend the results obtained via simulation.

4.2.1 Model reduction

Proposition 4.6

Let $\varepsilon =c\varepsilon '$, with $c=O(1)$, and let $U':=\varepsilon '\mu '$ be $\varepsilon '$-independent. Then, let us consider the following system:

$$\begin{aligned} \frac{{\text {d}}\bar{D}^A}{{\text {d}}{\tau }}&=\left( \frac{\mu (b\varepsilon +\varepsilon '{\bar{D}}^A)(u_A + {\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2 + \alpha {\bar{D}}^R_2}\right) {\bar{D}}^R_2-\left( \frac{(\varepsilon +\varepsilon '{\bar{D}}^R_2)({\bar{u}}^R_2 + \alpha {\bar{D}}^R_2)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2 + \alpha {\bar{D}}^R_2}\right) {\bar{D}}^A\nonumber \\ \frac{{\text {d}}\bar{D}^R_2}{{\text {d}}{\tau }}&=\left( \frac{(\varepsilon +\varepsilon '{\bar{D}}^R_2)({\bar{u}}^R_2 + \alpha {\bar{D}}^R_2)}{{\bar{u}}^A+{\bar{D}}^A+u^R_2 + \alpha {\bar{D}}^R_2}\right) {\bar{D}}^A-\left( \frac{\mu (b\varepsilon +\varepsilon '{\bar{D}}^A)(u_A + {\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_2 + \alpha {\bar{D}}^R_2}\right) {\bar{D}}^R_2, \end{aligned}$$

(31)

with ${\bar{D}}^A + {\bar{D}}^R_{2}=1$. Furthermore, let $({\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_{12}(\tau ))=M({\bar{D}}^A(\tau ),{\bar{D}}^R_2(\tau ))$ represent the system’s unique integral manifold. Then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau )$, ${\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))$ of (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_1(0),{\bar{D}}^R_{12}(0))-M({\bar{D}}^A(0),{\bar{D}}^R_2(0))|$ is sufficiently small, we have that, for $\varepsilon '$ sufficiently small, the solution of (31), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_2(\tau ))$, is such that

$$\begin{aligned} \begin{aligned} ({\bar{D}}^A(\tau ),{\bar{D}}^R_2(\tau ))&=({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_2(\tau ))+\zeta _1(\tau ),\\ ({\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_{12}(\tau ))&=M({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_2(\tau ))+\zeta _2(\tau ), \end{aligned} \end{aligned}$$

(32)

in which $\zeta _i(\tau )=O(e^{-(\theta /\varepsilon ')\tau })$, with $i=1,2$, $\theta >0$, and $\tau \ge 0$.

Proof

Let us introduce $U'=\varepsilon '\mu '$ in system (4), obtaining

$$\begin{aligned} \varepsilon '\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}&= ({\bar{u}}^A+{\bar{D}}^A){\bar{D}}-\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}^R_{12}}{\hbox {d}\bar{\tau }}&=(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^R_2+(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}^R_1\nonumber \\&\quad -(\varepsilon '\mu (b c+ \bar{D^A})+U'(\beta c+ \bar{D^A})){\bar{D}}^R_{12}\nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}^R_1}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\nonumber \\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+U'(\beta c+ \bar{D^A})){\bar{D}}^R_1 \nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}+U'(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\nonumber \\&\quad -(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+\varepsilon '\mu (b c+ \bar{D^A})){\bar{D}}^R_2\nonumber \\ \varepsilon '\frac{\hbox {d}\bar{D}}{\hbox {d}\bar{\tau }}&=(U'(\beta c+ \bar{D^A}){\bar{D}}^R_1+\varepsilon '\mu (b c+ \bar{D^A}){\bar{D}}^R_2)\nonumber \\&\quad +\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A\nonumber \\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+{\bar{u}}^R_1\nonumber \\&\quad +\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A){\bar{D}}. \end{aligned}$$

(33)

Now, let us define x, $y_2$, $f_1$ and $f_2$ as

$$\begin{aligned} x&= \begin{pmatrix} {\bar{D}}^A\\ {\bar{D}}^R_2\\ \end{pmatrix},y_2=\begin{pmatrix} {\bar{D}}^R_{12}\\ {\bar{D}}^R_1\\ {\bar{D}}\\ \end{pmatrix},f_1=\begin{pmatrix} f_{11}\\ f_{12}\\ \end{pmatrix},f_2=\begin{pmatrix} f_{21}\\ f_{22}\\ f_{23}\\ \end{pmatrix},\\ f_{11}&=({\bar{u}}^A+{\bar{D}}^A){\bar{D}}-\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^A,\\ f_{12}&=({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}+U'(\beta c+ \bar{D^A}){\bar{D}}^R_{12}\\&\quad -(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+\varepsilon '\mu (b c+ \bar{D^A})){\bar{D}}^R_2\\ f_{21}&=(u^R_{10}+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}^R_2+(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})){\bar{D}}^R_1\\&\quad -(\varepsilon '\mu (b c+ \bar{D^A})+U'(\beta c+ \bar{D^A})){\bar{D}}^R_{12},\\ f_{22}&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})){\bar{D}}+\varepsilon '\mu (b c+ {\bar{D}}^A){\bar{D}}^R_{12}\\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1+{\bar{D}}^R_{12})+U'(\beta c+ \bar{D^A})){\bar{D}}^R_1,\\ f_{23}&=(U'(\beta c+ \bar{D^A}){\bar{D}}^R_1+\varepsilon '\mu (b c+ \bar{D^A}){\bar{D}}^R_2)\\&\quad +\varepsilon '(c+({\bar{D}}^R_1+{\bar{D}}^R_{12})+({\bar{D}}^R_2+{\bar{D}}^R_{12}){\bar{D}}^A\\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{\alpha }}({\bar{D}}^R_1\nonumber \\&\quad +{\bar{D}}^R_{12})+{\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+{\bar{D}}^R_{12})+{\bar{u}}^A+{\bar{D}}^A){\bar{D}}. \end{aligned}$$

Now, it is possible to calculate that $\phi (x)=(0,0,\phi _{12})$, with $\phi (x)$ defined in C1 and with $\phi _{12}=\frac{(u^R_{10}+\alpha '{\bar{D}}^R_2){\bar{D}}^R_2}{U'(\beta c+{\bar{D}}^A)-\alpha '{\bar{D}}^R_2}$. Function $\phi _{12}$ is inversely proportional to $U'$. Furthermore, the matrix A, defined in (6), with ${\bar{D}} = {\bar{D}}^R_1=0$, ${\bar{D}}^R_{12}=\phi _{12}$ and $\varepsilon '=0$ can be written as

$$\begin{aligned} A(x,y_2=\phi (x),t,0)= \left( \begin{array}{cc} {\bar{A}}_{2,2}&{}\qquad {\bar{A}}_{2,3}\\ {\bar{A}}_{3,2}&{}\qquad {\bar{A}}_{3,3}\\ \end{array}\right) \end{aligned}$$

(34)

with

$$\begin{aligned} {\bar{A}}_{2,2}&= \left( {\begin{matrix} 0&{}\qquad 0\\ U'\phi _{12}&{}\qquad -\alpha ' {\bar{D}}^R_2-(u^R_{10}+\alpha '({\bar{D}}^R_2+\phi _{12}))\\ \end{matrix}}\right) ,\\ {\bar{A}}_{2,3}&= \left( {\begin{matrix} 0&{}\qquad 0&{}\qquad ({\bar{u}}^A+{\bar{D}}^A)\\ U'(\beta c+{\bar{D}}^A)-\alpha '{\bar{D}}^R_2&{}\qquad 0&{}\qquad {\bar{u}}_2+\alpha ({\bar{D}}^R_2+\phi _{12})+{\bar{\alpha }} \phi _{12}\\ \end{matrix}}\right) ,\\ {\bar{A}}_{3,2}&= \left( {\begin{matrix} -U'\phi _{12}&{}\qquad \alpha ' {\bar{D}}^R_2+(u^R_{10}+\alpha '({\bar{D}}^R_2+\phi _{12}))\\ 0&{}\qquad 0\\ 0&{}\qquad 0\\ \end{matrix}}\right) ,\\ {\bar{A}}_{3,3}&= \left( {\begin{matrix} {\hat{A}}_{1,1}&{}\qquad {\hat{A}}_{1,2}\\ 0_{2,1}&{}\qquad {\hat{A}}_{2,2}\\ \end{matrix}}\right) ,\;\; ({\hat{A}}_{1,1} {\hat{A}}_{1,2}) = \left( {\begin{matrix} -U'(\beta c+{\bar{D}}^A)+\alpha '{\bar{D}}^R_2 &{} u^R_{20}+\alpha ({\bar{D}}^R_2+\phi _{12})+{\bar{\alpha }} \phi _{12}&{} 0\\ \end{matrix}}\right) ,\\ {\hat{A}}_{2,2}&= \left( {\begin{matrix} {\tilde{A}}_{1,2} \\ {\tilde{B}}_{1,2}\\ \end{matrix}}\right) ,\;\; {\tilde{A}}_{1,2}= \left( {\begin{matrix} -(u^R_{20}+\alpha ({\bar{D}}^R_2+\phi _{12})+{\bar{\alpha }} \phi _{12}+U'(\beta c +{\bar{D}}^A))&-(u^R_{10}+\alpha '({\bar{D}}^R_2 + \phi _{12})) \end{matrix}}\right) ,\\ {\tilde{B}}_{1,2}&= \left( {\begin{matrix} U'(\beta c +{\bar{D}}^A)&-({\bar{u}}_2 + \alpha ({\bar{D}}^R_2+\phi _{12})+{\bar{\alpha }} \phi _{12}+{\bar{u}}_1 + \alpha ' ({\bar{D}}^R_2+\phi _{12})+{\bar{u}}^A + {\bar{D}}^A) \end{matrix}}\right) . \end{aligned}$$

The matrix (34) is singular, and this implies that the system (33) is singular singularly perturbed (Def. 3.2). Specifically, matrix A has a twofold zero eigenvalue, and two linearly independent eigenvectors associated with them, and matrix $B={\bar{A}}_{3,3}$, with B defined as in (7). When there are no external inputs ($u^A=u^R_1=u^R_2=0$ and then ${\bar{u}}^A=u^A_0$, ${\bar{u}}^R_1=u^R_{10}$, and ${\bar{u}}^R_2=u^R_{20}$), matrix B has three eigenvalues with negative real part if $u^R_{10}, u^R_{20}, u^A_0\ge l$ with $l>0$. This implies that we can apply Theorem 3.1 to reduce our system. To do that, let us first introduce the variable ${\hat{D}}^R_{12}={\bar{D}}^R_{12}-\phi _{12}$, and then the variables ${\tilde{D}}={\bar{D}}/ \varepsilon '$, ${\tilde{D}}^R_1={\bar{D}}^R_1/ \varepsilon '$ and ${\tilde{D}}^R_{12}={\hat{D}}^R_{12}/ \varepsilon '$ in (33):

$$\begin{aligned} \varepsilon '\frac{\hbox {d}\tilde{D}^R_1}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\tilde{D}}+\mu (b c+ {\bar{D}}^A)(\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})\nonumber \\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}\nonumber \\&\quad +\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+U^{'}(\beta c+ \bar{D^A})){\tilde{D}}^R_1\nonumber \\ \varepsilon '\frac{\hbox {d}\tilde{D}}{\hbox {d}\bar{\tau }}&=U'(\beta c+ \bar{D^A}){\tilde{D}}^R_1+\mu (b c+ \bar{D^A}){\bar{D}}^R_2\nonumber \\&\quad +(c+(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12}+{\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\bar{D}}^A\\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\tilde{D}}\nonumber \\&\quad -({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+{\bar{u}}^A+{\bar{D}}^A){\tilde{D}}\nonumber \\ \varepsilon ' \frac{\hbox {d}\tilde{D}^R_{12}}{\hbox {d}\bar{\tau }}&=-\frac{\hbox {d}\phi _{12}}{\hbox {d}\bar{\tau }}+(u^R_{20}+\alpha ({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\tilde{D}}^R_1\nonumber \\&\quad +\alpha '{\bar{D}}^R_2 {\tilde{D}}^R_{12}- \mu (b c +{\bar{D}}^A)(\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})-U'(\beta c + {\bar{D}}^A){\tilde{D}}^R_{12}\nonumber \\ \frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}&=({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\tilde{D}}\nonumber \\&\quad +U'(\beta c+ \bar{D^A}){\tilde{D}}^R_{12}-(\alpha '{\tilde{D}}^R_{12}+\mu (bc+{\bar{D}}^A)){\bar{D}}^R_2\nonumber \\ \frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}&= ({\bar{u}}^A+{\bar{D}}^A){\tilde{D}}-(c+(\varepsilon ^{' }{\tilde{D}}^R_1+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})+({\bar{D}}^R_2+\varepsilon ' {\tilde{D}}^R_{12}+\phi _{12})){\bar{D}}^A.\nonumber \end{aligned}$$

(35)

with

$$\begin{aligned} \frac{\hbox {d}\phi _{12}}{\hbox {d}\bar{\tau }}&=\frac{\partial \phi _{12}}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial \phi _{12}}{\partial {\bar{D}}^R_2}\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}=-\frac{U'(u^R_{10}+\alpha '{\bar{D}}^R_2){\bar{D}}^R_2}{(U'(\beta c+{\bar{D}}^A)-\alpha '{\bar{D}}^R_2)^2}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}\\&\quad +\frac{(u^R_{10}+2\alpha '{\bar{D}}^R_2)(U'(\beta c+{\bar{D}}^A)-\alpha '{\bar{D}}^R_2)+(u^R_{10}+\alpha '{\bar{D}}^R_2)\alpha ' {\bar{D}}^R_2}{(U'(\beta c+{\bar{D}}^A)-\alpha '{\bar{D}}^R_2)^2}\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}. \end{aligned}$$

Now, similarly to what we did in Sect. 4.1, in order to determine the integral manifold $M({\bar{D}}^A,{\bar{D}}^R_2)=({\bar{D}},{\bar{D}}^R_1,{\bar{D}}^R_{12})$, we evaluate the asymptotic expansion of ${\tilde{D}}$, ${\tilde{D}}^R_1$ and ${\tilde{D}}^R_{12}$, that can be written as follows:

$$\begin{aligned} {\tilde{D}}&= h_0({\bar{D}}^A, {\bar{D}}^R_2, \varepsilon ')= h_{00}({\bar{D}}^A, {\bar{D}}^R_2) + \varepsilon 'h_{01}({\bar{D}}^A, {\bar{D}}^R_2)+O(\varepsilon ^{'^2}),\\ {\tilde{D}}^R_1&= h_1({\bar{D}}^A, {\bar{D}}^R_2, \varepsilon ')= h_{10}({\bar{D}}^A, {\bar{D}}^R_2) + \varepsilon 'h_{11}({\bar{D}}^A, {\bar{D}}^R_2)+O(\varepsilon ^{'^2}),\nonumber \\ {\tilde{D}}^R_{12}&= h_2({\bar{D}}^A, {\bar{D}}^R_2, \varepsilon ')= h_{20}({\bar{D}}^A, {\bar{D}}^R_2) + \varepsilon 'h_{21}({\bar{D}}^A, {\bar{D}}^R_2)+O(\varepsilon ^{'^2}).\nonumber \end{aligned}$$

(36)

To this end, we plug (36) into the first three equations of (35), obtaining

$$\begin{aligned} \varepsilon '\frac{\hbox {d}h_0}{\hbox {d}\bar{\tau }}&=\varepsilon '\left( \frac{\partial h_0}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_0}{\partial {\bar{D}}^R_2}\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}\right) \\&=U'(\beta c+ \bar{D}^A)h_1+\mu (b c+ \bar{D}^A){\bar{D}}^R_2\\&\quad +(c+(\varepsilon ^{' }h_1+\varepsilon ' h_2+\phi _{12}+{\bar{D}}^R_2+\varepsilon ' h_2+\phi _{12})){\bar{D}}^A\\&\quad -({\bar{u}}^R_2+\alpha ({\bar{D}}^R_2+\varepsilon ' h_2+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+\varepsilon ' h_2+\phi _{12}))h_0\\&\quad -({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+\varepsilon ' h_2+\phi _{12})+{\bar{u}}^A+{\bar{D}}^A)h_0\\ \varepsilon '\frac{\hbox {d}h_1}{\hbox {d}\bar{\tau }}&=\varepsilon '(\frac{\partial h_1}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_1}{\partial {\bar{D}}^R_{2}}\frac{\hbox {d}\bar{D}^R_{2}}{\hbox {d}\bar{\tau }})\\&=({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2+\varepsilon ^{' }h_2+\phi _{12}))h_0+\mu (b c+ {\bar{D}}^A)(\varepsilon ^{' }h_2+\phi _{12})\\&\quad -(u^R_{20}+\alpha ({\bar{D}}^R_2+\varepsilon ^{' }h_2+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+\varepsilon ^{' } h_2+\phi _{12})+U'(\beta c+ \bar{D^A}))h_1\\ \varepsilon '\frac{\hbox {d}h_2}{\hbox {d}\bar{\tau }}&=\varepsilon '\left( \frac{\partial h_2}{\partial {\bar{D}}^A}\frac{\hbox {d}\bar{D}^A}{\hbox {d}\bar{\tau }}+\frac{\partial h_2}{\partial {\bar{D}}^R_2}\frac{\hbox {d}\bar{D}^R_2}{\hbox {d}\bar{\tau }}\right) \\&=-\frac{\hbox {d}\phi _{12}}{\hbox {d}\bar{\tau }}+(u^R_{20}+\alpha ({\bar{D}}^R_2+\varepsilon ' h_2+\phi _{12})+{\bar{\alpha }}(\varepsilon ^{' }h_1+\varepsilon ' h_2+\phi _{12}))h_1\\&\quad +\alpha '{\bar{D}}^R_2 h_2-\mu (b c +{\bar{D}}^A)(\varepsilon ' h_2+\phi _{12})-U'(\beta c + {\bar{D}}^A)h_2. \end{aligned}$$

Now, we can obtain $h_{i0}$ and $h_{i1}$, with $i=0,1,2$, by equating the terms of the right and left hand side of the equations above multiplied by the same power of $\varepsilon '$. Specifically, since $\frac{\partial h_{i0}}{\partial {\bar{D}}^R_2}$ and $\frac{\partial h_{i0}}{\partial {\bar{D}}^A}$ are bounded for any $i=0,1,2$ (i.e., $\varepsilon '\frac{\partial h_{i0}}{\partial {\bar{D}}^R_2},\varepsilon '\frac{\partial h_{i0}}{\partial {\bar{D}}^A}\ll 1$ for sufficiently small $\varepsilon '$) and since for $U'\gg 1$ we have that $\phi _{12}\ll 1$ and $d\phi _{12}/d\bar{\tau }\approx (u^R_{10}+2\alpha '{\bar{D}}^R_2)h_2$, we can rewrite the expressions for $h_{i0}$ and $h_{i1}$, $i=0,1,2$, as follows:

$$\begin{aligned} h_{00}&=\frac{\mu (bc + {\bar{D}}^A){\bar{D}}^R_2+(c+{\bar{D}}^R_2){\bar{D}}^A}{{\bar{u}}^R_2+\alpha {\bar{D}}^R_2+{\bar{u}}^A+{\bar{D}}^A},\nonumber \\ h_{10}&=\frac{({\bar{u}}^R_1+\alpha '({\bar{D}}^R_2 + \phi _{12}))h_{00}}{U'(\beta c+{\bar{D}}^A)},\;\;\;\;\;\; h_{20}=\frac{(u^R_{20}+\alpha {\bar{D}}^R_2)h_{10}}{U'(\beta c + {\bar{D}}^A)+(u^R_{10}+2\alpha '{\bar{D}}^R_2)},\nonumber \\ h_{01}&=\frac{U'(\beta c + {\bar{D}}^A)h_{11}+(h_{10}+2h_{20}){\bar{D}}^A-((\alpha +{\bar{\alpha }}+\alpha ')h_{20}+{\bar{\alpha }} h_{10})h_{00}}{{\bar{u}}^R_2+\alpha {\bar{D}}^R_2 + {\bar{u}}^R_1 + \alpha '{\bar{D}}^R_2+{\bar{u}}^A+{\bar{D}}^A},\\ h_{11}&=\frac{({\bar{u}}^R_1+\alpha '{\bar{D}}^R_2)h_{01}+\alpha 'h_{20}h_{00}+\mu (bc+{\bar{D}}^A)h_{20}-({\bar{\alpha }} h_{10}+(\alpha + {\bar{\alpha }})h_{20})h_{10}}{U'(\beta c+{\bar{D}}^A)},\nonumber \\ h_{21}&=\frac{({\bar{\alpha }} h_{10}+(\alpha +{\bar{\alpha }})h_{20})h_{10}+(u^R_{20}+\alpha {\bar{D}}^R_2)h_{11}-\mu (b c+{\bar{D}}^A)h_{20}}{U'(\beta c+{\bar{D}}^A)+(u^R_{10}+2\alpha '{\bar{D}}^R_2)}.\nonumber \end{aligned}$$

(37)

Then, by considering the limiting condition $U' \rightarrow \infty$, the expressions for $h_{i0}$ and $h_{i1}$ can be approximated as

$$\begin{aligned}&h_{00}=\frac{\mu (bc + {\bar{D}}^A){\bar{D}}^R_2+(c+{\bar{D}}^R_2){\bar{D}}^A}{{\bar{u}}^R_2+\alpha {\bar{D}}^R_2+{\bar{u}}^A+{\bar{D}}^A},\;\;\;h_{10}= h_{20}=h_{01}=h_{11}=h_{21}=0. \end{aligned}$$

(38)

Now, by plugging the asymptotic expansion of ${\tilde{D}}$, ${\tilde{D}}^R_1$ and ${\tilde{D}}^R_{12}$ (36) with the expressions for $h_{i0}$ and $h_{i1}$, $i=0,1,2$, given in (38), into the last two ODEs of (35), we obtain the reduced system in (31), in which we have re-introduced the original time variable $\tau ={\bar{\tau }} / \varepsilon '$. It is possible to notice that the sum of the two ODEs in (31) is equal to zero, implying that $\bar{D}^A+\bar{D}^R_{2}=$ constant. Since the conservation law $\bar{D}^A+\bar{D}^R_{12}+\bar{D}+\bar{D}^R_1+\bar{D}^R_2=1$ holds and, for sufficiently small $\varepsilon '$ and sufficiently large $U'$, $\bar{D}=\varepsilon ' {\tilde{D}} \approx 0$, $\bar{D}^R_1=\varepsilon ' \tilde{D}^R_1 \approx 0$, and $\bar{D}^R_{12}=\varepsilon ' \tilde{D}^R_{12} \approx 0$, then ${\bar{D}}^A + {\bar{D}}^R_2$ can be approximately set equal to 1 for sufficiently small values of $\varepsilon '$ and sufficiently large $U^{'}$.

Furthermore, given that the integral manifold obtained, $M({\bar{D}}^A,{\bar{D}}^R_2)=(\bar{D},{\bar{D}}^R_1,{\bar{D}}^R_{12})$, is exponentially attractive for a sufficiently small $\varepsilon '$ (see Theorem 3.1), then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_2(\tau ),{\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_{12}(\tau ))$ of the original system (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_1(0),{\bar{D}}^R_{12}(0))-M({\bar{D}}^A(0),{\bar{D}}^R_2(0))|$ is sufficiently small, we have a solution of the reduced system (31), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_2(\tau ))$, that satisfies (32) (see Remark 3.1). $\square$

Now, multiplying both sides of the ODEs in (31) by $D_\mathrm{{tot}}(k^A_ED_\mathrm{{tot}})$ and defining ${\bar{k}}^A_W=k^A_{W0}+k^A_W$, and ${\bar{k}}^2_W=k^2_{W0}+k^2_W$, system (31) can be rewritten as follows:

$$\begin{aligned} \dot{D}^A&=\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta +\bar{k}^R_E+k^R_ED^A)}{({\bar{k}}^A_W+k^A_MD^A)+({\bar{k}}^2_{W}+k^R_MD^R_2)}\right) D^R_2\nonumber \\&\quad -\left( \frac{({\bar{k}}^2_{W}+k^R_MD^R_2)(\delta +\bar{k}^A_E+k^A_ED^R_2)}{({\bar{k}}^A_W+k^A_MD^A)+({\bar{k}}^2_{W}+k^R_MD^R_2)}\right) D^A\\ \dot{D}^R_2&=\left( \frac{({\bar{k}}^2_{W}+k^R_MD^R_2)(\delta +\bar{k}^A_E+k^A_ED^R_2)}{({\bar{k}}^A_W+k^A_MD^A)+({\bar{k}}^2_{W}+k^R_MD^R_2)}\right) D^A\nonumber \\&\quad -\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta +\bar{k}^R_E+k^R_ED^A)}{({\bar{k}}^A_W+k^A_MD^A)+({\bar{k}}^2_{W}+k^R_MD^R_2)}\right) D^R_2.\nonumber \end{aligned}$$

(39)

This reduced system can be represented with the following chemical reactions:

$$\begin{aligned} {\mathrm{D^A} \mathop \rightarrow \limits ^{k_{AR}} \mathrm{D^R_2} },\;\;\;\mathrm{D^R_2} \mathop \rightarrow \limits ^{k_{AR}} \mathrm{D^A} \end{aligned}$$

(40)

with reaction rate coefficients $k_{AR}$ and $k_{RA}$ given by

$$\begin{aligned} k_{AR}&=\left( \frac{({\bar{k}}^2_{W}+k^R_MD^R_2)(\delta +\bar{k}^A_E+k^A_ED^R_2)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_{W}+k^R_MD^R_2}\right) ,\\ k_{RA}&=\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta +\bar{k}^R_E+k^R_ED^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^2_{W}+k^R_MD^R_2}\right) . \end{aligned}$$

It is important to point out that this reduced system includes histone modifications only, whose cooperative and competitive interactions are shown in the diagram in Fig. 1c.

4.2.2 Mathematical analysis of the stochastic properties

The stochastic behavior of the reduced chemical reaction system (40) can be represented by a one-dimensional Markov chain with state $x=n_{D^R_{2}}\in [0,{{D_\mathrm{{tot}}}}]$. Furthermore, for any state x, the rate associated with the transition to the next higher state ($x\rightarrow x+1$), $\lambda _x$, and the rate associated with the transition to the next lower state ($x\rightarrow x-1$), $\gamma _x$, for this Markov chain can be written as follows:

$$\begin{aligned} \begin{aligned} \lambda _x&=\left( \frac{({\bar{k}}^2_W+\frac{k^R_M}{\Omega }x)(\varepsilon +\varepsilon '\frac{x}{{{\textrm{D}_\mathrm{{tot}}}}})}{({\bar{u}}^A+\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}})+({\bar{u}}^R_2+\alpha \frac{x}{{{\textrm{D}_\mathrm{{tot}}}}})}\right) ({{\textrm{D}_\mathrm{{tot}}}}-x),\\ \gamma _x&=\left( \frac{({\bar{k}}^A_W+\frac{k^A_M}{\Omega }({{\textrm{D}_\mathrm{{tot}}}}-x))\mu (b\varepsilon +\varepsilon '\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}})}{({\bar{u}}^A+\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}})+({\bar{u}}^R_2+\alpha \frac{x}{{{\textrm{D}_\mathrm{{tot}}}}})}\right) x. \end{aligned} \end{aligned}$$

(41)

Let us first evaluate the stationary probability distribution $\pi (x)$. In particular, since this Markov chain is irreducible and reversible, we can exploit the expression for the stationary distribution $\pi (x)$ provided in Proposition 4.2 (Eq. 24), with transition rates $\lambda _x$ and $\gamma _x$ as defined in (41). Now, let us evaluate $\pi (x)$ for $\varepsilon \ll 1$:

Proposition 4.7

When $\varepsilon \ll 1$, the stationary distribution $\pi (x)$ associated with the one-dimensional Markov chain with rates $\lambda _x$ and $\gamma _x$ as defined in (41) can be approximated by

$$\begin{aligned} \pi _{\varepsilon \ll 1}(x)\approx {\left\{ \begin{array}{ll} \frac{1}{1+P} &{} {\text {if }} x=0 \\ 0 &{} {\text {if }} x \ne 0, {{\textrm{D}_\mathrm{{tot}}}}\\ \frac{P}{1+P} &{} {\text {if }} x = {\textrm{D}}_{\textrm{tot}}\end{array}\right. } \end{aligned}$$

(42)

with

$$\begin{aligned} P=\frac{({\bar{u}}_A+{\bar{u}}^R_2+\alpha )({\bar{u}}^R_2)}{({\bar{u}}_A+{\bar{u}}^R_2+1)({\bar{u}}_A)b}\cdot \prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\left( \frac{{\bar{u}}^R_2+\alpha \frac{i}{{{\textrm{D}_\mathrm{{tot}}}}}}{{\bar{u}}_A+\frac{{{\textrm{D}_\mathrm{{tot}}}}-i}{{{\textrm{D}_\mathrm{{tot}}}}}}\right) \cdot \left( \frac{1}{\mu }\right) ^{{{\textrm{D}_\mathrm{{tot}}}}}, \end{aligned}$$

with ${\bar{u}}_A=u_{A0}+u_A$, ${\bar{u}}^R_2=u^R_{20}+u^R_2$.

Proof

Given that for the $\lambda _x$ and $\gamma _x$ defined in (41) the product $\prod _{i=1}^{x} (\lambda _{x-1})/(\gamma _x)=O(\varepsilon )$ for any $x\ge 1$ except for $x={{\textrm{D}_\mathrm{{tot}}}}$, then the stationary probability distribution $\pi (x)$, provided in Proposition (4.2) when $\varepsilon \ll 1$ can be rewritten as done in (42). $\square$

It is possible to notice that when $\varepsilon \ll 1$, the distribution has two modes in correspondence to the fully active state $x=0$ and fully repressed state $x={{\textrm{D}_\mathrm{{tot}}}}$, and the probability of having the system in the intermediate states is approximately equal to zero. In contrast to what was observed for (25), by decreasing $\varepsilon$, P does not change; that is, the distribution does not shift toward either $x=0$ or $x={{\textrm{D}_\mathrm{{tot}}}}$. Qualitatively, when $\mu '$ is large, DNA methylation is erased quickly enough that its cooperation with the repressive histone modifications becomes less effective. This also implies that when $\mu '$ is sufficiently large, the chromatin modification circuit (Fig. 1b) can be well approximated by a circuit that takes into account histone modifications only (Fig. 1c).

Then, expression (42) also implies that when $\varepsilon \ll 1$, a system starting at $x={{\textrm{D}_\mathrm{{tot}}}}$ or at $x=0$ will tend to remain at that state, qualitatively implying that $\varepsilon$ controls the temporal extent of memory even when $\mu '$ is large. To make this statement mathematically precise, we evaluate how $\varepsilon$ affects the time to memory loss of the fully repressed chromatin state $x={{\textrm{D}_\mathrm{{tot}}}}$, $\tau _{{{\textrm{D}_\mathrm{{tot}}}}}^0=\mathbb {E}(t_{{{\textrm{D}_\mathrm{{tot}}}}}^0)$, and the time to memory loss of the fully active chromatin state $x=0$, $\tau _{0}^{{{\textrm{D}_\mathrm{{tot}}}}}=\mathbb {E}(t_{0}^{{{\textrm{D}_\mathrm{{tot}}}}})$. To this end, we can use the formulas provided in Proposition 4.4 (Eqs. 26, 27) and plug into them the transition rates defined in (41). Now, let us focus on the regime $\varepsilon \ll 1$:

Proposition 4.8

Assuming $\varepsilon '\ne 0$ and normalizing the time to memory loss with respect to $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ $({\bar{\tau }} = \tau \frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega })$, the normalized time to memory loss of the repressed and active state in the regime $\varepsilon \ll 1$ can be, respectively, approximated as follows:

$$\begin{aligned} \begin{aligned} {\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0\approx \frac{H_R}{\mu \varepsilon }\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{H^x_R}{h^x_1(\mu )}\right) ,\;\;{\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}\approx \frac{H_A}{\varepsilon }\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{h^x_2(\mu )}{H^x_A}\right) , \end{aligned} \end{aligned}$$

(43)

in which $h^x_1(\mu )$ and $h^x_2(\mu )$ are increasing functions of $\mu$ with $h^x_1(0)=0$ and $h^x_2(0)=0$, respectively, and $H_R$, $H^x_R$, $H_A$ and $H^x_A$ functions independent of $\varepsilon$ and $\mu$.

Proof

By multiplying by $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ the expressions for times to memory loss given in Prop. 4.4, with $\lambda _x$ and $\gamma _x$ defined in (41), we obtain the normalized expressions for times to memory loss. Then, by approximating them with their dominant term, that is the term of order $1/\varepsilon$ for both ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$, we obtain the expressions (43). $\square$

Both ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ are inversely proportional to $\varepsilon$. Therefore, also in this case lower $\varepsilon$ is critical to extend the temporal duration of the memory of both the active and repressed chromatin states. However, in contrast to what was observed in the previous case, both ${\bar{\tau }}_{{\textrm{D}_\mathrm{{tot}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ are $O(1/\varepsilon )$. This implies that a reduction of $\varepsilon$ has a similar effect on the memory of the repressed and active chromatin state and this is because large $\mu '$ leads to a fast erasure of DNA methylation, compared to the erasure of the other chromatin modifications, and then its cross-catalysis with the repressive histone modifications becomes less effective.

Now, let us also determine how the difference between the erasure rates of repressive and activating histone modifications, encapsulated by the non-dimensional parameter $\mu$, affect the duration of memory. From the expression for $\pi (x)$ in (42), it is possible to notice that lower $\mu$ leads to higher $\pi _{\varepsilon \ll 1}({{\textrm{D}_\mathrm{{tot}}}})$ and lower $\pi _{\varepsilon \ll 1}(0)$. This implies that the height of the peak in correspondence to the repressed state increases to the detriment of the height of the peak in correspondence to the active state. In accordance with this result, if we decrease $\mu$, then ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ increases, while ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ decreases.

4.2.3 Computational analysis

Also in this case, the analytical results were obtained using a deterministic quasi-steady-state approximation [20]. Then, in order to validate the trends obtained in Sect. 4.3.2 for the full reaction system and to extend the validity of these results to a broader parameter regime than $\varepsilon '$ sufficiently small and $\varepsilon =c\varepsilon '$, with $c=O(1)$, we conduct a computational study. In particular, we use the stochastic simulation algorithm (SSA) [12] to study via simulation the behavior of the original chemical reaction system (Fig. 1a, b) for large $\mu '$.

The effect of $\varepsilon$ and $\mu$ on the stationary distribution $\pi (x)$ (Fig. 3a) is in agreement with the results obtained by studying the analytical expression for $\pi (x)$, (24). The trend with $\varepsilon '$ is analogous to what we obtained for the previous case study (Fig. 3b). Now, let us study the effect of $\varepsilon$ on the switching time of the system temporal trajectories (Fig. 3c, d). In agreement with our analytical findings, if $\varepsilon$ is reduced, then the time that the system spends at the active state before switching to the repressed state (and vice versa) increases.

Finally, we determine via simulation the effect of $\mu$ on the reactivation time (Fig. 3e). As obtained for the previous case in which we did not consider large $\mu '$ (Fig. 2e), the time trajectories show a switch-like behavior. Furthermore, the time at which a trajectory switches to the active state after an activating input is applied is more variable for lower $\mu$.

Overall, comparing these results to the ones obtained in Sect. 4.1, it is possible to conclude that DNA methylation and its cooperation with repressive histone modifications extend the duration of memory of the repressed chromatin state.

4.3 Limiting behavior for large $\mu$

In this section, we analyze the stochastic behavior of the chromatin modification circuit for the other parameter regime of interest, that is, when $\mu$ is large (i.e., the erasure of repressive histone modification is much faster than the erasure of the other modifications). This study allows us to understand how the absence of H3K9me3 affects the stationary probability distribution and time to memory loss of chromatin states. Since we are interested in the limiting behavior for large $\mu$, we conduct again the model reduction, but now by assuming that the product $\varepsilon '\mu$ does not vanish as $\varepsilon '$ approaches zero. This leads to a different reduced model compared to the previous ones. To this end, we define the $\varepsilon '$-independent parameter $U:=\varepsilon '\mu$. Then, we introduce $U$ in the original ODE model (4), perform the model reduction with $\varepsilon '$ as a small parameter and consider the limiting case $U\rightarrow \infty$. We then conduct an analytical study to determine the stochastic behavior of the reduced system, and then a computational study to validate and extend the analytical findings.

Fig. 3

Computational analysis of the chromatin modification circuit, shown in Fig. 1b, for large $\mu '$, using SSA. a The stationary probability distribution, $\pi$, for the chromatin modification circuit represented in Fig. 1b, whose reactions are listed in Fig. 1a. The parameter values used to generate the plots are in Table S2 of S1 File. In particular, in the left-side plots $\varepsilon =0.32,0.16$, $\mu =1,0.85$ and $\varepsilon '=1$ and in the right-side plots $\varepsilon =0.32,0.16$, $\mu =1,0.85$ and $\varepsilon '=0.4$. In all plots $n_{D^A}$ and $n_{D^R}=n_{D^R_1}+n_{D^R_2}+n_{D^R_{12}}$ represent the number of nucleosomes with activating and repressive modifications. b The stationary distribution for the chromatin modification circuit for different values of $\varepsilon '$. The parameter values considered are listed in Table S2 of S1 File. In particular, $\varepsilon =0.32$ and $\varepsilon '=1,0.01$. c Time trajectories of $n_{D^A}$ and $n_{D^R}$ starting from the fully active state $n_{D^A}=50$, $n_{D^R}=0$ (left) and repressed state $n_{D^A}=0$, $n_{D^R}=n_{D^R_{12}}=50$ (right) for $\varepsilon '=1$ and different values of $\varepsilon$. d Time trajectories of $n_{D^A}$ and $n_{D^R}$, as described in c, but with $\varepsilon '=0.4$. e Time trajectories of the system starting from $n_{D^R}=n_{D^R_{12}}=50,n_{D^A}=0$ and with an input $u^A$ that, at steady state, leads to a unimodal distribution near the active state $n_{D^A}\approx 50$. Each trajectory is represented with a different color. In particular, we set $\varepsilon =0.16$, $\varepsilon '=1$, and $\mu =0.8,0.38$. In c–e, the time is normalized ($\tau =t\frac{k^A_M}{\Omega }{{\textrm{D}_\mathrm{{tot}}}}$, with $\Omega$ the reaction volume) and the parameter values are listed in Table S2 of S1 File

4.3.1 Model reduction

The result of the model reduction can be summarized by the following proposition:

Proposition 4.9

Let $\varepsilon =c\varepsilon '$, with $c=O(1)$, and let $U:=\varepsilon '\mu$ be $\varepsilon '$-independent. Then, let us consider the following system:

$$\begin{aligned} \frac{{\text {d}}\bar{D}^A}{{\text {d}}{\tau }}&=\left( \frac{\mu '(\beta \varepsilon +\varepsilon '{\bar{D}}^A)(u_A + {\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_1}\right) {\bar{D}}^R_1-\left( \frac{{\bar{u}}^R_1(\varepsilon +\varepsilon '{\bar{D}}^R_1)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_1}\right) {\bar{D}}^A\nonumber \\ \frac{{\text {d}}\bar{D}^R_1}{{\text {d}}{\tau }}&=\left( \frac{{\bar{u}}^R_1(\varepsilon +\varepsilon '{\bar{D}}^R_1)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_1}\right) {\bar{D}}^A-\left( \frac{\mu '(\beta \varepsilon +\varepsilon '{\bar{D}}^A)(u_A + {\bar{D}}^A)}{{\bar{u}}^A+{\bar{D}}^A+{\bar{u}}^R_1}\right) {\bar{D}}^R_1, \end{aligned}$$

(44)

with ${\bar{D}}^A + {\bar{D}}^R_1=1$. Furthermore, let $({\bar{D}}(\tau ),{\bar{D}}^R_2(\tau ),{\bar{D}}^R_{12}(\tau ))=M({\bar{D}}^A(\tau ),{\bar{D}}^R_1(\tau ))$ represent the system’s unique integral manifold. Then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_{12}(\tau )$, ${\bar{D}}(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}^R_2(\tau ))$ of (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_2(0),{\bar{D}}^R_{12}(0))-M({\bar{D}}^A(0),{\bar{D}}^R_1(0))|$ is sufficiently small, we have that, for $\varepsilon '$ sufficiently small, the solution of (44), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_1(\tau ))$, is such that

$$\begin{aligned} \begin{aligned} ({\bar{D}}^A(\tau ),{\bar{D}}^R_1(\tau ))&=({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_1(\tau ))+\zeta _1(\tau ),\\ ({\bar{D}}(\tau ),{\bar{D}}^R_2(\tau ),{\bar{D}}^R_{12}(\tau ))&=M({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_1(\tau ))+\zeta _2(\tau ), \end{aligned} \end{aligned}$$

(45)

in which $\zeta _i(\tau )=O(e^{-(\theta /\varepsilon ')\tau })$, with $i=1,2$, $\theta >0$, and $\tau \ge 0$.

Proof

The steps of the proof are similar to the ones in the proof of Proposition 4.6. In this case, we define x and y as $({\bar{D}}^A\;\; {\bar{D}}^R_1)^{T}$ and $({\bar{D}}^R_{12}\;\; {\bar{D}}^R_2\;\; {\bar{D}})^\mathrm{{T}}$, respectively. Then, we verify that the system is singular singularly perturbed, so that we can apply Theorem 3.1 to reduce it, and then we consider the limiting condition $U\rightarrow \infty$. Given that the integral manifold obtained, $M({\bar{D}}^A,{\bar{D}}^R_1)=(\bar{D},{\bar{D}}^R_2,{\bar{D}}^R_{12})$, is exponentially attractive for a sufficiently small $\varepsilon '$ (see Theorem 3.1), then, for any solution $({\bar{D}}^A(\tau ),{\bar{D}}^R_1(\tau ),{\bar{D}}(\tau ),{\bar{D}}^R_2(\tau ),{\bar{D}}^R_{12}(\tau ))$ of the original system (4) with initial conditions such that $|({\bar{D}}(0),{\bar{D}}^R_2(0),{\bar{D}}^R_{12}(0))-M({\bar{D}}^A(0),{\bar{D}}^R_1(0))|$ is sufficiently small, we have a solution of the reduced system (44), $({\bar{D}}^{A*}(\tau ),{\bar{D}}^{R*}_1(\tau ))$, that satisfies (45) (see Remark 3.1). See Section S.1 of S1 File for detailed derivation. $\square$

If we multiply both sides of the ODEs in (44) by $D_\mathrm{{tot}}(k^A_ED_\mathrm{{tot}})$ and define ${\bar{k}}^A_W=k^A_{W0}+k^A_W$, and ${\bar{k}}^1_W=k^1_{W0}+k^1_W$, system (44) can be rewritten as follows:

$$\begin{aligned} \dot{D}^A&=\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) D^R_1-\left( \frac{{\bar{k}}^1_W(\delta +\bar{k}^A_E+k^A_ED^R_1)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) D^A\\ \dot{D}^R_1&=\left( \frac{{\bar{k}}^1_W(\delta +\bar{k}^A_E+k^A_ED^R_1)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) D^A-\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) D^R_1, \end{aligned}$$

with ${\bar{k}}^A_W$ and ${\bar{k}}^1_W$ defined as done for the ODEs (20). This reduced system can be represented with the following chemical reactions:

$$\begin{aligned} {\mathrm{D^A} \mathop {\rightarrow }\limits ^{k_{AR}} \mathrm{D^R_1} },\;\;\;{\mathrm{D^R_1} \mathop \rightarrow \limits ^{k_{AR}} \mathrm{D^A} } \end{aligned}$$

(46)

with reaction rate coefficients defined as

$$\begin{aligned} k_{AR}&=\left( \frac{{\bar{k}}^1_W(\delta +\bar{k}^A_E+k^A_ED^R_1)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) ,\;\; k_{RA}=\left( \frac{({\bar{k}}^A_W+k^A_MD^A)(\delta ^{'}+k^{'}_T+k^{'*}_TD^A)}{{\bar{k}}^A_W+k^A_MD^A+{\bar{k}}^1_{W}}\right) . \end{aligned}$$

As opposed to what we obtained in the reduction done in Sect. 4.2.1, this system does not include repressive histone modifications, but only DNA methylation and activating histone modifications, whose interactions are shown in the diagram in Fig. 1d.

4.3.2 Mathematical analysis of the stochastic properties

The state of the one-dimensional Markov chain associated with the reduced system (46), x, represents the number of ${\textrm{D}^R_1}$, that is, $x=n_{D^R_1}\in [0,{{\textrm{D}_\mathrm{{tot}}}}]$. Furthermore, the rates associated with the transitions to the next higher and lower states, $\lambda _x$ and $\gamma _x$, respectively, can be written as

$$\begin{aligned}{} & {} \begin{aligned} \lambda _x=\left( \frac{{\bar{k}}^1_W(\varepsilon +\varepsilon '\frac{x}{{{\textrm{D}_\mathrm{{tot}}}}})}{{\bar{u}}^A+\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}}+{\bar{u}}^R_1}\right) ({{\textrm{D}_\mathrm{{tot}}}}-x), \end{aligned} \nonumber \\{} & {} \begin{aligned} \gamma _x=\left( \frac{({\bar{k}}^A_W+\frac{k^A_M}{\Omega }({{\textrm{D}_\mathrm{{tot}}}}-x))\mu '(\beta \varepsilon +\varepsilon '\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}})}{{\bar{u}}^A+\frac{({{\textrm{D}_\mathrm{{tot}}}}-x)}{{{\textrm{D}_\mathrm{{tot}}}}}+{\bar{u}}^R_1}\right) x. \end{aligned} \end{aligned}$$

(47)

Now, in order to study how large $\mu$ affects the memory of the chromatin states, we first derive the expression for the stationary probability distribution $\pi (x)$ and then the ones for the time to memory loss of the active and repressed states. Then, in the next section, we validate the theoretical predictions against stochastic simulations of the full set of chemical reactions (Fig. 1a).

Concerning the stationary distribution, also in this case we can exploit the expression for $\pi (x)$ introduced in Proposition 4.2 (Eq. 24), plugging into the transition rates $\lambda _x$ and $\gamma _x$ as defined in (47). Now, let us consider the regime $\varepsilon \ll 1$:

Proposition 4.10

When $\varepsilon \ll 1$, the stationary distribution $\pi (x)$ associated with the one-dimensional Markov chain with rates $\lambda _x$ and $\gamma _x$ as defined in (47) can be approximated by

$$\begin{aligned} \pi _{\varepsilon \ll 1}(x)\approx {\left\{ \begin{array}{ll} \frac{1}{1+P} &{} {\text {if }} x=0 \\ 0 &{} {\text {if }} x \ne 0, {{\textrm{D}_\mathrm{{tot}}}}\\ \frac{P}{1+P} &{} {\text {if }} x = {\textrm{D}}_{\textrm{tot}}\end{array}\right. } \end{aligned}$$

(48)

with

$$\begin{aligned} P=\frac{({\bar{u}}^A+{\bar{u}}^R_1)({\bar{u}}^R_1)}{({\bar{u}}^A+{\bar{u}}^R_1+1)({\bar{u}}^A)\beta }\cdot \prod _{i=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\left( \frac{{\bar{u}}^R_1}{{\bar{u}}_A+\frac{{{\textrm{D}_\mathrm{{tot}}}}-i}{{{\textrm{D}_\mathrm{{tot}}}}}}\right) \cdot \left( \frac{1}{\mu '}\right) ^{{{\textrm{D}_\mathrm{{tot}}}}}, \end{aligned}$$

with ${\bar{u}}^A=u^{A0}+u^A$, ${\bar{u}}^R_1=u^R_{10}+u^R_1$.

Proof

Given that for the $\lambda _x$ and $\gamma _x$ defined in (47) the product $\prod _{i=1}^{x} (\lambda _{x-1})/(\gamma _x)=O(\varepsilon )$ for any $x\ge 1$ except for $x={{\textrm{D}_\mathrm{{tot}}}}$, then, for the Markov chain considered here, the stationary probability distribution $\pi (x)$, provided in Proposition (4.2), can be rewritten as done in (48) when $\varepsilon \ll 1$. $\square$

From (48) it is possible to notice that, when $\varepsilon \ll 1$ the distribution is bimodal and a further reduction of $\varepsilon$ does not shift the distribution toward the repressed state. This is because considering large $\mu$ implies that repressive histone modifications (H3K9me3) are erased fast enough that their cooperation with DNA methylation becomes negligible. This confirms that when $\mu$ is sufficiently large, the chromatin modification circuit (Fig. 1b) can be well approximated by a circuit that takes into account only DNA methylation and activating histone modifications (Fig. 1d). Furthermore, comparing $\pi _{\varepsilon \ll 1}(x)$ for the large $\mu$ case, (48), with the one obtained for the large $\mu '$ case, (42), the main difference to notice is that in the large $\mu '$ case the expression for P, (43), has the $\alpha$ term, while in the large $\mu$ case the expression for P, (49), does not have it. This is because of the presence of the auto-catalytic loop for repressive histone modifications, but not for DNA methylation (Fig. 1c, d). As a consequence, when no external inputs are applied ($u^A=u^R_1=0$ and then ${\bar{u}}^A=u^A_0$, ${\bar{u}}^R_1=u^R_{10}$, with $u^A_{0}=u^R_{10}=u_0$), then the lower $u_0$, the lower P and then the more $\pi (0)$ increases to the detriment of $\pi ({{\textrm{D}_\mathrm{{tot}}}})$, that is, the distribution shifts toward the active state $x=0$.

On the contrary, in the large $\mu '$ case, even if the effect of cross-catalysis is negligible as in the large $\mu$ case, since we still have the auto-catalytic loop for $\textrm{D}^R_2$ with associated rate constant $\alpha$ (see expression for $\pi _{\varepsilon \ll 1}(x)$, (42), and P, (43)), then low $u_0$ does not have a critical effect on varying the relative height between the peaks (the values of $\pi (0)$ and $\pi ({{\textrm{D}_\mathrm{{tot}}}})$).

Now, let us evaluate how $\varepsilon$ affects $\tau _{{{_\mathrm{{tot}}}}}^0=\mathbb {E}(t_{{{\textrm{D}_\mathrm{{tot}}}}}^0)$ and $\tau _{0}^{{{\textrm{D}_\mathrm{{tot}}}}}=\mathbb {E}(t_{0}^{{{\textrm{D}_\mathrm{{tot}}}}})$, that is, the time to memory loss of the fully repressed and fully active chromatin state, respectively. To do that, we can exploit the formulas provided in Proposition 4.4 (Eqs. (26) and (27)) and plug into them the transition rates defined in (47). Now, focusing on the regime $\varepsilon \ll 1$, these expressions can be approximated as shown in the following proposition:

Proposition 4.11

Assuming $\varepsilon '\ne 0$ and normalizing the time to memory loss with respect to $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ (${\bar{\tau }} = \tau \frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$), the normalized time to memory loss of the repressed and active state in the regime $\varepsilon \ll 1$ can be, respectively, approximated as follows:

$$\begin{aligned} \begin{aligned} {\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0\approx \frac{{\tilde{L}}_R}{\mu '\varepsilon }\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{{\tilde{L}}^x_R}{{\tilde{l}}^x_1(\mu ')}\right) ,\;\;{\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}\approx \frac{{\tilde{L}}_A}{\varepsilon }\left( 1+\sum _{x=1}^{{{\textrm{D}_\mathrm{{tot}}}}-1}\frac{{\tilde{l}}^x_2(\mu ')}{{\tilde{L}}^x_A}\right) , \end{aligned} \end{aligned}$$

(49)

in which $l^x_1(\mu ')$ and $l^x_2(\mu ')$ are increasing functions of $\mu '$ with $l^x_1(0)=0$ and $l^x_2(0)=0$, respectively, and $L_R$, $L^x_R$, $L_A$ and $L^x_A$ functions independent of $\varepsilon$ and $\mu '$.

Proof

By multiplying by $\frac{k^A_M{{\textrm{D}_\mathrm{{tot}}}}}{\Omega }$ the expressions for times to memory loss given in Prop. 4.4, with $\lambda _x$ and $\gamma _x$ defined in (47), we obtain the normalized expressions for times to memory loss. Then, by approximating them with their dominant term, that is the term of order $1/\varepsilon$ for both ${\bar{\tau }}_{{{\textrm{D}_\mathrm{{tot}}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$, we obtain the expressions (49). $\square$

As for the large $\mu '$ case, both ${\bar{\tau }}_{{\textrm{D}_\mathrm{{tot}}}}^0$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}$ are $O(1/\varepsilon )$. This implies that lower $\varepsilon$ extends in a similar way the memory of both the active and repressed chromatin states, in contrast to what was observed for the original study case in Sect. 4.1, in which ${\bar{\tau }}_{{\textrm{D}_\mathrm{{tot}}}}^0=O(1/\varepsilon ^2)$ and ${\bar{\tau }}_{0}^{{\textrm{D}_\mathrm{{tot}}}}=O(1/\varepsilon )$. Concerning the effect of $\mu '$ (the non-dimensional parameter encapsulating the asymmetry between the erasure rates of DNA methylation and activating histone modifications) on the memory of the chromatin states, it is possible to notice that its trend on the stationary distribution and time to memory loss is the same as the one that $\mu$ has in the $\mu '$ case study (Sect. 4.2.2).

Fig. 4

Computational analysis of the chromatin modification circuit, shown in Fig. 1b, for large $\mu$, using SSA. a The stationary probability distribution, $\pi$, for the chromatin modification circuit represented in Fig. 1b, whose reactions are listed in Fig. 1a. The parameter values used to generate the plots are in Table S3 of S1 File. In particular, in the left-side plots $\varepsilon =0.16,0.08$, $\mu '=1,0.185$ and $\varepsilon '=1$ and in the right-side plots $\varepsilon =0.1,0.06$, $\mu '=1,0.15$ and $\varepsilon '=0.4$. In all plots $n_{D^A}$ and $n_{D^R}=n_{D^R_1}+n_{D^R_2}+n_{D^R_{12}}$ represent the number of nucleosomes with activating and repressive modifications. b The stationary distribution for the chromatin modification circuit for different values of $\varepsilon '$. The parameter values considered are listed in Table S3 of S1 File. In particular, $\varepsilon =0.16$ and $\varepsilon '=1,0.01$. c Time trajectories of $n_{D^A}$ and $n_{D^R}$ starting from the fully active state $n_{D^A}=50$, $n_{D^R}=0$ (left) and repressed state $n_{D^A}=0$, $n_{D^R}=n_{D^R_{12}}=50$ (right) for $\varepsilon '=1$ and different values of $\varepsilon$. d Time trajectories of $n_{D^A}$ and $n_{D^R}$, as described in c, but with $\varepsilon '=0.4$. e Time trajectories of the system starting from $n_{D^R}=n_{D^R_{12}}=50,n_{D^A}=0$ and with an input $u^A$ that, at steady state, leads to a unimodal distribution near the active state $n_{D^A}\approx 50$. Each trajectory is represented with a different color. In particular, we set $\varepsilon =0.16$, $\varepsilon '=1$, and $\mu '=0.8,0.38,0.08$. In c–e, the time is normalized ($\tau =t\frac{k^A_M}{\Omega }{{\textrm{D}_\mathrm{{tot}}}}$, with $\Omega$ the reaction volume) and the parameter values are listed in Table S3 of S1 File

4.3.3 Computational analysis

We use the stochastic simulation algorithm (SSA) [12] to study via simulation the original chemical reaction system (Fig. 1a, b) for large $\mu$. We can first notice that the trend with which $\varepsilon$ and $\mu '$ affect the stationary distribution $\pi (x)$ is in agreement with the analytical findings (Fig. 4a). That is, smaller $\varepsilon$ leads to more concentrated peaks, and reducing $\mu '$ increases the height of the peak for the repressed state to the detriment of the height of the active state peak. It is important to point out that in this case, when $\mu '=1$ (DNA methylation and activating histone modifications have the same erasure rate), the distribution is shifted toward the active state (Fig. 4a). This bias is given by the presence of the auto-catalytic loop characterizing the histone modification dynamics, but not the DNA methylation dynamics (Fig. 1c, d). Furthermore, the effect of $\varepsilon '$ is similar to what was observed for the previous case studies (Fig. 4b). We then consider a parameter regime in which the system displays a bimodal distribution and study the effect of $\varepsilon$ on the switching time of the temporal trajectories (Fig. 4c, d). It is possible to notice that smaller values of $\varepsilon$ increase the time that the system spends at the active state before switching to the repressed state, and vice versa. These results are in agreement with the ones obtained by studying the analytical expression for the time to memory loss of the repressed and active state (49).

Finally, concerning the reactivation time of this system (Fig. 4e), it is possible to notice that the absence of repressive histone modifications, compared to the case in which we do not have DNA methylation (large $\mu '$ case), leads to shorter reactivation time, unless $\mu '$ is sufficiently small. However, even for lower $\mu '$ and then slower reactivation, the time needed to switch to the active state is less variable compared to the previous case studies, suggesting that large $\mu$, i.e., the absence of repressive histone modifications, could reduce the stochasticity of gene reactivation.

Overall, similarly to what was obtained in the previous section, this analysis shows that when repressive histone modifications are erased quickly enough that their cooperation with DNA methylation becomes less effective, the duration of the repressed chromatin state memory decreases. However, in contrast to what we obtained for the previous limiting case study, the reactivation process of the system characterized by large $\mu$ is less stochastic compared to the reactivation process of the full system (Fig. 2e). This result highlights that repressive histone modifications contribute to the highly stochastic latency of state reactivation.

5 Discussion and conclusion

In this work, we considered a chromatin modification circuit including both histone modifications and DNA methylation [8] (Fig. 1b) to single out the specific contributions of DNA methylation and histone modifications to the duration of the active and repressed chromatin states memory. For this purpose, we first proved that system (3), with $\varepsilon '$ as a small parameter, is singular singularly perturbed and then exploited a proper reduction approach proposed in [9] to obtain a one-dimensional model suitable for analytical study. We performed this model reduction for the full chromatin modification system (Sect. 4.1) and for two limiting cases: DNA methylation almost completely absent (Sect. 4.2) and repressive histone modifications almost completely absent (Sect. 4.3).

Our analysis showed that the coexistence and interaction between DNA methylation and repressive histone modifications biases the system stationary distribution toward the repressed stated and, accordingly, strengthens memory of the repressed chromatin state (Fig. 2). When $\mu '$ is large enough to have a negligible amount of DNA methylation in the system and then the interplay between repressive chromatin marks does not have a relevant effect on the system dynamics (Fig. 1c), the bias in the stationary distribution and the asymmetry between the active and repressed state memory are reduced (Fig. 3a, c, d). However, the latency of state reactivation remains highly stochastic, especially for low $\mu$ (Fig. 3e). Different results can be obtained when $\mu$ is large (repressive histone modifications almost completely absent). The reason is that not only the cooperative interactions among repressive modifications can be neglected, but here the remaining repressive mark (DNA methylation) does not have the positive feedback loop associated with the auto-catalytic process (Fig. 1d). This implies that the state reactivation latency becomes less stochastic (Fig. 4e).

These results suggest then the removal of repressive histone modifications and of the positive feedback loops associated with the auto- and cross-catalysis could reduce the stochasticity associated with the reactivation of a silenced gene, shown in earlier experimental studies [21]. As future work, we are planning to develop theoretical tools to derive analytical expressions for stationary distribution and times to memory loss for our original reaction system and then to obtain a quantitative characterization of the original system. Furthermore, we will also conduct experimental investigations in order to validate the theoretical results obtained in this paper.