Different responses to cell crowding determine the clonal fitness of p53 and Notch inhibiting mutations in squamous epithelia

The growth and competition of cells in epithelial tissues plays an important role in both tissue homeostasis and the robustness of normal tissues to pre-cancer mutation. Whilst wild-type cells compete neutrally for dominance in the un-mutated tissue, naturally occurring mutations in individual cells may lend them a fitness advantage that can allow tissue colonisation. In mouse oesophageal epithelia, the growth of p53 mutants and a dominant negative mutant of the Notch downstream target Maml1 (DN_Maml1) have been shown to have different colonisation properties despite strong quantitative similarities in the growth of individual clones. Here we show that in order to recapitulate these behaviours whilst maintaining tissue turnover models need to take account of the response of cells to increased areal density in the tissue colonised by mutant cells. We demonstrate that p53 mutant clone growth approximates a logistic curve, but that without including limitations on mutation induced expansion the overall proliferation rate of the tissue drops due to space restrictions. In contrast, the ability of DN_Maml1 mutations to displace the wild-type population reflects a feedback that effects both mutant and wild-type cells equally. We go on to show how these distinct feedbacks are consistent with the distribution of mutations observed in human datasets.

Adult tissues rely on progenitor cells to maintain lifelong proper function of organs. A 2 major subject of study in biological research is the dynamics of progenitor cells in 3 squamous epithelia. These are rapidly regenerative tissues covering the external surface 4 of the body, the mouth, and the oesophagus and organised in layers of keratinocytes. 5 Progenitor cells communicate and interact with one another and their environment in a 6 tightly coordinated manner to respond to tissue needs and ensure homeostasis. 7 Progenitor cell dynamics has been shown to be accurately described by a simple 8 mathematical model, the Single Progenitor (SP) model [1][2][3][4]. This proposes that a 9 single, equipotent progenitor cell population maintains the tissue. In mouse stratified 10 squamous epithelia, proliferation is restricted in the basal layer where progenitor cells 11 stochastically divide or differentiate through stratifying into the upper layers of the 12 tissue before eventually being shed (Fig 1a). Division probabilities are balanced, 13 allowing homeostasis to be achieved across the progenitor population (Eq 1). 14 March 30, 2020 1/17 where A represents the basal layer progenitor cells, B the basal cells committed to 15 differentiate and C the suprabasal layer cells. Progenitor cells divide regularly with an 16 overall division rate λ and give rise to either two progenitor daughters (AA), two 17 differentiating daughters (BB) or one daughter of each type (AB) with fixed 18 probabilities. Given the fact that AA symmetric division leads to clone expansion and 19 BB symmetric division tends towards to clone extinction, the two symmetric division 20 rates should be equal in order for a steady state in terms of number of cells to be 21 maintained across the progenitor clone population. The probabilities of symmetric and 22 asymmetric divisions are r and 1 − 2r respectively with 0 < r < 0.5. Differentiating 23 daughters in the basal layer stratify to the suprabasal layer at rate Γ and supra basal 24 cells, C, are shed at rate µ. 25 In view of the fact that most common human cancers appear in epithelia [5], 26 understanding the mechanisms of cell growth and competition within the spatial 27 constraints of these tissues is fundamental to explain not only healthy tissue growth and 28 maintenance but also the processes of mutagenesis and cancer. Within the continuous 29 sheet of epithelium, cell populations constantly compete for space, minimizing 30 redundant or insufficient cell production and sustaining homeostasis. This may be 31 accomplished by mechanical interactions, e.g. through the mechanosensing ion channel 32 Piezo1 [6,7], and molecular signalling, e.g. through the Notch signalling pathway [8].

33
Whilst wild-type cells compete neutrally for space in the un-mutated tissue, naturally 34 occurring mutations in individual cells may lend them a fitness advantage that can 35 allow tissue takeover. Thus, non neutral competition between cells with different fitness 36 levels leads to the dominance of the fitter population ("winners") at the expense of the 37 less fit population ("losers"). This is a mechanism by which oncogenic mutations 38 colonize epithelial tissue and potentially drive preneoplasia and tumour formation. 39 Mutations in the tumour suppressor gene p53 and Notch pathway have been shown 40 to exhibit non neutral growth dynamics in epithelial tissues. Maml1 is required in the 41 canonical Notch pathway as it forms a complex with the Notch intracellular domain, 42 thus enabling the transcription of target genes [22]. An increasing body of work 43 suggests that p53 mutations are highly frequent in squamous epithelial cancers [9][10][11]. 44 Furthermore, recent studies report high incidence of p53 and Notch1 mutated genes in 45 normal human skin and oesophagus [12][13][14]. This highlights the importance of 46 studying the process of accumulation and interaction of these mutations within tissues 47 in order to understand the mutational landscape from which cancer develops. two-dimensional, hexagonal lattice was used to model the basal layer of the epithelium, 67 reflecting the observation that each oesophageal basal cell has six neighbours on average 68 (Fig 1d). Basal epithelial cells were simulated on a lattice originally containing 10,000 indicating a crowded region. A division event can lead to three potential outcomes: two 76 proliferating cells, two differentiating cells or one daughter of each cell type. The 77 neighbourhood in the SP CA model is defined by the six adjacent places. Division and 78 stratification events were considered as two independent processes determined solely by 79 r, λ and Γ parameters. This cell-autonomous approach could lead to cases where a cell 80 division event occurs at a region with no available neighbouring vacant space. In this 81 case, the two daughter cells were placed on the same grid space, indicating an increased 82 cell density area. Analogously, cases where an empty space generated by a recently 83 stratified B cell is not rapidly replaced by a nearby newly born cell might be observed, 84 representing a low cell density area. Considering the above, each lattice site could have 85 one of the following seven potential states: A, B, D AA , D AB , D BA , D BB and "empty" 86 (∅), where D AA , D AB , D BA , D BB correspond to double occupancies (Fig 1b,c). Thus, 87 the SP model (Eq 1) was extended to explicitly include space as follows: where A ∅ denotes a type A cell neighbouring a vacant lattice site and AX denotes 89 a type A cell neighbouring either a type A or type B cell, thus indicating that there is occupancies. The CA model was developed in NetLogo [15]. We used a Markovian 92 stochastic simulation algorithm where the basal layer was simulated as an asynchronous 93 CA. The algorithm included the following steps:  probabilities. Assign the division type as a next event for the selected cell. If a B 104 cell is selected, assign stratification as a next event for the selected cell. 105 5. If the next event is division, all neighbouring places are checked for empty space. 106 In the case of an existing neighbouring space, one new born cell will replace the 107 mother cell and the other will occupy the empty neighbouring space. If there is no 108 empty neighbouring space available then both will remain at the mother cell's 109 space (creating a "double state" cell), until a neighbouring space is released. If 110 stratification is the next event, B cell stratifies, leaving an empty space, which 111 allows potential neighbouring "double state" daughters to be released. The study of p53 and DN Maml1 mutant clone dynamics in mouse epithelia indicated 115 that the mutant progenitor clones are not in homeostasis in a mixed tissue, and have a 116 fitness advantage over their wild-type counterparts [9,16]. This has been shown to be 117 achieved by having a bias towards the production of proliferating progeny. Such bias 118 results in a gradual expansion in the proliferating population over time as there are less 119 chances for the mutant clones to be lost by differentiation. Thus, mutant clonal 120 expansion appears to be consistent with a SP model including a cell fate imbalance: where δ denotes the tilt in cell fate. Therefore, δ = 0 corresponds to homeostasis, δ 122 = 1 implies absence of symmetric differentiation, leading to persistence and δ = -1 123 implies absence of symmetric division, leading to extrusion.

124
To simulate mutant clonal dynamics in two-dimensional space, we used the spatial 125 SP model including a fate imbalance, considering both wild-type and mutant epithelial 126 cells. The choice of the grid architecture, the neighbourhood, the number of states and 127 spatial rules of the CA model were implemented as described in the homeostatic spatial 128 SP model. However, a new mutation-status property was introduced to distinguish 129 mutant cells to wild-type cells. In the case of p53 cells, mutants have the exact same 130 properties as the wild-type cells and the only thing that distinguishes the two cell 131 populations is that p53 mutant cells have an innate bias towards the production of 132 proliferating cells (δ) (Fig 1b,c). In the case of DN Maml1 mutant cells, mutants will 133 also have distinct λ, Γ and r values, as inferred by [9]. Considering the above, the 134 spatial SP model, initially described in Eq 2 was modified as follows, in order to 135 accommodate the mutant cell population: The simulation steps were modified as follows:    6. If the next event is division, all neighbouring spaces are checked for empty space. 151 In the case of an existing neighbouring space, one new daughter cell will replace 152 the mother cell and the other will occupy the empty neighbouring space. If there 153 is no empty neighbouring space available then both will remain at the mother 154 cell's space (creating a "double state" cell), until a neighbouring space is released. 155 If stratification is the next event, B cell stratifies, leaving an empty space, which 156 allows a pair of potential neighbouring "double state" daughters to be released. consisting of more than a defined number of cells would be considered as "crowded" 163 whereas fewer neighbouring cells than the defined crowding threshold would indicate an 164 underpopulated region ("empty"). In the former case, δ is negative rendering the fate 165 of dividing cells tilted towards symmetric differentiation to release crowding. In the 166 latter case, a positive δ is chosen, favouring symmetric division to fill the empty sites. 167 March 30, 2020 6/17 When the neighbourhood is neither "crowded" nor "empty", this indicates that the 168 local cell density is homeostatic and therefore δ = 0.
As each cell on the grid has six neighbours in normal conditions, i.e. no 170 overcrowding, no gaps, the local cell density crowding cut-off was set to six cells. This 171 reflects the topology of the healthy tissue. Furthermore, as wild type cells do not have a 172 δ, they would solely respond to the cell density bias (δ ), whereas in the case of mutant 173 cells the cell density bias (δ ) and their innate fate bias (δ) would be counterbalanced 174 (Eq 5). To avoid large increases in overall tissue cell density, an additional rule was 175 applied to every cell in the grid that had an overcrowded neighbourhood (n > 8 cells) . 176 In such case the probability of that cell to undergo symmetric division was minimized 177 (AA = 0, BB = 2r).

178
To perform simulations with the downstream feedback mechanism, no additional cell 179 density fate bias parameter (δ ) was used but mutants were set to lose their innate bias 180 towards the production of proliferating progeny, (δ), in response to crowding in their  The simulation outputs revealed that taking into consideration the spatial 212 constraints imposed by the lattice do not cause the model to deviate from the 213 experimental data taken from mouse oesophagus (Fig 2). The spatial model is able to 214 reproduce the characteristic features of the stochastic birth-death process (the hallmarks 215 of a single population of progenitor cells), as described formally in [3]. Homeostasis (a 216 constant number of proliferating cells) is maintained by a balance between cell 217 production and loss (Fig 2a). Over time, due to stochastic divisions resulting in a pair of 218 differentiated daughter cells, clones are lost following a simple relationship determined 219 by r and λ. To maintain a constant population, the average size of the persisting clones 220 rises linearly (with slope rλ/ρ). We observe that the simulated tissue demonstrates 221 both of these properties (Fig 2b,c) with growth dictated as expected by input model 222 parameters. Finally, whilst the clone size distribution becomes broader over time, the 223 shape remains constant once scaled by the average clone size (Fig 2d,e). Taken together, 224 these findings demonstrate that the average behaviour and the clone size distributions 225 of the epithelial basal cells are consistent with the SP model of epithelial homeostasis. 226 We further extended this testing to all available datasets, and find that the growth    further additions or fitting of the model, p53 clone growth and tissue takeover matched 246 experimental observations, showing an initial exponential growth slowing over time 247 (Fig 3a,b). Visualisation of the tissue showed increased packing of the mutant cells and 248 an increased abundance of progenitor cells, which over time lead effectively to the 249 growth becoming logistic. To confirm that the tissue size was acting as a limiting factor 250 and not the hexagonal grid, the model was re-implemented as a rule based model with a 251 carrying capacity. These simulations confirmed that whilst the presence of the grid 252 slowed growth slightly, the growth curves were ultimately logistic-like (Fig 3c). 253 Despite matching experimental data well, the simulations presented a new problem. 254 As the population of progenitor cells and double cells grows over time, the overall tissue 255 turnover slows as cells no longer have space to divide (Fig 3d). This is particularly 256 notable in the mutated regions which become packed earliest, and is not well supported 257 by biological evidence which suggests that the tissue is otherwise morphologically 258 normal. This has practical implications for the model; simulations can finish within the 259 lifetime of the mouse when the tissue becomes packed with double occupancies (Fig 3e), 260 whereas the tissue is found to have an increase in cell density of only ∼ 10% (Fig 3f).

261
To address these issues, we extended the model to explicitly enable feedbacks that 262 limit growth on the basis of crowding. Cell density increases in large areas replaced by 263 mutant cells, and is coincident with a return to neutrality in experiments [9]. Two 264 classes of feedback were tested based on the subpopulation of cells that respond to 265 crowding. In one class all cells respond to crowding by promoting differentation 266 (Fig 4a). In the second class only mutant cells respond to crowding by reducing their 267 propensity to symmetric division giving dividing cells (Fig 4b). If all cells respond to crowding, it implies that the genes responsible for coordinating this response either act 269 independently to or upstream of the mutated gene, that are active in the whole 270 population. Alternatively, if only mutated cells respond to crowding, it implies that the 271 genes responsible for coordinating the response act on targets of the mutated gene that 272 are not active in the wild-type subpopulation. We therefore term these mechanisms 273 "upstream" and "downstream" responses to crowding respectively to reflect the implied 274 relationship between the crowding responders and the mutated gene.

275
Both approaches were tested at a range of crowding thresholds and induction levels. 276 The upstream model enhanced tissue takeover, and no longer matched tissue takeover 277 by the end of the mouse lifetime. In contrast, the downstream model followed the 278 growth curves correctly and was broadly robust to the specific threshold selected 279 (Fig 4c,e). Furthermore, simulations no longer finished in the mouse lifetime and tissue 280 turnover remained roughly constant suggesting that the growth advantage of the 281 mutant p53 cells is sensitive to crowding, effectively acting as a negative feedback to the 282 mutant cell population.

283
Having shown that the p53 mutant cell growth can be described by a simple 284 "downstream" feedback, we revisited and retested DN Maml1 transgenic mutations in 285 mouse oesophagus. Testing with both feedbacks, we found that the "upstream" model 286 was required to capture the rapid takeover. This is consistent with prior observations 287 that showed through EdU labelling experiments that wild type cells neighbouring 288 mutant DN Maml1 clones were induced to differentiate and stratify [9]. Visual 289 examination of the underlying simulations suggested that the difference in growth 290 patterns between the upstream and downstream models arises as in the upstream model 291 the reduction in δ in all cells had the effect on ensuring that cells continued dividing, 292 but maintained the relative advantage of the mutant clone (Fig 4d,f). This effect, where 293 DN Maml1 mutation gives a constant fitness advantage whilst the p53 mutant 294 advantage is transient, was subsequently confirmed using a spatial Moran model. Here 295 each cell has a pre-defined fitness, and at each step of the simulation a cell has an 296 opportunity to replace its neighbour based on their relative fitnesses. We find that 297 whilst DN Maml1 mutations can be accurately modelled using the Moran model, p53 298 clone growth could not be accurately modelled with a single fitness (S2 Fig, S3 Fig). Having shown that both sets of mutations can be understood in terms of the 301 interactions of the underlying pathways and the tissue, we sought to assess how a single 302 model could enable both feedbacks behaved and how the mutations coexist and interact 303 in the tissue. In humans the prevalence of mutations varies between normal oesophagus 304 and the cancers that derive from it. p53 mutants are present the great majority of 305 cancers but are typically present in less than 10% of normal cells [12,13]. This implies 306 cancers develop from the p53 mutant population. In contrast NOTCH1 mutations are 307 several fold more common in normal oesophagus than in cancers, hinting that wild type 308 NOTCH1 may promote malignant transformation. To explore the interactions of the 309 mutations we designed a series of different mutational events to consider how the 310 different mutations interact; co-induction, timed inductions, and variable levels of 311 relative induction. As DN Maml1 represents a gene with control over several Notch 312 related genes, we modelled the effect of NOTCH1 mutations as an increase in the 313 imbalance combined with the measured increase in division rate, in order to examine 314 how the different growth behaviours in the tissue alter the competition. 315 We found that in all scenarios DN Maml1 mutant clones eventually took over the 316 tissue, excluding p53 clones from the tissue (Fig 5a,b). p53 clones that had grown for 317 extended periods before DN Maml1 mutations were introduced regressed more slowly, 318 as the larger clones were more slowly outcompeted but ultimately the constant fitness 319 advantage offered by DN Maml1 mutations eventually led to the loss of p53 clones.

320
However, a different p53 favouring competition scenario, where p53 mutants were 321 introduced at a substantially higher proportion, led to a higher rate of p53 clone 322 regression compared to scenarios of delayed DN Maml1 induction (Fig 5c).

323
In seeking to quantitatively describe potential distinct spatial properties of the two 324 competing mutant populations, we observed that the majority of p53 mutant cells were 325 consistently found in boundaries with DN Maml1 or wild-type populations, implying 326 that p53 mutant cells were not able to form coherent groups in space. On the contrary, 327 the proportion of boundary DN Maml1 mutant cells progressively decreased and was 328 almost minimized at later time points, as they colonized the entire grid (Fig 5d). p53 and DN Maml1 growth is explained by two distinct feedback mechanisms. The "upstream" (a) and "downstream" (b) feedback models simulating c: p53 mutant growth in back epidermis (SP parameters taken from [16]) and d: DN Maml1 growth in oesophagus (SP parameters taken from [9]). p53 mutant growth dynamics is appropriately recapitulated by the "downstream" feedback model whereas DN Maml1 by the "upstream" feedback model. e,f: Tissue takeover plots show the simulated and observed percentage of mutant cells. Data correspond to mean values across 100 simulations. Shaded areas correspond to SD. Simulation time lapses show p53 mutant growth with the "downstream" feedback (crowding threshold was set to 7 cells) and Maml1 growth with the "upstream" feedback (crowding threshold was set to 6 cells, δ = 1).
Furthermore, the proportion of p53 mutants' neighbours belonging to a different type 330 (cell mixing index) was consistently higher, indicating that p53 mutant cells were more 331 likely to share junctions with different cell types and therefore tended to be more 332 dispersed (Fig 5e). This is consistent with experimental observations [16].

333
The quantified distinct spatial features might be able to explain the behaviour 334 observed across different competition scenarios. When p53 mutants are introduced in 335 the grid in much higher numbers compared to their DN Maml1 competitors, they would 336 expand and form aggregates more easily and faster. However, as they cluster in space 337 and tend to appear in more crowded environments and encounter more p53 mutant 338 neighbours, they would start losing their innate proliferating advantage. The subsequent 339 contact with DN Maml1 mutant cells which start expanding substantially in the  Here, we explored the growth and competition of non neutral mutations in stratified 351 squamous epithelia using spatial models. Mutations in the tumour suppressor gene p53 352 and inhibition of the Notch pathway through DN Maml1 knock out were modelled. The 353 studied dynamics of p53 and DN Maml1 clones in mouse epithelia has been found to be 354 highly distinct. Whilst DN Maml1 knock outs are able to dominate the whole tissue 355 [9], p53 mutants initially expand but they eventually occupy no more than 30% of the 356 tissue [16]. We showed that in order to recapitulate these behaviours whilst 357 maintaining tissue turnover the spatial models need to take account of feedbacks 358 between neighbouring cells in the tissue. These are based on spatial regulatory rules 359 between surrounding cells in order to locally balance cell density. 360 We demonstrated that p53 clone growth approximates a logistic curve, but that 361 without including limitations on mutation induced expansion the overall proliferation 362 rate of the tissue drops due to space restrictions. In contrast, the ability of DN Maml1 363 mutations to take over the progenitor cell population reflects a feedback that affects 364 both mutant and wild-type cells equally. Thus, for reproducing p53 mutant behaviour, 365 we used a model where the feedback response was applied only to mutant cells, 366 indicating activation downstream to p53 mutations ("downstream" feedback). Such 367 model appears consistent with the p53 mutants' adaptable fate, suggested by their 368 observed incomplete tissue takeover. On the contrary, we modelled DN Maml1 369 dynamics using a feedback mechanism affecting all cell populations ("upstream" 370 feedback), consistent with the ability of DN Maml1 cells to evict their neighbours during 371 their early growth [9]. These spatial models of DN Maml1 and p53 mutant growth clustered but eventually appeared more dispersed and less cohesive, consistent with 384 experimental observations [16]. The consistent winning behaviour of Notch mutations 385 comes in contrast to a recent study, suggesting that the persistence either Notch or p53 386 mutant clones is determined by time [24].

387
Taken together, these findings suggest that the distinct tissue takeover outcomes 388 may be attributed to the differential response of p53 and DN Maml1 mutations to 389 crowding. Moreover, according to the competition spatial simulations, the sensitivity of 390 p53 mutant cells to increased cell density environments might be the cause of their loser 391 status. The crowding sensitivity as a hallmark of loser cells is also supported by 392 previous studies [7,17,25]. Furthermore, the increased dispersion of p53 mutants at 393 later time points in competition simulations, indicates that this spatial emergent 394 property might be a result of their mixing with DN Maml1 populations. This would be 395 consistent with previous studies in drosophila tissues, correlating the probability of loser 396 cell elimination with the surface of contact shared with winners [26]. Further 397 experimental work on p53 and DN Maml1 competition in mouse epithelia would shed 398 more light on how these two mutants interact in mammal tissues.

399
The two distinct feedback mechanisms that describe mutant dynamics may suggest 400 the distribution of mutations observed in human datasets. Recent studies identified high 401 incidence of many of the frequently found mutated genes in tumours such as p53 and 402 NOTCH1 in healthy human skin and oesophagus [12][13][14]. Strikingly, a higher 403 frequency of NOTCH1 mutations is reported for healthy human oesophagus compared 404 to oesophageal squamous cell carcinomas whilst the opposite applies to p53 mutations. 405 Considering the consistent winning behaviour of NOTCH1 over p53 mutants in our 406 competition simulations and given the paucity of NOTCH1 mutations in tumours, it is 407 tempting to speculate that the aggressive fitness of NOTCH1 may offer a 408 tumour-protective effect [27] under some circumstances. If however the mutations may 409 act combinatorially, p53 mutations that precede a NOTCH1 mutations may both