Mathematical analysis of the limiting behaviors of a chromatin modification circuit

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.


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 onedimensional 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].
The basic unit of the model is D, that is, the nucleosome with DNA wrapped around it. Then, we have D A , that is, nucleosome with H3K4me3/ac, D R 2 , that is, nucleosome with H3K9me3, D R 1 , that is, nucleosome with CpGme, and 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 0 , 1 , and 8 ). Then, the read-write mechanism, in which histone modifications recruit marks 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) of the same kind to nearby nucleosomes, generates auto-catalytic loops (reactions 2 , 3 ). 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 12 , 13 ). Finally, basal erasure and recruited erasure, wherein activating marks recruit repressive mark's eraser enzymes and vice versa, are represented by reactions 4 , 5 , 9 and 6 , 7 , 10 , 11 , 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 tot andD = n D /D tot , with D 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 D 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 = D tot / , with the reaction volume, and let us define the normalized inputs asū R We considerū R 1 ,ū R 2 andū 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 α = k M /k A M ,ᾱ =k M /k A M and α = 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 More precisely, μ represents the ratio between the erasure rates of repressive histone modifications and activating histone modifications and Table 1 Definitions and interpretations of ε, ε , μ, and μ Param.
Definition Interpretation Parameter that scales the ratio between the rate of the basal erasure and the one of the Auto-/cross-catalytic loop of each mark Parameter that scales the ratio between the rate of the recruited erasure and the one of the Auto-/cross-catalytic loop of each mark Ratio between the erasure rates of repressive and activating histone modifications Ratio between the erasure rates of DNA methylation and activating histone modifications μ represents the ratio between the erasure rates of DNA methylation and activating histone modifications. Furthermore, based on the previous definitions, we have that This implies that the dimensionless parameter ε 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 R E /k A M = με and k * T /k A M = μ ε , the dimensionless parameter ε 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 τ = tk A M D tot , the ODEs associated with the full chromatin modification circuit are withD +D A +D R 1 +D R 2 +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 ε = cε , with c = O (1). Then, introducing the new time variableτ = τ ε in the ODEs (3), the system of equations can be rewritten as follows: Here, ε 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 μ , we define the parameter U := ε μ and assume that it is ε -independent, so that the product ε μ does not vanish as ε approaches zero. Similarly, in order to study the behavior of the system for large μ, we will define U := ε μ and assume that it is ε -independent, so that ε μ does not vanish as ε approaches zero. Each case will lead to a different reduced model, as shown in Sect. 4.

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 dx dt = f (x, y, t) with x ∈ R n , t ∈ R, let us define a smooth surface S in R n × 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 ∈ R m and y 2 ∈ R n , let us introduce the system with functions f 1 and f 2 sufficiently smooth, and let us define the matrix If A(x, y 2 , t, 0) is singular on some subspace of R m × R n × 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 = φ(x, t) with t ∈ R and x ∈ R m ; • C2: the matrix A, defined in (6), with y 2 = φ(x, t) and ε = 0 has a kernel of dimension m and m corresponding linearly independent eigenvectors, and the matrix has n eigenvalues λ i (x, t) : Re(λ i ) ≤ −2θ , with θ > 0; • C3: defining the domain X as X = {(x, y 2 , t, ε )|x ∈ R m , ||y 2 − φ(x, t)|| ≤ ρ, t ∈ R, 0 ≤ ε ≤ ε 0 }, the functions f 1 , f 2 and the matrix A are continuously differentiable (k + 2) times in X , with k ≥ 0 for some positive ε 0 and ρ.

Remark 3.1
Given the fact that the integral manifold is exponentially attractive for a sufficiently small ε , then, for any solution (x(t), y 1 (t)) of (8) with initial conditions | is sufficiently small, we have a solution of the reduced system (10) such that 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).

Results
In this section, we obtain reduced versions of the full chromatin modification system (4) in the limit where ε approaches zero. We first show that, considering ε 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 ε -independent parameter U := ε μ , so that ε μ does not vanish as ε approaches zero in the model reduction, and then we consider the limiting behavior as U → ∞. Similarly, in order to single out the contribution of DNA methylation, we repeat the model reduction introducing the εindependent parameter U := ε μ, so that ε μ does not vanish as ε approaches zero, and then we consider the limiting behavior as U → ∞.

Behavior of the full chromatin modification circuit with " as small parameter
Let us summarize the results of the model reduction in the following proposition: 12 (0))| is sufficiently small, we have that, for ε sufficiently small, the solution of (13), (D A * (τ ),D R * 12 (τ )), is such that Proof Let us consider the ODE model (4), and let us define x, y 2 , f 1 and f 2 as follows: Now, it is possible to show that φ, defined in C1, is equal to φ(x) = (0, 0, 0), and that matrix A, defined in (6), withD =D R 1 =D R 2 = 0 and ε = 0 can be written as follows: , , 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 =Ā 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 This implies that we can apply Theorem 3.1 to reduce our system. To this end, let us first introduce the new variablesD =D/ε , (4): Now, to determine the integral manifold M( , let us find the expression for the asymptotic expansion ofD,D R 1 andD R 2 : To this end, let us substitute (17) in the first three equations of (16) to obtain 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 ε . More precisely, given that ∂h i0 1 for sufficiently small values of ε , we can write h i0 and h i1 , with i = 0, 1, 2, as follows: Then, by introducing in the last two equations of (16) the asymptotic expansion of D,D R 1 andD 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 τ =τ /ε . The sum of the equations in (13) is equal to zero, implying that can be approximately set equal to 1 for sufficiently small ε . Furthermore, given that the integral manifold obtained, , is exponentially attractive for a sufficiently small ε (see Theorem 3.1), then, for any solution Now, if we multiply both sides of the ODEs in (13) , we can rewrite system (13) as follows: This reduced system can be represented by the following chemical reactions: with reaction rate coefficients k AR and k R A given by

Mathematical analysis of the stochastic properties
Since the reduced chemical reaction system (21) is characterized by the conservation law D R 12 + D A ≈ D tot , its stochastic behavior can be approximately represented by a one-dimensional Markov chain with state x = n D R 12 ∈ [0, D tot ]. For any state x, the rate associated with the transition to the next higher state (x → x + 1), λ x , and the rate associated with the transition to the next lower state (x → x − 1), γ x , can be written as follows: 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 λ x and γ x represent the rate associated with the transition x → x + 1and the rate associated with the transition x → x − 1, respectively. Then, the stationary distribution associated with the one-dimensional Markov chain with rates λ x and γ x can be written as in which D tot x=0 π(x) = 1. Proof By using detailed balance, we can write λ x−1 π(x − 1) = γ x π(x), for any x ∈ [1, D tot ]. Then, these equalities can be combined to write π(x) = x i=1 (λ i−1 /γ i )π(0). Finally, exploiting D tot x=0 π(x) = 1, we obtain the formula in (24). Proposition 4.3 When ε 1, the stationary distribution π(x) associated with the one-dimensional Markov chain with rates λ x and γ x as defined in (23) can be written as with P given by , .
By studying the expression for π ε 1 (x) in (25), it is possible to notice that if ε 1, the only states in which π(x) does not vanish are the fully active state x = 0 and the fully repressed state x = D tot . More precisely, the stationary distribution for ε 1 is bimodal, with two modes in correspondence to x = 0 and x = D tot , and the probability of having the system in one of the intermediate states is approximately zero. Furthermore, when ε decreases, P increases and accordingly π ε 1 (D tot ) increases to the detriment of π ε 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 ε 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 ε 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: where x(t) is the Markov chain described above. The time to memory loss of the fully repressed chromatin state is defined as τ 0 D tot = E(t 0 D tot ). Similarly, the time to memory loss of the active state is defined as τ D tot 0 = E(t D tot 0 ).

Proposition 4.4 The time to memory loss of the repressed chromatin state is given by
in which s x = λ 1 λ 2 ...λ x γ 1 γ 2 ...γ x and λ x and γ x are defined in (23). The time to memory loss of the active chromatin state is given by in Proof By using first step analysis [11], we can write the following equations: Then, by solving system (28), we obtain the formula for τ 0 D tot as in (26). A similar approach can be used to obtain the formula for τ D tot 0 as in (27).
in which g x 1 (μμ ) and g x 2 (μμ ) are increasing functions of μμ with g x 1 (0) = 0 and g x 2 (0) = 0, respectively, and G x R , G R , G A and G x A are functions that do not depend on ε, μ and μ. = O(1/ε). This implies that decreasing ε 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 nondimensional parameters μ and μ . From the expression for the stationary distribution in (25) it is possible to notice that, by reducing μ or μ (i.e., reducing the erasure rates of the repressive marks compared to the erasure rate of the active marks), π ε 1 (D tot ) increases to the detriment of π ε 1 (0), i.e., the stationary distribution shifts toward the repressed state. In agreement with these results, when μ or μ decreases,τ 0 D tot increases, whileτ D tot 0 decreases, that is, the temporal duration of the memory of the repressed state increases, while the duration of memory of the active state decreases.

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 ε 1 and ε = cε , 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 ε and μ affect the stationary distribution of the original system π(x) is in agreement with the results obtained from the analytical study in Sect. 4.1.1. The parameter ε does not significantly vary the way in which ε, μ , μ affect the stationary distribution. However, decreasing ε compared to ε leads to less concentrated peaks in the bimodal stationary distribution and, by further decreasing ε , 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 ε (Fig. 2c, d).
In particular, in agreement with our analytical findings, it is possible to notice that lowering ε 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 μ 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 μ . This implies that for low values of μ the reactivation of the gene becomes a very stochastic process.

Limiting behavior for large
In this section, we determine the stochastic behavior of the chromatin modification circuit when μ 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 Fig. 2 Computational analysis of the full chromatin modification circuit, shown in Fig. 1b, using SSA. a The stationary probability distribution, π , 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 ε = 0.28, 0.14, μ = 0.675, 0.35 and ε = 1 and in the right-side plots ε = 0.28, 0.14, μ = 0.625, 0.35 and ε = 0.4. In all plots n D A and n D R = n D R  Table S1 of S1 File. In particular, ε = 0.28 and ε = 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  Table S1 of S1 File 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 μ , we conduct again the model reduction by assuming that the product ε μ does not vanish as ε 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 := ε μ and assume that it is ε -independent. Then, we introduce U in the ODE model (4), perform the model reduction with ε as a small parameter and consider the limiting case U → ∞. Then, we conduct an analytical study of the stochastic behavior of the reduced system, and validate and extend the results obtained via simulation.

the system's unique integral manifold. Then, for any solution (D
| is sufficiently small, we have that, for ε sufficiently small, the solution of (31), (D A * (τ ),D R * 2 (τ )), is such that Proof Let us introduce U = ε μ in system (4), obtaining Now, let us define x, y 2 , f 1 and f 2 as Now, it is possible to calculate that φ(x) = (0, 0, φ 12 ), with φ(x) defined in C1 and Function φ 12 is inversely proportional to U . Furthermore, the matrix A, defined in (6), withD =D R 1 = 0,D R 12 = φ 12 and ε = 0 can be written as with ,Ã 1,2 = ( −(u R 20 +α(D R 2 +φ 12 )+ᾱφ 12 +U (βc+D A )) −(u R 10 +α (D R 2 +φ 12 )) ) , 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 =Ā 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 , matrix B has three eigenvalues with negative real part if u R 10 , u R 20 , u A 0 ≥ 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 variableD R 12 =D R 12 − φ 12 , and then the variablesD =D/ε ,D R Now, similarly to what we did in Sect. 4.1, in order to determine the integral manifold M(D A ,D R 2 ) = (D,D R 1 ,D R 12 ), we evaluate the asymptotic expansion ofD, D R 1 andD R 12 , that can be written as follows: To this end, we plug (36) into the first three equations of (35), obtaining 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 ε . Specifically, since ∂h i0 ∂D R 2 and ∂h i0 ∂D A are bounded for any i = 0, 1, 2 (i.e., ε ∂h i0 1 for sufficiently small ε ) and since for U 1 we have that φ 12 1 and dφ 12 /dτ ≈ (u R 10 + 2α D R 2 )h 2 , we can rewrite the expressions for h i0 and h i1 , i = 0, 1, 2, as follows: , .
Then, by considering the limiting condition U → ∞, the expressions for h i0 and h i1 can be approximated as Now, by plugging the asymptotic expansion ofD,D R 1 andD 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 τ =τ /ε . It is possible to notice that the sum of the two ODEs in (31) is equal to zero, implying thatD A +D R 2 = constant. Since the conservation law D A +D R 12 +D +D R 1 +D R 2 = 1 holds and, for sufficiently small ε and sufficiently can be approximately set equal to 1 for sufficiently small values of ε and sufficiently large U . Furthermore, given that the integral manifold obtained, , is exponentially attractive for a sufficiently small ε (see Theorem 3.1), then, for any solution )| is sufficiently small, we have a solution of the reduced system (31), (D A * (τ ),D R * 2 (τ )), that satisfies (32) (see Remark 3.1). Now, multiplying both sides of the ODEs in (31) by D tot (k A E D tot ) and defininḡ k A W = k A W 0 + k A W , andk 2 W = k 2 W 0 + k 2 W , system (31) can be rewritten as follows: This reduced system can be represented with the following chemical reactions: with reaction rate coefficients k AR and k R A given by 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.

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 ∈ [0, D tot ]. Furthermore, for any state x, the rate associated with the transition to the next higher state (x → x + 1), λ x , and the rate associated with the transition to the next lower state (x → x − 1), γ x , for this Markov chain can be written as follows: Let us first evaluate the stationary probability distribution π(x). In particular, since this Markov chain is irreducible and reversible, we can exploit the expression for the stationary distribution π(x) provided in Proposition 4.2 (Eq. 24), with transition rates λ x and γ x as defined in (41). Now, let us evaluate π(x) for ε 1: Proposition 4.7 When ε 1, the stationary distribution π(x) associated with the onedimensional Markov chain with rates λ x and γ x as defined in (41) can be approximated by Proof Given that for the λ x and γ x defined in (41) the product x i=1 (λ x−1 )/(γ x ) = O(ε) for any x ≥ 1 except for x = D tot , then the stationary probability distribution π(x), provided in Proposition (4.2) when ε 1 can be rewritten as done in (42).
It is possible to notice that when ε 1, the distribution has two modes in correspondence to the fully active state x = 0 and fully repressed state x = D 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 ε, P does not change; that is, the distribution does not shift toward either x = 0 or x = D tot . Qualitatively, when μ is large, DNA methylation is erased quickly enough that its cooperation with the repressive histone modifications becomes less effective. This also implies that when μ 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 ε 1, a system starting at x = D tot or at x = 0 will tend to remain at that state, qualitatively implying that ε controls the temporal extent of memory even when μ is large. To make this statement mathematically precise, we evaluate how ε affects the time to memory loss of the fully repressed chromatin state x = D tot , τ 0 D tot = E(t 0 D tot ), and the time to memory loss of the fully active chromatin state x = 0, τ D tot 0 = E(t D tot 0 ). 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 ε 1: 1 can be, respectively, approximated as follows:  (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/ε for bothτ 0 D tot andτ D tot 0 , we obtain the expressions (43).
Bothτ 0 D tot andτ D tot 0 are inversely proportional to ε. Therefore, also in this case lower ε 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τ 0 D tot andτ D tot 0 are O(1/ε). This implies that a reduction of ε has a similar effect on the memory of the repressed and active chromatin state and this is because large μ 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 μ, affect the duration of memory. From the expression for π(x) in (42), it is possible to notice that lower μ leads to higher π ε 1 (D tot ) and lower π ε 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 μ, thenτ 0 D tot increases, whileτ D tot 0 decreases.

Computational analysis
Also in this case, the analytical results were obtained using a deterministic quasisteady-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 ε sufficiently small and ε = cε , 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 μ . The effect of ε and μ on the stationary distribution π(x) (Fig. 3a) is in agreement with the results obtained by studying the analytical expression for π(x), (24). The trend with ε is analogous to what we obtained for the previous case study (Fig. 3b). Now, let us study the effect of ε on the switching time of the system temporal trajectories (Fig. 3c, d). In agreement with our analytical findings, if ε 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 μ on the reactivation time (Fig. 3e). As obtained for the previous case in which we did not consider large μ (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 μ.
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.

Limiting behavior for large
In this section, we analyze the stochastic behavior of the chromatin modification circuit for the other parameter regime of interest, that is, when μ 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 μ, we conduct again the model reduction, but now by assuming that the product ε μ does not vanish as ε approaches zero. This leads to a different reduced model compared to the previous ones. To this Fig. 3 Computational analysis of the chromatin modification circuit, shown in Fig. 1b, for large μ , using SSA. a The stationary probability distribution, π , 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 ε = 0.32, 0. 16 Table  S2 of S1 File end, we define the ε -independent parameter U := ε μ. Then, we introduce U in the original ODE model (4), perform the model reduction with ε as a small parameter and consider the limiting case U → ∞. 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.

Model reduction
The result of the model reduction can be summarized by the following proposition: Proposition 4.9 Let ε = cε , with c = O(1), and let U := ε μ be ε -independent. Then, let us consider the following system: | is sufficiently small, we have that, for ε sufficiently small, the solution of (44), (D A * (τ ),D R * 1 (τ )), is such that 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 (D AD R 1 ) T and (D R 12D R 2D ) 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 → ∞. If we multiply both sides of the ODEs in (44) by (44) can be rewritten as follows: withk A W andk 1 W defined as done for the ODEs (20). This reduced system can be represented with the following chemical reactions: with reaction rate coefficients defined as 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.

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 D R 1 , that is, . Furthermore, the rates associated with the transitions to the next higher and lower states, λ x and γ x , respectively, can be written as Now, in order to study how large μ affects the memory of the chromatin states, we first derive the expression for the stationary probability distribution π(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 π(x) introduced in Proposition 4.2 (Eq. 24), plugging into the transition rates λ x and γ x as defined in (47). Now, let us consider the regime ε 1: Proposition 4.10 When ε 1, the stationary distribution π(x) associated with the one-dimensional Markov chain with rates λ x and γ x as defined in (47) can be approximated by Given that for the λ x and γ x defined in (47) the product x i=1 (λ x−1 )/(γ x ) = O(ε) for any x ≥ 1 except for x = D tot , then, for the Markov chain considered here, the stationary probability distribution π(x), provided in Proposition (4.2), can be rewritten as done in (48) when ε 1.
From (48) it is possible to notice that, when ε 1 the distribution is bimodal and a further reduction of ε does not shift the distribution toward the repressed state. This is because considering large μ implies that repressive histone modifications (H3K9me3) are erased fast enough that their cooperation with DNA methylation becomes negligible. This confirms that when μ 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 π ε 1 (x) for the large μ case, (48), with the one obtained for the large μ case, (42), the main difference to notice is that in the large μ case the expression for P, (43), has the α term, while in the large μ 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ū A = u A 0 ,ū 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 π(0) increases to the detriment of π(D tot ), that is, the distribution shifts toward the active state x = 0.
On the contrary, in the large μ case, even if the effect of cross-catalysis is negligible as in the large μ case, since we still have the auto-catalytic loop for D R 2 with associated rate constant α (see expression for π ε 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 π(0) and π(D tot )). Now, let us evaluate how ε affects τ 0 tot = E(t 0 D tot ) and τ D tot 0 = E(t D tot 0 ), 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 ε 1, these expressions can be approximated as shown in the following proposition: 1 can be, respectively, approximated as follows: in which l x 1 (μ ) and l x 2 (μ ) are increasing functions of μ 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 ε and μ .
Proof By multiplying by k A M D tot the expressions for times to memory loss given in Prop. 4.4, with λ x and γ 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/ε for bothτ 0 D tot andτ D tot 0 , we obtain the expressions (49).

Fig. 4
Computational analysis of the chromatin modification circuit, shown in Fig. 1b, for large μ, using SSA. a The stationary probability distribution, π , 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 ε = 0. 16 Table S3 of S1 File. In particular, ε = 0. 16  Concerning the effect of μ (the nondimensional 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 μ has in the μ case study (Sect. 4.2.2).

Computational analysis
We use the stochastic simulation algorithm (SSA) [12] to study via simulation the original chemical reaction system (Fig. 1a, b) for large μ. We can first notice that the trend with which ε and μ affect the stationary distribution π(x) is in agreement with the analytical findings (Fig. 4a). That is, smaller ε leads to more concentrated peaks, and reducing μ 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 μ = 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 ε 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 ε on the switching time of the temporal trajectories (Fig. 4c, d). It is possible to notice that smaller values of ε 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 μ case), leads to shorter reactivation time, unless μ is sufficiently small. However, even for lower μ 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 μ, 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 μ 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.

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 ε 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 μ 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 μ (Fig. 3e). Different results can be obtained when μ 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.