Bivalent chromatin as a therapeutic target in cancer: An in silico predictive approach for combining epigenetic drugs

Tumour cell heterogeneity is a major barrier for efficient design of targeted anti-cancer therapies. A diverse distribution of phenotypically distinct tumour-cell subpopulations prior to drug treatment predisposes to non-uniform responses, leading to the elimination of sensitive cancer cells whilst leaving resistant subpopulations unharmed. Few strategies have been proposed for quantifying the variability associated to individual cancer-cell heterogeneity and minimizing its undesirable impact on clinical outcomes. Here, we report a computational approach that allows the rational design of combinatorial therapies involving epigenetic drugs against chromatin modifiers. We have formulated a stochastic model of a bivalent transcription factor that allows us to characterise three different qualitative behaviours, namely: bistable, high- and low-gene expression. Comparison between analytical results and experimental data determined that the so-called bistable and high-gene expression behaviours can be identified with undifferentiated and differentiated cell types, respectively. Since undifferentiated cells with an aberrant self-renewing potential might exhibit a cancer/metastasis-initiating phenotype, we analysed the efficiency of combining epigenetic drugs against the background of heterogeneity within the bistable sub-ensemble. Whereas single-targeted approaches mostly failed to circumvent the therapeutic problems represented by tumour heterogeneity, combinatorial strategies fared much better. Specifically, the more successful combinations were predicted to involve modulators of the histone H3K4 and H3K27 demethylases KDM5 and KDM6A/UTX. Those strategies involving the H3K4 and H3K27 methyltransferases MLL2 and EZH2, however, were predicted to be less effective. Our theoretical framework provides a coherent basis for the development of an in-silico platform capable of identifying the epigenetic drugs combinations best-suited to therapeutically manage non-uniform responses of heterogenous cancer cell populations. Author summary Heterogeneity in cancer cell populations is one of the main engines of resistance to targeted therapies, as it induces nonuniform responses within the population that clears the sensitive subpopulation, whilst leaving unaffected the non-responsive cells. Although this is a well-known fact, few successful approaches have been proposed aimed at both quantifying the variability associated to cell heterogeneity, and characterising strategies that circumvent its drug-resistance inducing effects. Here we present a computational approach that addresses these issues in the particular context of targeting epigenetic regulators (specifically, chromatin modifiers), which have been proposed as therapeutic targets in several types of cancer and also in ageing-related diseases. Our model predicts that the more successful combinations involve modulators of demethylase activity (specifically, KDM5/6 and UTX). By contrast, strategies involving EZH2 activity are predicted to be less effective. Our results support the use of our framework as a platform for in silico drug trials, as it accounts for non-homogeneous response of cell populations to drugs as well as identifying which subpopulations are more likely to respond to specific strategies.

Heterogeneity is a primary cause of de novo and acquired resistance to targeted 2 treatments in cancer therapeutics [1][2][3]. There is ample evidence of genetic diversity 3 among cancer cells both within and across tumours have been shown to exist [4][5][6], with 4 genomic instability being the main engine of genetic heterogeneity [7]. However, it is 5 now apparent that genetic variability alone does not account for the whole variation in 6 responses to targeted anti-cancer drugs [8,9]. Non-genetic heterogeneity has two broad 7 sources, epigenetic and stochastic, the latter being due to intrinsic factors such as noise 8 in gene expression [10][11][12][13] and asymmetric cell division [14,15], and it can produce 9 phenotype diversity even in genetically identical cells [8,9,16]. Not surprisingly, the 10 analysis of the consequences of non-genetic cell-to-cell variability in the cellular response 11 to drugs and its potential impact for the treatment of human diseases including cancer 12 has become a crucial issue to understanding drug resistance phenomena and developing 13 new targeted agents [1,3,17]. 14 Beyond its well-documented role in development [18][19][20], bivalent chromatin is an 15 emerging epigenetic trait of cancer that has been identified as a druggable therapeutic 16 target [19,21]. Genomic studies have found that there exists significant overlap between 17 bivalent domains in pluripotent stem cells and regions in cancer exhibiting abnormal 18 patterns of methylation, which pinpoints bivalency as a promising marker of 19 tumourigenesis. Furthermore, different bivalent chromatin regulators have been 20 identified as targets for cancer therapeutics [22][23][24]. An outstanding example of such 21 regulator is EZH2 (a histone methyl transferase which is the catalytically active 22 component of the PRC2 complex), which has found been to be mutated or overexpressed 23 in a number of cancers [23]. Specifically, somatic EZH2 mutations have a major role in 24 promoting lymphoid transformation in germinal-center leukemias [21]. Another family 25 of epigenetic regulators that has emerged as a potential target are the so-called lysine 26 demethylases (KDMs), whose expression has been found to be dysregulated in several 27 types of neoplasms [22,24]. Flavahan et al. [25] have recently proposed bivalent 28 chromatin regulators as key factors in epigenetic plasticity, i.e. the epigenetically-driven 29 alteration of the stability of cellular states (phenotypes). Variations in the activity of 30 certain epigenetic enzymes such as EZH2 can lead to repressive states, where the 31 epigenetic barriers are raised resulting in cells locked in a certain state, regardless of the 32 presence of signalling cues instructing them to exit such states. By contrast, such 33 October 11, 2020 2/26 alterations can also produce permissive states in which the epigenetic barriers are 34 lowered so that cellular states change by the mere presence of random noise. The latter 35 states often lead to pathological plasticity and cell reprogramming [25][26][27]. 36 The so-called Waddington or epigenetic landscape, a conceptual framework allowing 37 the integration of both genetic and non-genetic variability, was proposed over seven 38 decades ago [28]. In analogy to the potential energy landscape of physical systems, 39 Waddington put forward a representation of a cellular state in the form of an effective 40 potential landscape where the local minima correspond to cellular phenotypes. 41 Phenotypic switch transitions occur when cells transits to other minima by jumping 42 over the barriers separating adjacent basins. Such transitions are driven by intrinsic or 43 extrinsic fluctuations [29,30]. Noteworthy, the topography of such landscapes depends 44 on a complex set of biochemical interactions and, in particular, on the values of the 45 corresponding set of rate constants [17,[31][32][33][34]. These constants strongly depend on 46 protein structure (e.g. the accessibility of a binding domain), which is encoded within 47 the DNA [32]. Thus, an epigenetic landscape emerges from a particular genotype [17], 48 and therefore can be interpreted as a representation of the genotype-phenotype 49 map [35][36][37][38][39][40]. Beyond their dependence on DNA sequence, in the case of biochemical 50 reactions involving enzymatic catalysis, the associated rate constants depend also on the 51 concentration of metabolic cofactors [33,34], which may differ among cells thus adding a 52 further source of cell-to-cell heterogeneity. 53 In order to shed some light on the role of bivalent chromatin as a therapeutic target, 54 we propose a mathematical model of a simple gene regulatory circuit with a bivalent 55 promoter (see Figures 1(a) and (b)). Specifically, our model addresses mutual nonlinear 56 feedbacks between marked nucleosomes and transcription factors (TFs), which drive not 57 only the rates of both addition of new marks and TF synthesis [41][42][43][44][45], but also the 58 heterogeneity of cellular populations [17,33,34,46,47]. 59 The determination of epigenetic regulatory mechanisms has triggered an interest in 60 developing mathematical models regarding both epigenetic regulation (ER) of gene 61 expression [27, 33, 34, 41-43, 45, 50-53] and epigenetic memory [50][51][52][54][55][56][57][58]. 62 Specifically, Sneppen et al. [41] were the first to propose a computational model to 63 address the coupling between transcription factor and epigenetic regulation and how it 64 produces ultrasensitive behaviour, but they did not consider bivalent chromatin. Such a 65 coupling has been further studied in the context of bivalent chromatin [42,43,45]. 66 Although these works, similarly to the model presented here, describe bistability in gene 67 regulatory systems with TF self-regulation, it should be noted that they consider a 68 self-activating TF with high intrinsic cooperativity (n = 3) [42,43]. Moreover, 69 noise-induced phenomena in gene regulatory networks have been extensively 70 studied [12,[59][60][61][62][63]. In particular, Biancalani and Assaf [12] have studied noise-induced 71 bistability in a regulatory network (toggle switch) of monomeric, self-activating, 72 mutually-inhibiting genes, but they did not consider ER. 73 Our model (whose main features are summarised in Figure 1) proves that bivalent 74 chromatin regulation can induce complex behaviour even in the simplest circuits of gene 75 regulation, which would not exhibit such features in the absence of bivalent regulation. 76 We consider a model of a self-regulatory, monomeric TF with a bivalent promoter (see 77 Figure 1(a)), characterised by a number of both positive and negative feedbacks 78 identified in [42,43] (and schematically illustrated in Figure 1(b)). The resulting model, 79 in spite its apparent simplicity, exhibits considerable complexity regarding its 80 nonlinearities and multiplicity of time scales. In order to make analytical progress, we 81 exploit the latter and, by means of an asymptotic model reduction and singular 82 perturbation analysis, we reduce the model to a one variable stochastic system ( Figure 83 1(c)). By doing so, we have been able to determine, both numerically and analitycally, 84 that the system exhibits three different behaviours: monostability, bistability, and noise-induced bistability [64] (Figure 1(d)). Specifically, we show that the activity of the 86 epigenetic enzymes directly controls the effective cooperativity of the gene regulatory 87 circuit, which impinges on its ability to produce bistable behaviour (Figure 1(e)).

88
Our analysis of a bivalent TF is based on the explicit calculation of the Waddington 89 landscape for the system under consideration, which allows us to exhaustively classify 90 the set of possible landscape topographies. We account for the effects of heterogeneous 91 behaviour in the epigenetic-regulatory system [46,47] by means of an ensemble approach 92 previously developed [33,34], in which we take advantage of the direct connection 93 between landscape topography and the rate constants that characterise the underlying 94 biochemical networks from which it emerges. After validating our ensemble approach 95 against quantitative experimental data [20,21,65], we finally address the issue of what 96 epigenetic strategies are more effective regarding their ability to overcome heterogeneity. 97 We conclude that those involving tackling the activity of demethylases are more efficient 98 than those that targeting the activity of methyltransferases.  3-methyl, to be referred as "me3", epigenetic mark to two specific lysine residues in the 103 tail of the H3 histone (H3K4 and H3K27, see Figure 1(a)). Each of these enzymatic enzymatic catalysis (see [66] and S1 File for details specific to our model). We apply 106 this model to address how specific enzymes catalyse the addition of the me3 mark to 107 each residue (MLL2 to H3K4 and EZH2 to H3K27 [25]), and the removal of me3 from 108 Figure 1 (preceding page). Summary of the model formulation and model results. Panel (a) shows a schematic representation of the epigenetic-regulatory modifications considered in our model (i.e. trimethylation of H3K4 and H3K27), which affect promoter/enhancer regions. We focus on a transcription factor (TF) which is self-regulatory and monomeric. This panel also shows the different states of bivalent chromatin (open (closed ) if H3K4me3 (H3K27me3) predominates, poised if both marks are present) and how they affect TF expression. Plot (b) depicts a scheme showing the feedbacks between ER modifications and TF (full details given in Section S.1.1, S1 File). Plot (c) displays a schematic representation of the model reduction approximation, based on the presence of multiple time scales. Such a separation of time scales allows for a hierarchical elimination of the faster variables (full details given in Section S.1.2, S1 File). Plot (d) shows results of direct simulation of the process (carried out by means of the Gillespie algorithm [48,49]). These simulations illustrate the four type of behaviours afforded by the stochastic model: in descending order, inactive TF, bistability, noise-induced bistability, and active TF. Finally, Plot (e) illustrates one of our main results, namely, that the TF exhibits emergent, non-linear behaviour (e.g. bistability) which could not be observed in the absence of bivalent chromatin. Specifically, we show that such behaviour is associated to the tunability of the effective cooperativity, n ef f , of the TF regulatory circuit: by altering the activity of the chromatin-modifying enzymes we can control the quantity n ef f . In the particular case shown in the figure, by varying K 11 (which corresponds to the K M of the KDM that removes H3K27me3), the effective cooperativity changes from n ef f > 1 (upper plot) to n ef f ≤ 1 (lower plot). The response curve, R(x), and the n ef f of the TF are defined by x n ef f λ ef f +x n ef f .  [41,54]. Modified residues enhance recruitment of epigenetic enzymes, providing 112 a positive feedback mechanism where residue modification increases the rate of further 113 me3 addition (see Figure 1(b)). Following previous models of ER [33,34], we account for 114 such positive feedbacks by taking the rate constants associated to each of the 115 me3-addition MM reactions to change linearly with the number of corresponding 116 modified residues (see Table S.1 in S1 File). Furthermore, recruitment of epigenetic 117 enzymes to H3K4 is also enhanced by the presence of transcription factors [41][42][43] (see 118 Figure 1(b)). Such positive feedback is also accounted for by assuming the rate 119 constants of the MM reaction associated with H3K4→H3K4me3 to change linearly with 120 the number of TF molecules. By contrast, TF presence hinders recruitment of 121 me3-addition to H3K27 [42,43] (see Figure 1(b)). To take into account this negative 122 feedback, we have assumed that the rate corresponding to the back reaction of the MM 123 reaction associated to H3K27→H3K27me3 depends linearly on the number of TF 124 molecules (see Table S.1 in S1 File).

125
Regarding the gene regulatory system, we consider a simple self-regulatory gene (see 126 Figure 1(a)), i.e. a gene such that its protein stimulates its own expression, modelled by 127 a simple Hill model (with Hill constant α = 1, since we are considering a monomeric 128 TF) [69] with natural decay and basal production of the TF molecule. As shown in 129 Figure 1, the feedback between TF dynamics and ER is accounted for by assuming that 130 the rate of TF binding to the gene's promoter region is proportional to the number of 131 H3K4me3, whilst the rate TF unbinding is proportional to the number of H3K27me3. A 132 fully detailed account of the model formulation is provided in Section S.1 in S1 File.
133 Asymptotic analysis and model reduction. The model resulting from the 134 previous discussion is rather complicated. In order to make analytical progress, we 135 exploit the presence of separation of time scales. In particular, we resort to asymptotic 136 model reduction techniques and singular perturbation analysis (see Figure 1(c) for a 137 schematic representation). A fully detailed derivation is presented in Section S.1 in S1 138 File, which we summarise here. We proceed with our analysis in two steps. First, by 139 considering the usual assumptions regarding separation of time scales and the associated 140 quasi-steady state approximations (QSS approximation) normally involved in the 141 analysis of MM and Hill systems (see [66,69], and S1 File, Section S.3, for the stochastic 142 version of the QSS approximation), one can considerably simplify the system.

143
Specifically, under assumptions considered in detail in S1 File, we can reduce the 144 dynamics of the whole system to that of three variables, namely, the number of TF 145 molecules and the number of H3K4me3 and H3K27me3 residues. In the case of the 146 mean-field limit, it reduces to the following system of three differential equations (see 147 Section S.1 in S1 File for a detailed derivation): where Also, , and D I = Λ I λ I K 12 q I q I + K 11 .

150
Note that, stemming from the assumption that 2 = Ω S Ω T 1, where Ω S and Ω T are 151 characteristic scales for the number of H3K4 and H3K27, and TF molecules, 152 respectively, the resulting reduced system, Eqs. (1)-(3), still exhibits slow-fast features 153 (see Section S.1 in S1 File). The same holds true for the stochastic version of the model, 154 which is described by the corresponding Master Equation (ME): where the operators H x , H A , and H I are given by The step operators are defined in the usual manner: [70,71].

160
In order to proceed forward we propose that the PDF Φ(x, p A , p I , τ ) be expanded in powers of 2 [64,72]. In S1 File, we show that, at the lowest order, Φ 0 can be written , with φ A and φ I such that H A φ A = 0 and H I φ I = 0. The WKB solution to the latter equations is given by with J = A, I. However, as usual with singularly-perturbed problems such as Eq. (4), 161 the lowest order approximation does not provide any information regarding φ x (x, τ ) and 162 we need to resort to the next order [72,73]. At steady state, by considering the 163 appropriate solvability conditions, we can show that, to the lowest order, where the effective potential, U, is given by: A fully detailed derivation of the asymptotic calculation is given in S1 File, Section 166 S.1.2.

167
October 11, 2020 7/26 The ensemble approach: Modelling heterogeneity and parameter sensitivity analysis. Our ensemble approach serves two distinct but not unrelated purposes, namely, to perform a parameter sensitivity analysis and to model heterogeneity observed in epigenetic regulatory systems [46,47]. Regarding the former, given the complexity of our model, we need to consider carefully robustness to variations in parameter values. We consider such issue using an ensemble approach introduced in [33,34]. Briefly, we generate an ensemble of parameter sets. A parameter set is described by means of a vector, θ, defined as (see Table S.2 in S1 File): Each component of θ is independently generated by sampling from a uniform prior which parameters present statistically significant differences between sub-ensembles, we 184 can identify which of them are key to produce the behaviour associated to a particular 185 sub-ensemble [33,34]. We refer the reader to S1 File, Section S.1.4.2, for a full account 186 of the details of the classification procedure and further statistical analysis.

187
Experimental evidence of heterogeneity in ER systems comes from a number of 188 sources. It has been directly observed in single-cell CHIP-seq experiments [47]. Other 189 studies have provided evidence that the de novo reprogramming potential is higher 190 within selected subpopulations of cells and that such pre-existing epigenetic 191 heterogeneity can be tuned to make cells more responsive to reprogramming stimuli [46]. 192 Furthermore, Beguelin et al. [21] have reported that the efficiency of drugs that target 193 chromatin regulators, specifically EZH2 inhibitors, is very heterogeneous: their IC 50 Beguelin et al. [21].

205
Validation of the ensemble model. We first validate our approach against 206 qualitative and quantitative results [20,21,25,65,74]. Figures 3(a) H3K27me3-marked TFs. We also observe that TF expression is largely suppressed in 244 bivalent systems, their levels being only slightly larger than that of H3K27me3-marked 245 TFs. The pattern of TF expression within the High GE expression shown in Figure 3(o) 246 fits that of DCs: H3K4me3-marked systems exhibit higher TF expression than that of 247 H3K27me3-marked and bivalent systems, but the overall levels of expression are much 248 more uniform than in the case of the Bistable sub-ensemble. The pattern of TF

261
These data can be compared with our ensemble model, as shown in Figures 3(d) and (e). 262 Figure 3(d) shows the sub-ensemble average of the number of H3K27me3-marked 263 residues (see S1 File, Section S.1.5 for details) for the Bistable, High GE, and Low GE 264 sub-ensembles. The predicted levels of H3K27me3 in the Bistable sub-ensemble doubles 265 those in the High GE. This is consistent with the relative increase of H3K27me3 in the 266 Y641F and Y641N mutants relative to WT (although in absolute numbers, our model 267 prediction is slightly lower than those measured experimentally, figure 3b in [21]). The 268 Low GE sub-ensemble exhibits a 4-fold increase in H3K27me3 with respect to the High 269 GE sub-ensemble, which is out of the scale observed experimentally.  3)) for each of the four sub-ensembles (see S1 File, Section 277 S.1.4 for full analysis). By inspection, whilst there are obvious differences between the 278 Low GE PDFs and those corresponding to the High GE and the Bistable sub-ensembles, 279 the difference between the latter ones is much less obvious. In fact, quantitative analysis 280 (see S1 File) shows statistically significant differences, but with p-values very close to 281 the significance level. This observation is consistent with direct in vitro measurements 282 of the Michaelis-Menten constant, K M , and the catalytic constant, k cat , of the 283 wild-type EZH2 and Y641F and Y641N mutants by [65], who found that in all variants 284 of the enzyme, the K M s retained similar values within the error bars.

285
Taken together, these results validate two essential aspects of our model. First, the 286 ensemble variability appears to be consistent with the amount of heterogeneity observed 287 in the experiments carried out in [20,21,65]. Furthermore, these results allow us to 288 identify the systems within the Bistable sub-ensemble (and probably also the  Bivalency allows for bistability and tunable cooperativity. Our analysis 293 reveals that our model of a bivalent TF exhibits an array of complex behaviours, which 294 are inaccessible to the gene regulatory system in the absence of bivalent ER (see [41]), 295 and give rise to a variety of topographies of the epigenetic landscape, U ef f . Such behaviour is associated with effective cooperativity n ef f > 1 (see Figure 1(e)), in 300 contrast to the monomeric self-regulating gene with no bivalent ER (which has 301 n ef f = 1). By means of a sensitivity analysis, we determine that variation of key 302 parameters associated with the activity of ER enzymes, specifically, B 2 , B 8 , K 5 and 303 K 11 (see Eqs. (2)-(3)), provides us with ability to control n ef f thus allowing to suppress 304 bistable behaviour, since variations of n ef f can drive the system through a saddle-node 305 (s-n) bifurcation, thus removing pathological equilibria (see Figure 4). Such result is 306 consistent with current experimental evidence, whereby controlling the activity of the 307 ER enzymes, specifically of EZH2 [21,75], has the effect of reducing 308 pathologically-heightened epigenetic barriers and allow for physiological differentiation. 309 Specifically, Figure 4 shows the change of the steady state concentration of TF as B 2 310 and B 8 vary. The s-n bifurcation causes sharp transitions at both the lower and upper 311 yellow boundaries of Figure (4)(b), which enclose the region of bistability. The variety 312 of dynamical behaviours exhibited by our system is depicted by several examples shown 313 in the panels: (b1) monostability with TF at ON state (high level of expression); 314 (b2)-(b3) with bistability behaviour; and (b4) monostability with the TF at OFF state 315 (low level of expression).

316
In the case of bistable behaviour, it is interesting to analyse the variation of the where we present simulations of the system for different system sizes, which show that, 349 in agreement with our analytical approximation, noise-induced bistability ensues when 350 the system size is reduced below its predicted critical value. Furthermore, we have 351 ascertained that variations in ER enzyme activity induces noise-induced bistable 352 behaviour. Figure S.5 in S1 File shows that modification of parameters controlling ER 353 enzyme activity (specifically, B 2 , B 8 , and K 11 ) remove mean-field bistability but 354 preserve bimodality of the PDF (i.e. bistable behaviour is noise-induced). However, 355 noise-induced bistability remains a scarcely robust phenomenon within the bivalent TF 356 model: the percentage of noise-induced bistable θs within the generated ensemble is 357 ∼ 0.4% (see Table S.4 in S1 File) and the percentage of bistable θs that respond to 358 parameter change by becoming noise-induced bistable barely reaches a 10% of the size 359 of the bistable sub-ensemble (see Figure S.5 in S1 File).

360
Tackling heterogeneity. Heterogeneity in cellular populations is one of the main 361 obstacles regarding resistance to targeted therapies in cancer [1,3]. Since the effects of 362 drugs is not uniform on a heterogeneous cell population, those cells whose response is 363 poorer become a resistant subpopulation that compromises drug efficiency.

364
Heterogeneity has thus become a major factor in developing strategies to circumvent the 365 emergence of resistance [76]. ER enzymes such as EZH2 have been recently proposed as 366 targets to treat certain types of B cell leukemias [21,75]. Since recent single-cell studies 367 have revealed that heterogeneity also occurs in ER status [47], it is necessary to 368 evaluate its effect on therapies that target ER enzymes.

369
Following [33,34], where we put forward that heterogeneity in ER enzyme activity  Figure 6 382 (see also Section S.1.4 in S1 File). 383 We characterise heterogeneity within the (sub)ensemble using two statistics, namely, 384 the posterior marginal PDF for each of the parameters (see Figure S.3 in S1 File) and 385 the linear correlations between them. The former allows us to identify systematic, 386 statistical significant bias in the values of each of the parameters associated with an 387 specific behaviour (in this case, bistability) [33,34]. In Figure 6(a), we plot the p−value 388 corresponding to the Kolmogorov-Smirnov (KS) test comparing the marginal PDFs of 389 the bistable sub-ensemble with those of the Low GE, High GE and noise-induced 390 bimodal sub-ensembles. Provided that p is smaller than the significance level, the 391 smaller the p-value, the more statistically significant is the difference between the 392 marginal PDFs. Such information allows us to identify the subset of parameters which 393 are more likely to drive the behaviour of the system from one sub-ensemble (e.g. 394 bistable) into another (e.g. Low GE or High GE).

395
In order to perform our analysis, we choose the subset of parameters for which the 396 KS test gives smaller p-values for the comparison of the bistable sub-ensemble with the 397 Low GE and the High GE sub-ensembles: A 3 , B 2 , K 2 , B 8 , K 8 , K 11 , and K 12 (see 398 Figure 6(a)). The motivation for such a choice is that this set of parameters are the 399 ones that are more relevant to bistable behaviour [33,34]. We have analysed the 400 response of the bistable sub-ensemble by measuring the proportion of bistable θs that 401 become monostable when these parameters are modified (one at a time) by a factor ν 402 (see Figure 6(d)). According to such a metric, our results show that the parameters that 403 fare better concerning their ability to beat heterogeneity are the rescaled 404 Michaelis-Menten parameters for the KDMs (i.e. K 5 and K 6 for KDM2/KDM5 which 405 remove me3 from trimethylated H3K4 residues and K 11 and K 12 for KDM6 which 406 remove trimethyl from H3K27me3). By contrast, our analysis shows that targeting the 407 activity of MLL2 or EZH2 is a much poorer strategy regarding the aim of producing a 408 maximally homogeneous sub-ensemble response.

409
Single parameter variations have limitations of two kinds. First, when considering 410 ν > 1, the response curve is nonlinear but very quickly saturates to a response between 411 35%-65%, approximately (see Figures 6(c) and (d)). By contrast, if ν < 1, the response 412 curve is approximately linear. Whilst this scenario is more favourable than the ν > 1 413 one, ideally we would like to have a combination of both: a sigmoidal response curve 414 with saturation value close to 100%. Such scenario is achieved when considering 415 combinations of parameters (see Figure 6(e)). Specifically, of the combinations that we 416 have explored, the most efficient combination is K 6 -K 12 , i.e. those involving both 417 parameters regulating KDM activity. The reason for the enhanced efficiency of this 418 specific combination is that K 6 and K 12 are correlated (see Figure 6(b)). All other 419 combinations that we have tried involved uncorrelated parameters and perform worse 420 than the K 6 -K 12 one.  [1,3,76], few successful approaches [17] have been proposed aimed at 428 formulating strategies capable of quantifying the variability associated to individual 429 cancer cell heterogeneity and minimizing its undesirable impact on clinical outcomes.

430
Our current work accepts this challenge to provide a computational approach that 431 allows the rational design of combinatorial therapies involving epigenetic drugs against 432 cancer-driving chromatin modifiers [77][78][79].

433
To analyse the effects of targeted therapies which affect the enzymatic efficiency of 434 specific chromatin modifiers and combinations, thereupon we have formulated a 435 stochastic model of a bivalent transcription factor [42,43,45,58], for which we have been 436 able to derive the steady-state probability distribution function via singular bistable, and noise-induced bistability). In order to tackle heterogeneity, we 440 follow [33,34] and generate an ensemble of models which we then classify into four 441 sub-ensembles according to their qualitative behaviour. Comparison between analytical 442 results and experimental data allows us to determined that the so-called Bistable and 443 the High GE sub-ensembles show the same behaviour as undifferentiated and 444 differentiated cell types, respectively. By contrast, the Low GE sub-ensemble failed to 445 show any behaviour that could be identified as any cell type previously studied 446 in [20,21].

447
Since undifferentiated cells, i.e. locked in a self-renewing phenotype, have been 448 identified as aberrant phenotypes which, if left unchecked, will eventually give rise to a 449 tumour [21,25,33,34], we have focused on analysing the role of heterogeneity within the 450 Bistable sub-ensemble regarding their response to targeted epigenetic therapies. Such 451 therapies are assumed to affect the value of specific parameters which, in turn, alter the 452 enzymatic activity of the corresponding chromatin modifiers. Although our ensemble 453 approach provided a rationale to choose those parameters to which undifferentiated 454 (bistable) behaviour should be more sensitive [33,34], it was not surprising that 455 single-targeted strategies mostly failed to circumvent the therapeutic problems 456 represented by tumour heterogeneity whereas combinatorial strategies fared much 457 better. Specifically, those strategies involving the H3K4 and H3K27 methyltransferases 458 MLL2 and EZH2 were predicted to be notably less effective than those involving 459 modulators of the histone H3K4 and H3K27 demethylases KDM5 and KDM6A/UTX. 460 Accordingly, we are accumulating evidence that transcriptional differentiation programs 461 governed by KDM6A/UTX can enforce and safeguard cellular identity in several cancer 462 types whereas its pharmacological manipulation might impede tumour progression to 463 highly aggressive subtypes by blocking epigenetic roads to de-differentiation [80][81][82]. 464 Our theoretical framework provides a coherent basis for the development of an in 465 silico platform capable of identifying the epigenetic drugs combinations best-suited to 466 therapeutically manage non-uniform responses of heterogenous cancer cell populations. 467 The proposed paradigm might accelerate future tool developments on epi-drugs  . Panel (a) shows the equilibrium concentration of the TF (variable x) in the parameter space (B2, B8), the upper (resp., lower) sheet corresponding to the high (resp., low) level of expression, labelled as ON (resp., OFF) state. Panel (b) displays an enlarged view within the range contained in the cube in (a), namely 0 ≤ B 2 , B 8 ≤ 50. The region contained in the central orangish area corresponds to bistability. (b1-b4) Phase portraits showing the dynamics in the parameter space (x, q A ) for the combination of parameters B 2 , B 8 indicated with dashed lines in (b), computed using = 0.1. The orange and blue orbits correspond to the dynamics achieving the ON and OFF states, respectively. Note that both (b2) and (b3) fall into the region of bistability. Panels (c) and (d) show the relative size of the basin of attraction of the OFF equilibrium for parameter sets within the bistability region: (c) Parameter plane B 8 −B 2 , where K 11 is set to be 94.155549; (d) Parameter plane K 11 −B 2 , where B 8 is set to be 30. For the sake of visualisation, B 2 stands for a normalization of parameter B 2 , explained in Section S.2 in S1 File. Other constants are fixed according to the parameter set 2 (see Table S.7, Section S.2 in S1 File). Figure 5. Noise-induced transition to bistable behaviour.This figures shows stochastic simulation results supproting the analytical prediction that our system exhibits a noise-induced transition form mono-to bi-stable behaviour. The plot shows how, as Ω decreases (i.e. the level of noise increases), the PDF of the system undergoes a transition from mono-to bi-modality.