A multiscale model of epigenetic heterogeneity reveals the kinetic routes of pathological cell fate reprogramming

The inherent capacity of somatic cells to switch their phenotypic status in response to damage stimuli in vivo might have a pivotal role in ageing and cancer. However, how the entry-exit mechanisms of phenotype reprogramming are established remains poorly understood. In an attempt to elucidate such mechanisms, we herein introduce a stochastic model of combined epigenetic regulation (ER)-gene regulatory network (GRN) to study the plastic phenotypic behaviours driven by ER heterogeneity. Furthermore, based on the existence of multiple scales, we formulate a method for stochastic model reduction, from which we derive an efficient hybrid simulation scheme that allows us to deal with such complex systems. Our analysis of the coupled system reveals a regime of tristability in which pluripotent stem-like and differentiated steady-states coexist with a third indecisive state. Crucially, ER heterogeneity of differentiation genes is for the most part responsible for conferring abnormal robustness to pluripotent stem-like states. We then formulate epigenetic heterogeneity-based strategies capable of unlocking and facilitating the transit from differentiation-refractory (pluripotent stem-like) to differentiation-primed epistates. The application of the hybrid numerical method validated the likelihood of such switching involving solely kinetic changes in epigenetic factors. Our results suggest that epigenetic heterogeneity regulates the mechanisms and kinetics of phenotypic robustness of cell fate reprogramming. The occurrence of tunable switches capable of modifying the nature of cell fate reprogramming from pathological to physiological might pave the way for new therapeutic strategies to regulate reparative reprogramming in ageing and cancer. Author summary Certain modifications of the structure and functioning of the protein/DNA complex called chromatin can allow adult, fully differentiated cells to adopt a stem cell-like pluripotent state in a purely epigenetic manner, not involving changes in the underlying DNA sequence. Such reprogramming-like phenomena may constitute an innate reparative route through which human tissues respond to injury and could also serve as a novel regenerative strategy in human pathological situations in which tissue or organ repair is impaired. However, it should be noted that in vivo reprogramming would be capable of maintaining tissue homeostasis provided the acquisition of pluripotency features is strictly transient and accompanied by an accurate replenishment of the specific cell types being lost. Crucially, an excessive reprogramming to pluripotency in the absence of controlled re-differentiation would impair the repair or the replacement of damaged cells, thereby promoting pathological alterations of cell fate. A mechanistic understanding of how the degree of chromatin plasticity dictates the reparative versus pathological behaviour of in vivo reprogramming to rejuvenate aged tissues while preventing tumorigenesis is urgently needed, including especially the intrinsic epigenetic heterogeneity of the tissue resident cells being reprogrammed. We here introduce a novel method that mathematically captures how epigenetic heterogeneity is actually the driving force that governs the routes and kinetics to entry into and exit from a pathological pluripotent-like state. Moreover, our approach computationally validates the likelihood of unlocking chronic, unrestrained pluripotent states and drive their differentiation down the correct path by solely manipulating the intensity and direction of few epigenetic control switches. Our approach could inspire new therapeutic approaches based on in vivo cell reprogramming for efficient tissue regeneration and rejuvenation and cancer treatment.


Introduction
The correlation between ageing and cancer incidence rate is a well established empirical 26 fact. The currently accepted explanation for such a correlation is subsumed under the 27 multiple hit hypothesis or Knudson hypothesis [1,2]. This multi-mutation theory 28 considers cancer as the result of the successive accumulation of genetic mutations and, 29 therefore, time is necessary for cells to overcome a certain mutagenic threshold before 30 cancer can develop. An alternative paradigm has emerged in recent years, in which the 31 ability of the ageing process per se to interfere with the robustness of the epigenetic 32 regulation (ER) of differentiated phenotypes might suffice to generate malignant 33 transformation [3]. 34 Fully committed somatic cells can spontaneously reprogram to pluripotent stem-like 35 cells during the normal response to injury or damage in vivo [4]. Such cellular processes 36 involving dedifferentiation and cell-fate switching might constitute a fundamental 37 element of a tissue's capacity to self-repair and rejuvenate [5,6]. However, such 38 physiological (normal) cell reprogramming might have pathological consequences if the 39 acquisition of epigenetic and phenotypic plasticity is not transient. In response to 40 chronically permissive tissue environments for in vivo reprogramming, the occurrence of 41 unrestrained epigenetic plasticity might permanently lock cells into self-renewing By switching on or off different parts of the genome, ER is in fact responsible for the 49 variety of phenotypes in complex multicellular organisms (where all somatic cells are 50 genetically identical). Recent advances in experimental determination of ER 51 mechanisms have triggered an interest in developing mathematical models regarding 52 both ER of gene expression [13][14][15][16][17][18] and epigenetic memory [14][15][16][19][20][21][22]. 53

Fig 1.
Physiological and pathological cell fate reprogramming: A mathematical approach. Reprogramming-like phenomena in response to damage signalling may constitute a reparative route through which human tissues respond to injury, stress, and disease via induction of a transient acquisition of epigenetic plasticity and phenotype malleability. However, tissue regeneration/rejuvenation should involve not only the transient epigenetic reprogramming of differentiated cells, but also the committed re-acquisition of the original or alternative committed cell fate. Chronic or unrestrained epigenetic plasticity would drive ageing/cancer phenotypes by impairing the repair or the replacement of damaged cells; such uncontrolled phenomena of in vivo reprogramming might also generate cancer-like cellular states. Accordingly, we now know that chronic senescence-associated inflammatory signalling might lock cells in highly plastic epigenetic states disabled for reparative differentiation and prone to malignant transformation. We herein introduce a first-in-class stochastic, multiscale reduction method of combined epigenetic regulation (ER)-gene regulatory network (GRN) to mathematically model and computationally simulate how ER heterogeneity regulates the entry-exit mechanisms and kinetics of physiological and pathological cell fate reprogramming. (SAIS: Senescence-associated inflammatory signalling).
We are rapidly amassing evidence that, beyond the role of genetic alterations, 54 non-genetic stimuli such as inflammation or ageing, among others, can promote 55 epigenetic plasticity, namely, overly restrictive epigenetic states -capable of preventing 56 the induction of tumour suppression programmes or blocking normal differentiation -or 57 overly plastic epigenetic states -capable of stochastically activating oncogenic 58 programmes and non-physiological cell fate transitions including those leading to the 59 acquisition of stem cell-like states [18,23]. Disruption of normal epigenetic resilience has 60 the potential to give rise to each classic cancer hallmark [23]. An ideal therapy should 61 not only need to address chronic epigenetic plasticity of senescence-damaged tissues, but 62 also to "unlock" stem cell-like states to drive tissue regeneration. Unfortunately, most 63 current models and drug discovery strategies are largely biased on the prevailing, 64 mutation theory of the origin of cancer. Such a framework cannot capture the stochastic 65 aspects that drive the structure and dynamics of epigenetic plasticity and chromatin 66 structure that connects ageing and cancer, and might arise in the absence of initiating 67 genetic events. Schematic reprentation of the ER-GRN model and its multiscale reduction. (a): Gene regulatory network (GRN) of two self-activating, mutually-inhibitory genes with epigenetic regulation. In the GRN model, the gene product (denoted by X i in S2 Table) is its own transcription factor which, upon dimerisation, binds the promoter region of the gene thus triggering gene transcription. The transition rates corresponding to this GRN are given in S2 Table. For simplicity, we use an effective model in which the formation of the dimer and binding to the promoter region is taken into account in a single reaction, and the resulting number of promoter sites bound by two transcription factors is denoted by X ij (see S2 Table). Furthermore, depending on whether the epigenetic state is open (i.e. predominantly acetylated (A)) or closed (i.e. predominantly methylated (M)) the promoter region of the gene is accessible or inaccessible to the transcription factor, respectively. (b): Schematic representation of the time separation structure of the multiscale method developed to simulate the ER-GRN system. See text and S1 File for more details.
Identification of the molecular interactions controlling the transition between normal, 69 restrictive and permissive chromatin states is expected to have major impact in the 70 understanding and therapeutic management of the causative relationship between 71 ageing and cancer [9,24,25]. In this regard, our mathematical approach intends to 72 deconstruct and model the predictive power that epigenetic landscapes might have for 73 the susceptibility of cells to lose their normal identity. In order to determine the key 74 mechanisms underlying epigenetic plasticity and its connections with stem cell-locking, 75 we consider a gene network model that regulates the phenotypic switch between 76 differentiated and pluripotent states. Each gene within this regulatory system is acted 77 upon by epigenetic regulation which restricts/enables its expression capability (see Fig. 78 2(a) for a schematic representation). Crucially, the precise deconstruction of the highly 79 complex mechanisms through which epigenetic plasticity allows cells to stochastically 80 activate alternative regulatory programs and undergo pathological versus physiological 81 cell fate transitions requires the incorporation of a central role for epigenetic 82 heterogeneity in phenotypic plasticity [18,26] (see Fig. 1 for a schematic representation). 83 Pour et al. [26] have reported evidence that the potential to reprogram is higher within 84 select subpopulations of cells and that preexisting epigenetic heterogeneity can be tuned 85 to make cells more responsive to reprogramming factor induction. Along this line, we 86 recently presented a model of cell fate reprogramming, in which the heterogeneity of 87 epigenetic metabolites, which operates as a regulator of the kinetic parameters 88 promoting/preventing histone modification, stochastically drives phenotypic variability 89 capable of producing cell epistates primed for reprogramming [18]. However, how the  In the present work, we present a stochastic model of a coupled ER-GRN system 94 aimed at analysing the effects of epigenetic plasticity on cell-fate determination and 95 reprogramming driven by the heterogeneity of the ER system. Furthermore, we 96 introduce a stochastic model reduction analysis based on multiple scale asymptotics of 97 the combined ER-GRN system [27][28][29][30][31][32][33][34], which allows us to study a variety of different 98 behaviours due to heterogeneity in the epigenetic regulatory system [26] In this paper, we aim to study an ER-GRN model which can describe cell differentiation 121 and cell reprogramming. One of the simplest GRNs which allows to do this consists of 122 two genes, one promoting differentiation, and the other promoting pluripotency (see Fig. 123 2(a)). Nevertheless, in this section , we formulate our model and we analyse it 124 considering the most generic case, i.e. we assume to have an arbitrary number of genes 125 N G . By doing so, our theoretical analysis can be further applied to any ER-GRN model, 126 which implies a wide applicability of the derived formulation. However, when possible, 127 we try to relate the theory developed to our particular ER-GRN so as to keep track of 128 our case study.

147
In addition to TF regulation, we further consider that each gene is under epigenetic 148 regulation (ER). ER controls gene transcription by modulating access of TFs to the 149 promoter regions of the genes. In other words, in our model, ER is associated with an 150 upstream drive that regulates gene expression [44]. Such epigenetic control is often 151 related to alternative covalent modifications of histones. To address the high complexity 152 of ER, we focus on a simpler stochastic model of ER, based on that formulated 153 in [13,18] and [19]. Our model belongs to a wider class of models which consider that 154 single unmodified (U) chromatin loci can be modified so as to acquire positive (A) or 155 negative (M) marks. These positive and negative marks involve covalent modification of 156 histones. Of such modifications we consider methylation (associated with negative 157 marks) and acetylation (associated with positive marks) [19]. An illustrative example on 158 how epigenetic modifications, acetylation and methylation, alter the accessibility of TFs 159 to the promoter regions of the genes is shown in Fig. 2(a). Both modifications are 160 mediated by associated enzymes: histone methylases (HMs) and demethylases (HDMs), 161 and histone acetylases (HACs) and deacetylases (HDACs). For simplicity, we only 162 explicitly account for HDM and HDAC activity (see Fig. 2(a)). In our model, a positive 163 feedback mechanism is introduced whereby M marks help to both add more M marks 164 and remove A marks from neighbouring loci. The positive marks are assumed to be 165 under the effects of a similar positive reinforcement mechanism [16,19]. A full 166 description of the details of the ER model are given in Section S.2 of the S1 File (see 167 also [18]) and S3 Table, where the transition rates for the ER model are given.

168
Under suitable conditions, determined by the activity and abundance of 169 histone-modifying enzymes and co-factors, the positive reinforcement mechanism 170 produces robust bistable behaviour. In this bistable regime, the two possible ER stable 171 states are: an open epigenetic state where the levels of positive (negative) marks are 172 elevated (downregulated). In this case, the promoter of the gene is accessible to TFs and 173 transcription can occur. By contrast, in the absence (abundance) of positive (negative) 174 marks the gene is considered to be silenced, as TFs cannot reach the promoter. 175 An essential part of the stochastic dynamics of the ER system is the noise-induced 176 transitions between the open and silenced states. Escape from steady states is a 177 well-established phenomenon (see e.g. [45]) and thoroughly analysed within the theory 178 of rate processes [46] and large deviation theory [40,47,48]. As we will illustrate below, 179 these noise-induced dynamics are essential to classify the epiphenotypes of somatic 180 cells [18] and stem cells and unravel the mechanisms of reprogramming and locking.

181
Multi-scale analysis and model reduction 182 The system that results from coupling the ER and GRN models becomes rather 183 cumbersome and computationally intractable as the GRN grows. For this reason, in 184 order to analyse the behaviour of the resulting stochastic model, we take advantage of 185 intrinsic separation of time scales [27][28][29][30][31][32][33][34]. We exploit this time scale separation to 186 reduce our model by performing stochastic quasi-steady state approximations (QSSA) 187 by means of asymptotic analysis of the stochastic ER-GRN system as established 188 in [29][30][31]33] (see Fig. 2(b)). Especifically, we assume that the characteristic scale for  S1 Table for the definition of these variables).

192
Note that the assumption Y Z is exactly the Briggs-Haldane hypothesis for enzyme 193 kinetics [49] since the ER modification sites are the substrates for the ER enzymes (see 194 Section S.2 in the S1 File). The multiscale analysis is carried out in this Section, with 195 additional technical details provided in Sections S.3 and S.4 of the S1 File. Furthermore, 196 the corresponding numerical method is described in S1 Appendix. We show that, upon 197  Evolution to the silenced state of the pluripotency-regulating gene. Raw simulated data is generated by using the SSA on the model defined by the rates shown in S3 Table with parameter values given in Tables S.5 and S.6 in S1 File, for (a) and (b), respectively. Resulting fitted data correspond to the 100 ABC parameter sets that best fit the raw data.  The assumption that S Y allows for further simplification of the 208 model, as it allows to take the limit of S 1 in the stochastic equations for the TFs 209 monomers which leads to a piece-wise deterministic Markov description: the dynamics 210 of the number of TFs monomers is given by an ODE which is perturbed at discrete 211 times by a noise source [31]. 212 We present a summarised version of the asymptotic model reduction. Details of this 213 analysis are provided in Section S.3 of the S1 File. The starting point of our analysis is 214 the so-called Poisson representation of the stochastic process, which is equivalent to the 215 Master Equation, [30]: where X i denotes the product of gene i, X ij refers to the number of dimers of type j 217 bound to the promoter region of gene i and Y ij corresponds to the number of molecular 218 species of type j within the ER model of gene i. We also use the notation 219 X = (X 1 , . . . , X N G , X 11 , . . . , X N G N G ) and Y i = (Y i1 , . . . , Y i7 ). P(λ) ∼ Poisson(λ), i.e. 220 P(λ) is a random number sampled from a Poisson distribution with parameter λ [30], 221 R G and R E denote the total number of reactions in the GRN model (see S2 Table) and 222 in the ER model (S3 Table), respectively, with W ik and V ik denoting the transition 223 rates corresponding to the GRN model and the ER model (see S2 Table and S3 Table, 224 respectively). The stoichiometries r G ik , r G ijk and r E ijk denote the change in number of 225 molecules that reaction k has on X i , X ij and Y ij , respectively. Eqs. (1) and (2) are 226 associated with the the stochastic dynamics of the GRN (see S2 Table), which are 227 regulated by the ER part of the model. Eq. (3) describes the dynamics of the ER 228 system, which drives the dynamics of the GRN (see S2 Table and S3 Table).

229
Under the appropriate conditions, separation of time scales can be made explicit by 230 re-scaling the random variables and the transition rates. Based on our previous 231 work [18,32,34], we propose the following rescaling: In Eq. (4), the scale factors S, E, Y , and Z are the characteristic number of protein 233 transcripts, promoter region binding sites, histone modification sites, and epigenetic 234 enzymes (HDMs and HDACs), respectively. For simplicity, we assume that these scales 235 are the same for all the genes involved in the GRN. We assume that S E Y Z. 236 We further define a re-scaled (dimensionless) time: τ = b 11 ESt.  can assume that the (rescaled) rates associated with the fast variables GRN-ER 245 dynamics (Eqs. (6) and (8)) are much larger than those corresponding to their slow 246 counterparts (Eqs. (5) and (7)). Under these conditions, the stochastic dynamics of the 247 fast variables reaches their (quasi-)steady states while the slow variables are effectively 248 frozen [27][28][29]50]. 249 We proceed with the asymptotic model reduction by first addressing the QSSA

250
PDFs of the fast variables (see Inner solution below). We then move on to study the 251 QSS approximation of the slow variables, in particular, the large-S asymptotics of the 252 protein concentration dynamics (see Outer solution below).

253
Inner solution The inner solution corresponds to the relaxation dynamics of the fast 254 variables onto their quasi-equilibrium state, while the slow variables remain unchanged. 255 The solution of the inner dynamics allows us to determine the QSSA PDFs of the fast 256 variables conditioned to fixed values of the slow variables. 257 We proceed by considering the following rescaling of the time variable T = −1 1 τ .

258
Upon such rescaling, it is straightforward that all the rates of the reactions affecting the 259 slow variables (Eqs. (5) and (7)) are now O( 1 ), which implies that the slow variables, 260 x i and y ij (for j = 1, 2, 3), can be considered to remain frozen whilst the fast variables 261 reach their quasi-equilibrium distribution according to the dynamics: where l = 4, 5, 6, 7, 2 = O(1) and the slow variables, x i and y ij , j = 1, 2, 3, are 263 considered to stay constant.

269
Since the number of binding sites is a constant, the (quasi-)steady state distribution of 270 bound TFs to each promoter is a multinomial (see Section 3.2 of the S1 File for a 271 detailed derivation of this result). Otherwise, if gene i is epigenetically closed, then 272 X ij (T ) = 0 for all j with probability one. Therefore, the random vector describing the 273 number of TFs bound at the promoter region of gene i, B i , whose components are where N = (X 1 , . . . , X N G ) is a vector containing the monomer gene product of all genes, 276 multinomial PDF, whose generating function is given by: (12) where e i denotes the number of binding sites at the promoter region of gene i, and other 279 parameters are defined in S2 Table (transition rates) and S4 Table (rescaled parameters). 280 The quantity η i is defined as: Eq. (10) describes the inner dynamics of the fast components (enzymes and 283 enzyme-substrate complexes) of the ER system for each gene i. The resulting stochastic 284 dynamics describes how the enzymes switch between their free state and their complex 285 state at constant rates (since the number of the different substrates is constant under 286 the hypothesis of time scale separation). Since the number of enzymes is conserved, the 287 (quasi-)steady distribution of the number of enzymes of each type in complex form is a 288 binomial (see Section 3.4 of the S1 File for a detailed derivation of this result). The 289 corresponding generating functions are given by: where i = 1, . . . , N G and κ ij are defined in S4 Table. The number of free HDM and 291 HDAC molecules is then obtained from the conservation equations Outer solution The outer solution, corresponding to the dynamical evolution of the 294 slow variables, is obtained by sampling the fast variables, whose values are needed to 295 compute the reaction rates for the slow variables, from their QSSA PDFs (see Eqs.

296
(12)- (14)): The QSSA PDFs of the fast variables are conditioned by the current value of the slow 298 variables. We complete our asymptotic analysis by looking at the large S behaviour of 299 the slow GRN variables (see Eq. (5)). We resort to a law of large numbers enunciated 300 and proved by Kurtz which states that S −1 P(Su) → u when S 1 [29][30][31]51]. We can 301 apply this result straightforwardly to Eq. (15), which eventually leads to the asymptotic 302 reduction of the full ER-GRN system: The resulting dynamics consists on a hybrid system where the dynamics of the TF 304 monomers, x i (τ ), Eq. (17), is described in terms of a piece-wise deterministic Markov 305 process [52,53], i.e. by a system of ODEs perturbed at discrete times by two random 306 processes, one corresponding to stochastic ER (Eq. (18)) and the other to TF dimers 307 binding to the promoter regions. The latter are sampled from their QSSA PDFs, Eq. 308 (11). The stochastic dynamics of the slow ER variables, Eq. (18) is in turn coupled to 309 the random variation of the associated fast variables (ER enzymes, HDM and HDAC, and complexes). The number of complexes, Y i5 and Y i7 , are sampled from their QSSA 311 PDFs, Eqs. (13) and (14). The corresponding numerical method used to simulate such 312 system is described in detail in S1 Appendix.  (24). Vertical blue (horizontal green) hatching denotes regions where the pluripotency (differentiated) state is stable. Diagonal pink hatching denotes regions where the undecided state is stable. Regions of the phase diagram where different hatchings overlap correspond to regions of bistability or tristability. In the labels in the plot, P stands for pluripotency, D stands for differentiation and U for undecided. This phase diagram was obtained using the methodology formulated in [54]. Parameters values: ω 11 = ω 21 = 4.0. Other parameter values as per Table S.12 in Section S.7 of the S1 File.

Transitions between ER states: minimum action path approach 314
Noise-induced transitions are essential to understand ER dynamics and their effect on 315 cell-fate determination [23]. Throughout the bistable regime, sufficiently large 316 fluctuations in the stochastic ER system will induce switching between the open and 317 silenced states. The rate at which such transitions occur can be described using 318 reaction-rate theory [46] and large deviation theory [47], which show that the waiting 319 time between transitions is exponentially distributed. The average switching time, τ s , 320 increases exponentially with system size, which in this case is given by the scale of ER 321 substrates, Y [47,48,55,56]: where C is a constant and S is the action of the stochastic switch. Eq. (19) is derived 323 from considering the probability distribution of the so-called fluctuation paths, ϕ(τ ), 324 which connect the mean-field steady states in a time τ . According to large deviation 325 theory [47,48], we have P (ϕ(τ )) ∼ e −Y A F W (ϕ(τ )) , which implies that the probability of 326 observing paths different from the optimal, i.e. the path ϕ * that minimises the action, 327 is exponentially supressed as system size, Y , increases. This means that, for large enough system size, the behaviour of the system regarding large fluctuations is 329 characterised by the optimal path, which is such that: An explicit form of the functional A F W (ϕ(τ )) can be given if the dynamics is given by 331 the corresponding chemical Langevin equation [57]: where B jt denotes a Wiener process, and the mean-field drift, f j (y i ), is . . , 7, and the rescaled rates v ik (y i ) are defined in Multi-scale analysis and model 336 reduction. In this case, the action functional A F W (ϕ(τ )) is the Freidlin-Wentzel (FW) 337 functional: The norm T is the 339 diffusion tensor. Using Eq. (22), the optimal value of the action, S, can be found by 340 numerical minimisation, which provides both the optimal or minimum action path 341 (MAP) and the rate at which the ER system switch state driven by intrinsic noise.

342
Details regarding implementation of the action-optimisation algorithm are given in 343 Section S.5 of the S1 File. A complete description of τ s requires to estimate the 344 pre-factor C, which is not provided by the FW theory, but can be easily estimated using 345 stochastic simulation. constants c ij (see S1 Table and S3 Table). Such an ensemble is generated using 352 approximate Bayesian computation (ABC) [58,59], whereby we generate an ensemble of 353 parameter sets θ i = (c ij , i = 1, . . . , N G , j = 1, . . . , 16) compatible with simulated data 354 for the epigenetic regulation systems.

355
Our approach follows closely that of [18], to which we refer the readers for a detailed 356 presentation of the implementation. To summarise, we start by generating synthetic 357 (simulated) data (denoted as "raw data" in Fig. 3) regarding the ER system of a genetic 358 network epigenetically poised for differentiation, i.e. open differentiation-promoting 359 genes and silenced pluripotency-promoting genes (see the example shown in Fig. 3).

360
This simulated data will play the role of the experimental data, x 0 , to which we wish to 361 fit our model. The data set consists of 10 realisations and 25 time points per realisation 362 for each of the N G epigene regulatory systems. For each time point, t i , we consider two 363 summary statistics: the mean over realisations,x(t i ), and the associated standard simulated) data to a subensemble average consisting of the 100 sets that best fit the 371 data, for the differentiation-and pluripotency-promoting genes, respectively.

372
The above procedure provides us with an ensemble of parameter sets that are 373 compatible with our raw data, i.e. such that they fit the data within the prescribed 374 tolerances. The heterogeneity associated with the variability within this ensemble has a 375 clear biological interpretation. The rates c ij are associated with the activity of the 376 different enzymes that carry out the epigenetic-regulatory modifications (HDMs, 377 HDACs, as well as, histone methylases (HMs) and histone acetylases (HACs)), so that 378 variation in these parameters can be traced back to heterogeneity in the availability of 379 cofactors, many of them of metabolic origin such as NAD+, which are necessary for 380 these enzymes to perform their function [18].

381
The generated kinetic rate constants are dimensionless, i.e. they are relative to a 382 given rate scale [18]. Such a feature implies that there is an undetermined time scale in 383 our system. This additional degree of freedom can be used to fit our model of epigenetic 384 (de-)activation to particular data. Since the global time scales associated with different 385 ER regulation systems may differ among them, our model has the capability of 386 reproducing different systems characterised by different time scales as previously shown 387 by Bintu et al. [14]. Different colours and black lines show the three clusters resulting from a k-means analysis discussed in Sections Co-factor heterogeneity gives rise to both pluripotency-locked and differentiation-primed states and Analysis of ensemble heterogeneity. Dots in plot (b) represent PERSs within the ensemble defined in Section ER-systems ensemble generation and analysis.

389
In order to focus our discussion, we illustrate the application of our model 390 formulation on the case we want to study: a gene regulatory circuit with two genes, one 391 whose product promotes differentiation and another one whose protein induces 392 pluripotency. These two genes are further assumed to interact through mutual competitive inhibition (see Fig. 2(a)). Although such a system may appear to be too 394 simplistic to describe realistic situations, there is evidence that mutual inhibition 395 between two key transcription factors controls binary cell fate decisions in a number of 396 situations [60,61]. Our results are straightforward to generalise to more complex 397 situations. 398 We proceed to analyse how ER sculpts the epigenetic landscape over the substrate of 399 the phase space given by the model of the gene regulatory network. The latter provides 400 the system with a variety of cell fates, corresponding to the stable steady states of the 401 dynamical system underpinning the model of gene regulatory network [62]. The 402 transitions between such cellular states, both deterministic and stochastic, depend upon 403 the ability of the cell regulatory systems to elevate or lower the barriers between them. 404 Epigenetic regulation is one of such mechanisms. Here, we examine how ER is affected 405 by ensemble variability associated with variations in the availability of the necessary 406 co-factors on which histone modifying enzymes (HMEs) depend to carry out their 407 function. In particular, we will show that such variability is enough to produce a variety 408 of behaviours, in particular differentiation-primed and stem-locked states.

409
The GRN model exhibits a complex phase space, including an 410 undecided regulatory state 411 We start our analysis by studying the phase space of the dynamical system underlying 412 our model of gene regulation, schematically illustrated in Fig. 2(a). Using the 413 methodology described in detail in Section S.3 of the S1 File, we have derived the 414 (quasi-steady state approximation) equations for the optimal path theory of the 415 stochastic model of the mutually inhibitory two-gene system [13]. Such equations 416 describe the most likely relaxation trajectories towards their steady states [55,56], under 417 conditions of time scale separations described in detail in Section S.3. For the two GRN, 418 they read as October 18, 2018 14/34 where q 1 and q 2 are the variables (generalised coordinates) associated with the number 420 of molecules of proteins, X 1 and X 2 . The re-scaled variables, q i and q ij , and the 421 re-scaled parameters, ω ij , β ij , and δ ij , are defined in S4 Table (see also Section S.3 in   422 the S1 File).

423
The multiscale analysis carried out in Section S.3 shows that the parameters p, p ∞1 424 and p ∞2 are such that p ∞1 p = e1 E and p ∞2 p = e2 E , where e 1 and e 2 are the number of 425 sites in the promoter regions of our two genes exposed to and available for binding by 426 TFs, which implies that p ∞1 p and p ∞2 p can be directly related to ER: p ∞i p → 0,   This phenomenon, so-called epigenetic plasticity, has been recently proposed as a major 464 driver for disrupting cell-fate regulatory mechanisms in cancer and aging [23]. We 465 further focus on the role of heterogeneity within the ensemble described in Section 466 Materials and methods (see also [18]).

467
In order to characterise robustness of the different ER systems within the ensemble, 468 we have focused on the analysis of the average transition times between the open and 469 closed ER states. We define τ 1+ (τ 1− ) as the average transition time for a differentiation 470 ER system-DERS-to switch from closed to open (open to closed). Similarly, the 471 quantities τ 2+ and τ 2− are analogously defined for the pluripotency ER systems-PERSs. 472 The results are shown in Fig. 5(a) and (b), where we present scatter plots of the 473 average transition times within the ensemble of DERSs (Fig. 5(a)) and PERSs (  Table 1. Minimum action values, S, corresponding to the optimal escape paths shown in Fig. S.6 of the S1 File (see Section Transitions between ER states: minimum action path approach and Section Co-factor heterogeneity gives rise to both pluripotency-locked and differentiation-primed states for details). Parameter values are given in Section S.7 of the S1 File.
ER   Table 1). Whilst the action value for PERS1 is about twice the value of 500 PERS2 in Fig. S.6(b), there is an over 8-fold increase when comparing the action values 501 of DERS1 and DERS2 in Fig. S.6(a). Similarly, when comparing the action S for the 502 open to closed optimal paths, we observe that the variability associated with the DERSs 503 ( Fig. S.6(c)) is also larger than the one in PERSs (Fig. S.6(d)). This property partly 504 explains the difference between Fig. 5(a) and Fig. 5 the heterogeneity within the DERS ensemble ( Fig. 5(a)). Blue cluster DERSs exhibit 507 optimal closed-to-open paths with larger value of the optimal action than that found in 508 their red cluster counterparts (see Fig. S.6(a), where DERS1 belongs to the red cluster 509 and DERS2 to the blue cluster). This property has the consequence that the 510 closed-to-open waiting time, τ 1+ , is longer for blue cluster DERSs.

511
To quantify the effects of bistable ER on the landscape related to the gene regulatory 512 system (see Fig. 4), we proceed to estimate the probability, Q, that the combined 513 activity of each pair of DERS and PERS within our ensemble produces a global 514 epigenetic regulatory state compatible with differentiation. DERS-PERS pairs with 515 high values of Q are associated with differentiation-primed states. By contrast, those 516 DERS-PERS combinations with low Q are identified with pluripotency-locked states. 517 We proceed forward with this programme by recalling that escape times from a 518 stable attractor in a stochastic multi-stable system are exponentially distributed [47,48]. 519 This implies that the PDFs for the escape times for both DERSs and PERSs are fully 520 determined by the corresponding values of τ 1± and τ 2± . We also assume that, for a 521 given ER-GRN system, the DERS and the PERS evolve independently of each other. 522 We consider the PDF of the waiting time associated with a scenario of full 523 reprogramming of the epigenetic landscape, τ P . Such a scenario assumes that the where 532 . P 1 (τ p ) and P 2 (τ p ) are the probabilities related to each of the landscape reprogramming 535 routes. The probability that the ER landscape has undergone reprogramming from 536 pluripotency-locked into differentiation-primed state within the time interval (0, τ P ], Q, 537 is thus given by: where in our case, τ P has been taken as the mean time of τ 1+ , that is, the mean time 539 for the differentiation ER systems (DERSs) to switch from the closed to the open state, 540 which is needed for a cell to differentiate. Furthermore, τ 1+ exhibits a larger range of 541 variability than the time for the pluripotency ER systems to switch from its open to its 542 closed state, which is also a necessary condition for differentiation to happen. 543 We investigate the DERSs belonging to the different clusters of Fig. 5(a) regarding 544 their likelihood to produce pluripotency-locked epigenetic landscapes (results shown in 545 Fig. 6). The analysis shows that DERSs within the red cluster (Fig. 6(c)) correspond to 546 differentiation-primed epigenetic landscapes (Q = 1) for all the DERS-PERS pairs. By 547 contrast, the blue cluster ( Fig. 6(a)) and the green cluster ( Fig. 6(b)) contain DERSs 548 associated with both differentiation-primed (large Q) and pluripotency-locked (small Q) 549 epigenetic landscapes. As discussed in the next section, the latter are more abundant 550 within the blue cluster.

551
Analysis of ensemble heterogeneity 552 We now proceed to analyse the patterns observed in our ensemble of ER systems 553 regarding both the differences between the three clusters observed in the ensemble of 554 DERSs (Fig. 5(a)) and the distinctive features that characterise pluripotency-locked 555 DERS-PERS pairs. In order to do this, we follow the methodology put forward by [18], 556 whereby ensemble statistics (cumulative distribution functions (CDFs)) of the 557 parameters c ij (see S1 Table) corresponding to the DERSs/PERSs associated with the 558 subensemble of systems exhibiting a particular behaviour are analysed. By comparing 559 such CDFs to either the general population (i.e. whole ensemble) or to different 560 subensembles, we can detect statistically significant biases, which allows us to identify 561 key parameters (and their biases) associated with the behaviour displayed by the focal 562 subensemble.  S3 Table). We 570 proceed to look for which c 1j there are statistically significant differences, by comparing 571 the CDF of each cluster with that corresponding to the whole DERS ensemble, and also 572 the CDFs of the clusters among them (see Fig. 7). Each of these two-sample 573 comparison is carried out by means of the Kolmogorov-Smirnov (KS) test. Statistically 574 significant differences were found in the cases we comment below. The p-values are 575 reported in Section S.6 of the S1 File.

576
Red cluster versus blue cluster. As discussed in the previous section, the differences 577 between DERSs within the blue and red clusters are essential to ascertain the main 578 features that distinguish differentiation-primed and pluripotency-locked systems. The DERSs exhibit both long τ 1− and τ 1+ (see Fig. 5(a)).

599
Significant differences between differentiation-primed and 600 pluripotency-locked ER landscapes 601 The quantity Q allows us to classify each pair DERS-PERS drawn from our ensemble 602 regarding their degree of resilience to switch into a state prone to differentiation. If Q is 603 larger than a threshold value T , the corresponding DERS-PERS pair is categorised as 604 differentiation-primed. By contrast, when Q < T , the DERS-PERS pair is classified as 605 pluripotency-locked. 606 We first proceed to compare within the whole population (without discriminating 607 between clusters) those DERSs such that Q ≥ T (differentiation-primed ER landscapes) 608 against those with Q < T (pluripotency-locked ER landscapes). We take T = 0. If we now restrict our analysis to those DERSs within the blue cluster (see Fig. 8), 623 we observe that the parameters whose CDFs differ significantly when splitted into  The results of the previous sections suggest a number of strategies to unlock resilient 635 pluripotency states which hinder differentiation. One of our main conclusions is that 636 such states of resilient pluripotency are mostly vinculated to DERS-PERS combinations 637 such that the DERS belongs to either the blue or the green cluster. In view of this, a 638 possible strategy in order to encourage differentitation-primed ER landscape consists on 639 changing a selected combination of parameter values according to a rationale provided 640 by the analysis carried out in the previous two sections. Our results are shown in Fig. 9. 641 One possible strategy consists on first transforming a blue cluster DERS into a green 642 cluster one, and then completing the reprogramming of the DERS by transforming the 643 resulting set into a red cluster DERS. A candidate strategy involves first changing a the green cluster. The second step is then to change a parameter that exhibits 646 significant difference between the green and red cluster. Taking the results of the 647 previous section into consideration, we consider the reduction of c 111 (unrecruited 648 deacetylation) and the increase of c 13 (unrecruited demethylation). The result of this 649 reprogramming strategy is shown in Fig. 9(a), where we show that a blue cluster DERS 650 is first transformed into a green cluster one (green square in Fig. 9(a)), and then, finally, 651 into a red cluster DERS (red square in Fig. 9(a)). The initial blue cluster DERS has 652 been chosen as the set with the largest value of c 111 , which has been shown to be a 653 significant difference when comparing the blue cluster to the red one, and the blue 654 cluster to the green one, leading to the idea that this property is linked to the blue 655 cluster (idea which is reinforced because c 111 is not significant when comparing the red 656 and the green cluster).

657
The efficiency of such a strategy to unlock resilient pluripotency and to encourage 658 differentiation is shown in Fig. 9(b) where we present statistics of the differentiation consists on increasing the value of c 13 (unrecruited demethylation). Such a strategy is 668 not obvious, since c 13 is not one of the parameters whose empirical CDF has significant 669 differences when DERS in the red cluster are directly compared with those in the blue 670 cluster. However, since the CDF of c 13 is significantly different when both the blue 671 cluster and the red cluster are compared to the green cluster, it is conceivable that 672 increasing c 13 without further intervention could reprogram blue cluster DERSs. The 673 result of this reprogramming strategy is shown in Fig. 9(a) (red diamond). Simulation 674 results shown in Fig. 9(b) (two step reprogramming) confirm the viability of this Besides variability associated with cofactor heterogeneity, our model allows us to 680 address the issue of variability regarding HME activity. HME activity is affected by 681 both normal physiological processes, such as ageing, and pathologies such as cancer. For 682 example, impaired activity of HDM and HDAC has been observed in relation to cancer 683 and ageing. Here, we analyse the impact of HDM and HDAC loss of activity on the 684 dynamics of differentiation. In particular, we simulate differentiation in our ER-GRN 685 model to obtain statistics of the differentiation time to assess the effect of loss of HME 686 activity. The simulations shown in this section have all been carried out using the 687 hybrid multiscale simulation algorithm described in the S1 Appendix.

688
In order to clarify the effect of loss of HME activity on the ER model, we first 689 consider the phase diagram of its mean-field limit in different situations (see [18] for

701
This property suggests that a possible strategy to promote differentiation would be 702 to decrease HDM activity, as this would drive the PERS into its closed phase whilst 703 allowing the DERS to remain within its bistability region. In order to assess this and 704 other scenarios, we consider a base-line scenario where the number of HMEs is exactly 705 equal to average, i.e. e HDM = e HDAC = Z. We then compare different scenarios 706 regarding the abundance of HDM and HDAC to the base-line scenario.

707
Contrary to what could be expected, simulation results show that the strategy of 708 reducing HDM activity alone beyond the PERS closing boundary further hinders 709 differentiation. As can be seen in Fig. 10(c), a decrease in HDM activity actually leads 710 to longer differentiation time (see also [13]). Similarly, Figs. 10(a) and (b), which show 711 statistics of the differentiation time, reveal that a decrease in both HDM and HDAC 712 activity also leads to an increment in differentiation times, that is, this strategy fails to 713 decrease the differentiation time below the base-line scenario. In both cases, such 714 hindrance of differentiation is the product of the increase in the opening times (τ 1+ ) of 715 the DERS. This effect occurs because, as HDM and HDAC activity is reduced, the 716 DERS is driven towards its closed-bistability boundary. Close to such a region, the 717 DERS closed state becomes more stable and thus the corresponding τ 1+ increases. By 718 contrast, further reduction of HDAC activity moves the DERSs system closer to their 719 open-bistable boundary, resulting in a reduction of the differentiation time. However, 720 since the differentation times remain above those corresponding to the base line HDM 721 and HDAC activity scenario, we conclude that loss of both HDM and HDAC activity 722 contributes towards hindering differentiation.

724
In this paper we have presented a model of epigenetic plasticity which has helped us to 725 uncover some of the details and mechanisms underlying epigenetic regulation of 726 phenotypic robustness, in particular regarding the robustness of pluripotent states. We 727 have further uncovered how epigenetic heterogeneity regulates the decision mechanisms 728 and kinetics driving phenotypic robustness in a stem-lock model of pathological 729 pluripotency. Our deconstruction of epigenetic plasticity and phenotypic malleability 730 provides crucial insights into how pathological states of permanently acquired 731 pluripotency can be therapeutically unlocked by exploiting epigenetic heterogeneity. 732 We have added an ER layer to previous approaches in which cell phenotypes were 733 associated with the attractors of complex gene regulatory systems and their robustness, 734 with the resilience of such attractors tuned by the presence of intrinsic noise, 735 environmental fluctuations, and other disturbances [35][36][37][38][39][40][41][42][43]. Our approach is based on 736 two main pillars: namely, a framework for the generation of the ensemble of ER systems, 737 and a multiscale asymptotic analysis-based method for model reduction of the stochastic 738 ER-GRN model (see Section Multi-scale analysis and model reduction). The ensemble 739 generation method allows the definition of epi-phenotypes based on epigenetic-regulatory 740 modes compatible with a given state of the whole ER-GRN system [18]. We initially  Fig. 3 and [18]). With such an ensemble generated, 746 we then proceeded to evaluate its hidden, intrinsic heterogeneity in terms of the physical 747 properties of the ER-GRN systems. This approach is closely related to the notion of 748 neutral networks formulated to analyse systems with genotype-phenotype maps [63][64][65]. 749 By making a number of assumptions regarding separation of characteristic scales, we 750 re-scaled both the variables and the parameters of the ER-GRN system, which allowed 751 us to discriminate the underlying separation of time scales and consequently construct 752 an asymptotic expansion, leading to a stochastic QSSA of the system. This 753 approximation reduces a rather complex stochastic system (Eqs. (1)-(3)) to a hybrid, 754 piece-wise deterministic Markov system (see Section Multi-scale analysis and model 755 reduction). Furthermore, our model reduction procedure gives rise to an efficient and 756 scalable, hybrid numerical method to simulate the ER-GRN system (see S1 Appendix). 757 Although the model reduction was formulated for a GRN with an arbitrary number of 758 mutually inhibiting genes, such a procedure is applicable to broader situations.

759
When analysing the behaviour of the mean-field limit of the GRN, it became 760 apparent that, even in the simplest case considered involving a gene regulatory circuit of 761 only two genes, the system exhibited a complex space that included several multi-stable 762 phases. We observed a regime of tri-stability where the expected stem-lock (pluripotent) 763 and differentiated steady-states coexist with a third state, the indecision state, in which 764 the expression level of both genes is very low. From a developmental perspective, the disengage self-renewal pathways [23]. The opposite situation of so-called overly 779 permissive epigenetic states is accompanied by lowered epigenetic barriers that allow the 780 promiscuous sampling of alternative cell states [23]. Yamanaka originally appreciated 781 the link between epigenetic heterogeneity and plasticity when aiming to explain the 782 extremely low efficiency of somatic cell reprogramming at the population level [66]. We 783 now know that an epigenetic predisposition to reprogramming fates exists in somatic 784 cells and, therefore, the potential to acquire stem cell-like traits might in part reflect a 785 pre-existing heterogeneity in cell states [26]. Furthermore, by perturbing the epigenetic 786 state of somatic populations via inhibition of some epigenetic enzymes (e.g., the histone 787 methyltransferase Ezh2, which catalyses repressive H3K27 methylation [67]), such 788 heterogeneity can be harnessed to fine-tune the cellular response to 789 reprogramming-to-pluripotency factors. Indeed, our findings support a scenario in 790 which a sub-ensemble of ER systems with higher reprogramming potential pre-exists 791 within the ensemble of ER systems compatible with a terminally differentiated cell state, 792 and that such a sub-ensemble could be harnessed by targeting chromatin-modifying 793 enzymes such as HDMs and HDACs.

794
A careful evaluation of our ensemble of ER systems (i.e., combinations of DERSs 795 and PERSs) concerning stem-locked systems associated with repressive ER states (see 796 Sections Co-factor heterogeneity gives rise to both pluripotency-locked and 797 differentiation-primed states, Analysis of ensemble heterogeneity, and Ensemble-based 798 strategies for unlocking resilient pluripotency) has concluded that DERS heterogeneity 799 had a stronger influence on such pluripotency-locked systems when compared with that 800 of PERSs. Accordingly, we found that the ensemble of DERSs can be divided into three 801 different clusters, with each one exhibiting distinct properties regarding stem locking.

802
The so-called red cluster appeared to generate differentiation-permissive ER systems 803 irrespective of their PERSs counterparts. By contrast, the so-called blue and green 804 clusters contained DERSs yielding pluripotency-locked ER systems irrespective of their 805 PERSs companions. In light of these findings, we conducted a detailed comparative 806 analysis to uncover the underlying, statistically significant differences between DERSs 807 within the differentiation-permissive sub-ensemble and those associated with the 808 differentiation-repressive epigenetic states (see Section Co-factor heterogeneity gives rise 809 to both pluripotency-locked and differentiation-primed states). This approach allowed us 810 to detect which kinetic ER parameters were key to determine whether a DERS within a 811 given system produces either permissive or repressive ER differentiation systems (see 812 Section Analysis of ensemble heterogeneity). Remarkably, the elucidation of the identity 813 of such critical regulators (see Fig. 11) would allow the formulation of strategies aimed 814 to unlock differentiation-repressive epigenetic states by solely changing the values of 815 such parameters (i.e., epigenetic cofactors). The feasibility of such strategies was verified 816 by direct simulation of the ER-GRN system using our hybrid simulation method.  (14) and (17)- (18) 848 provide the reduced version of the original stochastic model which allows for a far more 849 efficient numerical implementation of a complex ER-GRN stochastic system.

850
The asymptotic reduction of the full stochastic model provides the basis for a hybrid 851 numerical method with enhanced performance with respect to the stochastic simulation 852 algorithm (as illustrated in Figure S.4 in the S1 File). The current hybrid method is 853 based on that formulated in [68]. The numerical method proceeds through iteration of a 854 basic algorithm composed of the following steps:  3. Consider the stochastic dynamics of the slow ER variables, Eqs. (7). These 861 stochastic equations must be solved by numerical simulation using Gillespie's SSA. 862 We first set the corresponding time step, ∆τ , using the SSA.   Rate of transcription of gene i = 1, . . . , N G k i2 Degradation rate of the protein of type i = 1, . . . , N G b ij Binding rate of the homodimers of protein of type j onto the promoter region of gene i u ij Unbinding rate of the homodimers of protein of type j from the promoter region of gene i c ij Kinetic rate of the jth reaction corresponding to the ER system of gene i (see S3 Table) X i Number of transcription factor monomers of type i = 1, . . . , N G X ij Number of sites of the promoter region of gene i bound to a dimer of proteins of type j Y i1 Number of unmodified nucleosomes (U-nucleosome) associated with the ER system of gene i Y i2 Number of methylated nucleosomes (M-nucleosome) associated with the ER system of gene i Y i3 Number of acetylated nucleosomes (A-nucleosome) associated with the ER system of gene i Y i4 Number of free HDM enzyme molecules associated with the ER system of gene i Y i5 Number of methylated nucleosome-HDM enzyme complexes associated with the ER system of gene i Y i6 Number of free HDAC enzyme molecules associated with the ER system of gene i Y i7 Number of acetylated nucleosome-HDAC enzyme complexes associated with the ER system of gene i  Table 3. Transition rates associated with the stochastic dynamics of gene regulatory circuit. Note that the rate of binding of homodimers to the promoter region of gene i is modulated by the level of acetylation of gene i, Y i3 : if Y i3 is above the threshold Y 0 , gene i is open, i.e. the promoter is accessible to homodimers and TFs. By contrast, if the gene's acetylation levels decay, gene i is silenced and the promoter is inaccessible to gene-transcription regulatory dimers.

Rescaled variables
Dimensionless parameters Fig 7. Empirical CDFs for the whole ensemble of DERS parameter sets (magenta lines). This ensemble has been generated according to the methodology explained in Section ER-systems ensemble generation and analysis (see also [18]). We also show the partial empirical CDFs corresponding to each of the clusters from Fig. 5(a) (red, green, and blue lines). We analyse a total of 90 DERS parameter sets. The red cluster includes 31 sets, the green cluster contains 13 sets, and the blue cluster has 46 sets. For reference, we also show the CDF for a uniform distribution (black line). Empirical CDFs for the DERS parameter sets within the blue cluster. This ensemble has been generated according to the methodology explained in Section ER-systems ensemble generation and analysis (see also [18]). The DERSs within the blue cluster have been divided into two subsets: those such that Q < T (SC-locked, blue lines) and those such that Q ≥ T (non-SC-locked, orange lines), with T = 0.7. For comparison, we plot the CDFs of the whole DERS ensemble (magenta lines), and, for guidance the CDF corresponding to a uniformly distributed random variable (black lines).
(b) (a) Fig 9. Plots showing the effect of the different reprogramming strategies of blue cluster DERSs, as evaluated in terms of the statistics of the differentiation time (τ D ). (a) Two step reprogramming is illustrated by the green square (first step), which finally becomes the red square (second step). One step reprogramming is depicted as the red diamond (see Section Ensemble-based strategies for unlocking resilient pluripotency for details). (b) Comparison of τ D for the original DERS and the ones resulting from the reprogramming strategies. We consider a base-line scenario where the number of HMEs is exactly equal to average, i.e. e HDM = e HDAC = Z. We then compare the simulation results obtained for different scenarios regarding the different strategies to the base-line scenario. Parameter values: Z = 5 and Y = 15. Other parameter values given in Table  S.11, Section S.7 of the S1 File.  Tables S.7, S.8, S.9 and S.10, Section S.7 of the S1 File. Strategies to unlock pluripotent stem-like states in ageing and cancer. Epigenetic regulation heterogeneity of differentiation genes (DERS), but not that of pluripotency genes (PERS), was predominantly in charge of the entry and exit decisions of the pluripotent stem-like states (blue). The application of the hybrid numerical method validated the likelihood of epigenetic heterogeneity-based strategies capable of unlocking and directing the transit from differentiation-refractory to differentiation-primed (red) epistates via kinetics changes in epigenetic factors. (Note: The epigenetic parameters regulating the entry into robust epi-states throughout the entire ER-GRN system revealed a regime of tri-stability in which pluripotent stem-like (blue) and differentiated (red) steady-states coexisted with a third indecisive (green) state). (R: Recruited; U: Unrecruited).