Appendix. Stochastic method details

To simulate the stochastic dynamics of the whole system, the Gillespie algorithm has been used. Let us consider N = 40 histones and 2N reaction channels in total. Here, in the simulation, let us represent the state vector of histones at time t by s(t) =  and propensity functions for i-th histone, p2i − 1(t) for reaction between si(t) and si(t) + 1 and p2i(t) for reaction between si(t) and si(t) − 1 defined as

$$_=\left[_4\textrm}f\left(_i\right)+_4\textrm}\right]\left(_+_+_\right)+\left[_27\textrm}f\left(_i\right)+_4\textrm}\right]\left(_+_+_\right)$$

$$_=\left[_27\textrm}g\left(_i\right)+_27\textrm}\right]\left(_+_+_\right)+\left[_4\textrm}g\left(_i\right)+_4\textrm}\right]\left(_+_+_\right),$$

where j = 1, 2, . . , N, \(_i=\frac}\sum_\ \textrm\ i}_\), \(_i=\frac}\sum_\ \textrm\ i}_\), nnei is the number of histone neighbors for the global interaction that has been chosen N-1 at this study, \(f\left(_i\right)=\frac_f}\) and \(g\left(_i\right)=\frac_g}\), and \(_f=_g=\frac\), such that αk4me3 = γk4me3kk4me3, αk27me3 = γk27me3kk27me3, βk4me3 = γk4me3dk4me3, and βk27me3 = γk27me3dk27me3 where γk4me3 and γk27me3 are level of noise-to-feedback ratio. The j-th reaction channel induces a change on the state of the whole system which can be represented by a vector φj, j = 1, 2, …N such that if j is odd, then the \(\frac\)-th element of φj is 1 and the other N-1 elements are zeros and if j is even, then the \(\frac\)-th element of φjis -1, and the other N-1 elements are zeros.

The Gillespie algorithm is as follows:

1)

Initialization. Set the desired values to S (0) = 3.

2)

Calculate and store the values of 2N propensity functions pi(t), i = 1, 2, ..2N at current time t. Also, calculate the sum of the 2N values, denoted as p0.

3)

Generate two independent random numbers with uniform distribution over the interval [0, 1], r1and r2, and calculate the waiting time \(\tau =-\frac_1\right)}\) for the next reaction to occur. Determine the integer value ϑ as an index of reaction

Channel, such that the relation \(\sum_^_j(t)\le _2_o\le \sum_^_j(t)\) is satisfied.

4)

Update the time variable as τ + t and update the system state as s(t + τ) = s(t) + φϑ.

5)

Repeat the steps 2-4 until the terminating condition is satisfied. DNA replication during the cell cycle L=20 is also simulated in the study. If (t + τ) exceeded the time of DNA replication, i.e. t < nL ≤ t + τ, for some positive integer n, then we will stop the Gillespie algorithm and execute DNA replication. In particular, time variable t is updated as nL rather than t + τ, and each nucleosome will be replaced by an unmodified nucleosome (i.e., both s2k − 1(nL) and s2k(nL) ) are set to be zeros, for \(k=1,,2,..\frac\) with a probability of 0.5. After DNA replication, we need to resume by repeating the steps 2-4 of the Gillespie algorithm.