A mathematical framework for the emergence of winners and losers in cell competition

Cell competition is a process in multicellular organisms where cells interact with their neighbours to determine a “winner” or “loser” status. The loser cells are eliminated through programmed cell death, leaving only the winner cells to populate the tissue. Cell competition is context-dependent; the same cell type can win or lose depending on the cell type it is competing against. Hence, winner/loser status is an emergent property. A key question in cell competition is: how do cells acquire their winner/loser status? In this paper, we propose a mathematical framework for studying the emergence of winner/loser status based on a set of quantitative criteria that distinguishes competitive from non-competitive outcomes. We apply this framework in a cell-based modelling context, to both highlight the crucial role of active cell death in cell competition and identify the factors that drive cell competition. Highlights Emergent winner/loser status tests for competitive outcomes in cell-based models Differences in biomechanical properties alone are not sufficient for cell competition Winners have both higher tolerance and higher emission of death signals than losers

cell competition criteria, we tabulate these outcomes in a viability matrix ( Figure 2).
In this paper, we will assess the viability of a cell population based on its survival frequency, which is a statistic summarising cell population growth (or decline) in simulations of cell-based models. Later on, in Section 3.3, we introduce its analytical analogue, the survival probability. Suppose, for the sake of illustration, that we have two cell types, labelled A and B, and we want to determine whether they satisfy the cell competition criteria. As Figure 1 suggests, we need to run at least two homotypic simulations, one per cell type, and one heterotypic simulation in order to measure their viability in homotypic and heterotypic conditions. We compute the homotypic survival frequencies asλ for cell types A and B from their respective homotypic simulations. Similarly, we compute the heterotypic survival frequencies from a heterotypic simulation aŝ for cell types A and B, respectively. The simulations thus yield four survival frequencies: 64λ A ,λ B ,ξ A|B , andξ B|A . If a survival frequency is below one half, then the cell population

Homotypic viability
Heterotypic viability where µ is the friction coefficient. The force acting on vertex i is given by where ∇ i is the gradient of an energy function E with respect to the spatial coordinates 148 of vertex i. We use the energy function presented in [21], which describes three major 149 biomechanical properties: cell elasticity, cell contractility, and cell-cell adhesion; The first term represents cell elasticity, i.e. the cell's resistance against deformation.

151
The parameters K α and S 0 α are the elasticity constant and the target cell area of cell α, 152 respectively, while S α is the cell area of cell α. The second term models cell contractility, with Γ α and L α corresponding to the contractility constant and the cell perimeter of cell 154 α, respectively. The final term represents cell-cell adhesion, which is implemented as where β ∈ ∅ signifies that i, j is a boundary edge. Furthermore, we impose that each 181 cell division results in cells that are of the same type as the mother cell, i.e. a cell of 182 type A divides into two daughter cells of type A. 183 We implemented the mechanical model within Chaste, an open-source simulation 184 package for computational physiology and biology [24] that includes a range of cell-185 based models [25]. We refer the reader to the following GitHub repository for the code 186 of the mechanical model: https://github.com/ThomasPak/cell-competition. 187

188
After constructing the heterotypic mechanical model, we now determine whether it 189 can generate competitive outcomes. We first performed a systematic parameter grid 190 search varying the parameters of only one cell type, but we did not find any statis-  each parameter, as well as its default value, are given in Table 1.  (4)), respectively.
105 128 1 313 Table 2: Count of homotypic and heterotypic viability outcomes for the parameter sweep, summarised using the viability matrix ( Figure 2).

212
Out of 2 809 parameter sets, 23 resulted in simulation errors because the timestep 213 was too large. Since this only represented a tiny proportion of the parameter sweep, 214 we excluded these parameters from our analysis. We summarised the outcomes for the 215 remaining parameters using a viability matrix in Table 2.

216
The majority of parameter sets resulted in outcomes on the main diagonal, account- suggesting that cell competition generally depends on an active mechanism of cell death, 258 such as apoptosis [28], and that mechanical cell competition is no exception in this respect [15,16]. We note that these results do not exclude the possibility of mechanical 260 interactions playing a role in cell competition. They do strongly suggest, however, 261 that passive cell death alone is an insufficient mechanism for cell competition and 262 that mechanical interactions must be paired with an active mechanism of cell death 263 to produce robust competitive behaviour. The results of Section 2 suggest that cell competition requires an active and non-266 autonomous mechanism of cell death. This observation is also supported experimen-267 tally [7,8,28]. Therefore, the aim in this section is to develop a modelling framework 268 for cell competition implementing such an active and non-autonomous mechanism for 269 cell death. The core idea is that cells exchange "death signals" with their neighbours 270 and that these signals are accumulated by the cell into an abstract quantity called the 271 "death clock". When the death clock reaches a threshold value, apoptosis is triggered. 272 We do not yet attach the death signal to a concrete biological mechanism because there 273 are multiple competing hypotheses regarding the mode of intercellular communication 274 that underlies cell competition, and because the mode of communication may depend 275 on the specific type of cell competition under consideration. 276 We first discuss our biological assumptions and modelling choices in Section 3.1, 277 before introducing the death clock framework in Section 3.2. In Section 3.3, we define 278 the survival probability and derive its analytic expression for a given death signal.

279
Crucially, the survival probability enables us to analyse the death clock framework 280 from a theoretical perspective and make predictions on the viability of cell populations.

312
At division, we sample a stochastic G1 duration, denoted as t * , from the G1 dura- where C is subject to the constraints that (i) t * ∈ [0, ∞) and (ii) E(t * ) = t G1 , with 315 t G1 the autonomous G1 duration. If apoptosis is not triggered by the death clock, 316 the cell spends a duration t * in G1 phase and then transitions into G2 phase. After 317 spending a fixed duration, t G2 , in G2 phase, the cell divides and the process repeats for 318 each of the daughter cells. 319 We model the accumulation of death signals using an ordinary differential equation 320 (ODE) model in which the death clock, denoted by τ (t), evolves according to the where f (t) ≥ 0 is the death signal experienced by the cell. At birth, the death clock of a cell is initialised to zero, i.e. τ (t = 0) = 0. The apoptosis rule is then Cell is in G1 phase and τ (t) reaches T † ⇓ initiate apoptosis , where T † is the death threshold. We define the survival condition as We note that there are two potential sources of uncertainty in the death clock

Survival probability 331
In order to predict the viability of a cell population, we must determine the proba- Assuming that f (t) is a non-negative integrable function, we define such that the value of the death clock at time t * is F (t * ). This lets us write the survival 340 condition as F (t * ) < T † . We define the pseudoinverse function of F (t) as so that we can reformulate the survival condition as Substituting this into Equation (13), and denoting the cumulative distribution function for the distribution of t * as Ψ(t), we obtain As a special case, consider the constant death signal f (t) = c, where c > 0 is a positive 345 constant. We then have  In this paper, we implement the death clock mechanism in two particular cell-based  We provide the code for both models in the following GitHub repository: https: 359 //github.com/ThomasPak/cell-competition.
and summarise its contents in Table 3. The state of the system, denoted S(t), is then Symbol Description We evolve the death clock for each cell α as where f α (·) is the death signal function and x α (t) is the "input vector" representing the 367 extracellular environment. Since the cell population is well-mixed, this environment is 368 composed of every cell except itself, i.e. x α (t) = y 1 (t), . . . , y α−1 (t), y α+1 (t), . . . , y N (t) (t) .

369
In addition, we define two discrete operations: cell division and cell death. When a 370 cell's age reaches its total cell cycle duration, the division operation is triggered which

395
Concretely, the G2 death signal is defined as where g(t) is the proportion of neighbouring cells in G2 phase, i.e.
and c is a positive constant. 398 We first investigate the effect that the G2 death signal model has on homotypic pop-399 ulations in Section 4.1. This is done by deriving an expression for the homotypic sur-   given a fixed total cell cycle duration? In order to investigate this question, we denote 434 the total cell cycle duration as t G , and define β as the fraction of the cell cycle that is 435 spent, on average, in G1 phase, so that Even though cell cycle phases are stochastic in the G2 death signal model, we found 437 that the death signal is not only relatively stable, but also predictable. In particular, 438 we observe that the system is ergodic, in the sense that the average proportion of cells 439 in G2 phase relative to the population well approximates the average proportion of the 440 cell cycle spent in G2 phase. More precisely, we state that the system is ergodic if, on Furthermore, if the system is well-mixed, then we can approximate g(t) as Combining Equations (24) and (25), we have so that the death signal is Applying the methodology of Section 3.3, we use this result to derive the homotypic 446 survival probability, denoted λ, as For an exponential cell cycle model more specifically, this becomes .
which can be interpreted as a normalised death threshold. Hence, we write the homo-450 typic survival probability as a function of two dimensionless parameters: We validate this expression via simulation in Section S4 of the supplementary material.
This analysis therefore predicts that a population is viable for all η > ln(2)/4, and for 464 η ≤ ln(2)/4 it is viable for extreme values of β and nonviable otherwise (Figure 4(a)). 465 1 The astute reader may note the discrepancy between the definition of viability based on survival probability versus the definition based on survival frequency (Section 1.2): λ = 1/2 is considered nonviable, whereasλ = 1/2 is considered viable. This subtle distinction is rooted in the theory of birth-death Markov chains but bears no significance on our argument so we will not go into it further.  For each simulation k, we computed the homotypic survival frequency, denoted by 471λ k , using Equation (1). For every unique parameter set, we averaged the homotypic 472 survival frequency as where N sim is the number of simulations for the given parameter set. 474 We expect that nonviable populations tend to have a survival frequency below a 475 half, i.e.λ k < 1/2, and vice versa for viable populations.  The left-hand plot in Figure 4(b) shows that the observed border between nonviable 480 and viable regimes closely matches predictions for the well-mixed model. We see that for 481 small η values, the survival frequency is asymmetrical with respect to β, with higher 482 survival frequencies for β < 1/2 than β > 1/2. The reason for this discrepancy is 483 discussed in Section S4.5. In short, for low η values, the rate of apoptosis is so high that the limiting factor is the number of cells susceptible to apoptosis, rather than 485 the survival probability. For small β values, cells spend less time in G1 phase and are 486 therefore susceptible for a shorter amount of time.

487
The right-hand plot in Figure 4(b) also shows good agreement between theory and 488 simulations for the vertex-based model, although the border is less finely resolved than 489 in the well-mixed case. We also observe the same asymmetry for small η values as seen 490 in the well-mixed model.

505
After analysing these derived quantities, we are able to derive the heterotypic prolifer-506 ation regimes in Section 4.2.5, which we validate computationally using the well-mixed 507 and vertex-based models in Section 4.2.6. Finally, in Section 4.3, we pull together the 508 analyses from Section 4.1 and this section to characterise the competition regimes.

509
Similarly to Section 2, we create a heterotypic population in the G2 death signal  We demonstrated for homotypic populations that the proportion of the cell cycle spent in G1 phase, β, is an important nondimensional parameter in determining the survival probability. Hence, in analogy with Equation (23), we define β A and β B for heterotypic populations such that: Furthermore, we assume that the ergodic property holds for both cell types separately.
For cell type A, we have # A cells in G2 # A cells ≈ G2 duration of A cells cell cycle duration of A cells and an analogous expression can be derived for cell type B. We denote the number of 519 A-type and B-type cells with n A (t) and n B (t), respectively, so that we can write the 520 fraction of cells in G2 phase for the whole population as We substitute Equation (36) and its analogue for cell type B to obtain To simplify notation, we define the weighted average Assuming that the population is well-mixed, i.e. that Equation (25) For cell type A, the death signal is thus approximated as f where we use the symbol ξ A|B (t) to denote the instantaneous survival probability at In order to derive the instantaneous heterotypic survival probability for the exponential 536 cell cycle model in particular, we first define the dimensionless parameters in analogy with Equation (30). We can then derive that the instantaneous heterotypic survival probabilities are for cell types A and B, respectively. For brevity, we omit the word "instantaneous" 538 going forward and use the symbols 1 − β and ξ A|B instead of 1 − β (t) and ξ A|B (t), 539 except when we wish to emphasise their time dependence. Furthermore, in the rest of 540 the paper we will assume an exponential cell cycle model, unless stated otherwise.

541
Comparing the expressions for the heterotypic survival probability and the homo-

551
We define the heterotypic survival difference between cell types A and B as 552 The sign of the heterotypic survival difference tells us which cell type is at a proliferative  We define winners and losers here in a weak sense; if the population were to repro-558 duce indefinitely, the winner cells would come to dominate the heterotypic population.

559
It is not specified whether the loser population is eliminated. The classical definition 560 of winners and losers, however, is based on the stronger condition of loser elimination.

561
In Section 4.3, we will refine our terminology and differentiate winners and losers into 562 more precise categories, which include classical winners and losers. 563 We also note that this definition of winners and losers relies on the assumption that t G,A = t G,B , such that differences in survival probability alone determine relative 565 proliferative success. In the general case, however, differences in the total cell cycle 566 duration can also affect the dynamics of heterotypic populations. For instance, a cell 567 type with a lower survival probability may become more abundant than the competing 568 cell type by dividing more rapidly. However, for the sake of simplicity we do not consider 569 such cases in this paper, and instead characterise population dynamics solely in terms 570 of survival probabilities.

571
To obtain an expression for the sign of ∆ = A|B , we substitute Equations (44) and (45)   572 into Equation (46) and rearrange to give Since exp(·) is a monotonically increasing function, we have sgn(exp(x) − exp(y)) = sgn(x − y). Applying the sign function thus yields We validate this result using simulations of the well-mixed and vertex-based models in 587 Section S6 of the supplementary material.

594
In particular, we define the homotypic survival difference as for cell types A and B, respectively. The homotypic survival difference compares the 596 fitness of a cell type in a heterotypic environment to its fitness in a homotypic environ-597 ment.

598
The sign of the homotypic survival difference indicates whether a cell type is more 599 or less fit as a result of the heterotypic interaction, compared to homotypic conditions.

600
If ∆ = A|B > 0, then we say that cell type A is more fit when competing with cell type We expand 1 − β (t) to give .
The denominator of the right-hand side is strictly positive, so we only need to consider For cell type B, we derive an analogous expression: Comparing Equations (53) and (54), we derive the following identity: shows that the sign of the homotypic survival difference is determined by the difference    types of competitive interactions (Table 4). After accounting for the fact that cell type 651 labels are arbitrary, we can group these types into five distinct categories: potentially leading to a scenario where a previously nonviable loser cell type is "rescued" by the interaction with the winner cell type and becomes viable. In direct competition, 683 the winners become more fit and the losers less fit, potentially leading to a previously 684 viable loser cell type becoming nonviable as a result of the interaction, which is one of 685 the cell competition criteria. We therefore expect any competitive outcomes to be the 686 result of direct competition.

687
As discussed previously, we can partition cross sections of parameter space using the 688 coexistence curve and the neutral competition curve. In Figure 5, we plot these curves Hence, when considering the long-term behaviour of the population, we can substitute 1 − β A for 1 − β into the heterotypic survival probability for cell types A and B to obtain the asymptotic survival probabilities: Comparing Equation (59) so that we can write the asymptotic survival probabilities more succinctly as Conversely, if cell type B is the winner, i.e. ∆ = A|B < 0, we have where ξ ∞ A|B is defined analogously to Equation (61). 714 We can now use the asymptotic survival probability to characterise the viability of 715 competing cell types in a heterotypic population. Assuming that cell type A is the 716 winner, we distinguish between the following outcomes: We thus have three distinct proliferation regimes for ∆ = A|B > 0. Three analogous pro-727 liferation regimes exist for ∆ = A|B < 0, for a total of six proliferation regimes overall. We hypersurface, two winner viability hypersurfaces and two loser viability hypersurfaces. 741 We visualise the heterotypic proliferation regimes using cross sections for particular values of β B and η B in (β A , η A )-space. In these cross sections, the hypersurfaces become the following curves: A winner curve: B loser curve: A loser curve: The B winner viability hypersurface does not map onto a curve in (β A , η A )-space be- interacting with cell type B and are therefore "rescued" by the competitive interaction.

763
This is also present in Cross Section I, but it is more visible here. We note that this 764 region is contained within the indirect competition sector because only an indirect 765 competitive interaction can increase the fitness of loser cells. 766 We see three regimes above the coexistence curve. Below the A winner viability 767 curve, the winning A-type cells are nonviable, which renders both cell types nonviable.

768
Above this curve, the winner A-type cells are viable. In this subspace, the survival of  conditions. Further details are provided in Section S8 of the supplementary material.

785
To estimate the survival frequency for a particular parameter set, we averaged the 786 heterotypic survival frequencies across repeated simulations as The results for the well-mixed model are given in Figure 6 for cell type A that extend below their predicted limits in all cross sections.

798
In Figure 7(b), we see significant deviations from the predicted proliferation regimes.

799
When comparing the plots for cell type A with the results for the homotypic proliferation 800 regimes in Figure 4

812
The first condition of the cell competition criteria is that both cell types are homo-813 typically viable, i.e. λ A , λ B > 1/2. In order to satisfy λ A > 1/2, we only consider the 814 parameter space above the homotypic viability curve, as shown in Figure 8. To sat- The second condition is that only one cell type remains viable when the two cell types 820 compete. This implies a nonzero heterotypic survival difference, i.e. ∆ = A|B = 0, splitting 821 the homotypic viability regime into the coexistence regime and the competition regime The competition regime is further subdivided according to which cell type is the winner.

824
The G2 death signal model is symmetric with respect to swapping cell type labels, so 825 the choice of winner or loser is arbitrary. Therefore, for ease of notation, we henceforth label the winner cell type with W and the loser cell type with L, such that ∆ = W |L > 0 827 by construction.

828
As we saw in Section 4.2.5, the viability of the winner cell type is determined by its 829 homotypic viability, which is guaranteed by Equation (71). Therefore, we only need to 830 impose further that the loser cell type is heterotypically nonviable, i.e. ξ ∞ L|W ≤ 1/2. We 831 define the loser elimination regime as and the loser survival regime as The loser elimination regime satisfies the cell competition criteria and is non-empty for 834 the G2 death signal model. In addition, we have validated the predicted proliferation 835 regimes with computational simulations. We therefore conclude that the G2 death 836 signal model is capable of producing competitive outcomes. 837 We can further refine the competition regimes by considering, in addition, the type of 838 competitive interaction. Figure 8 shows that the neutral competition curve, defined by 839 ∆ = L|W = 0, runs through the loser survival regime. We define the neutral competition 840 regime as The neutral competition curve separates the loser survival regime into two subregimes 842 where ∆ = L|W < 0 and ∆ = L|W > 0, respectively. In the case of ∆ = L|W < 0, the fitness of 843 losers is reduced by the winners, but not enough to cause loser elimination. We define 844 this as the incomplete cell competition regime In addition, we can partition the loser elimination regime into the complete cell 846 competition regime and the critical cell competition regime which is the threshold regime between complete and incomplete cell competition. The 849 common feature of complete, critical, and incomplete cell competition is that the win-850 ners negatively impact the losers. We group these regimes under the cell competition Finally, on the other side of the neutral competition curve we have ∆ = L|W > 0, where 853 loser cells have a higher fitness than in homotypic conditions. We denote this as the 854 indirect competition regime We plot the competition regimes in Figure 8 for Cross Sections I and II, and summarise 856 them in Table 5.

857
The competition regimes let us discriminate between different types of winners and 858 losers. We define complete winners, critical winners, incomplete winners, neutral win- tal observations that spatial mixing is required for cell competition [42], and has been 892 replicated in other cell-based models of cell competition [13].

893
Our derivation of heterotypic proliferation regimes is based on the assumption that We define the tolerance and emission of cell type X, respectively denotedη X and 915 d X , as follows: We can formulate the homotypic viability condition, 1/2 < λ X , usingη X and d X by 917 substituting the homotypic survival probability (Equation (31)) and rearranging: The last inequality reads as the condition that cells must have a higher tolerance than 919 emission to be homotypically viable. The biological interpretation is that cells must 920 be capable of tolerating the death signal that they themselves emit in order to survive 921 as a group. The loser elimination condition, ξ ∞ L|W < 1/2, can also be expressed using tolerance and emission. Denoting the winner and loser cell types using the labels 923 W and L, respectively, we substitute the asymptotic survival probability of the loser 924 (Equation (61)) to obtain This means that winner cells must emit death signals at a rate that loser cells cannot 926 tolerate in order to eliminate the loser cell type from the tissue.

927
To satisfy the cell competition criteria, we require that both cell types are homo-928 typically viable, i.e. d W <η W and d L <η L , and that the loser is eliminated, i.e. 929η L < d W . Combining these expressions, we can summarise the conditions on the model 930 parameters such that the cell competition criteria are satisfied in a single statement:  Figure 9: Diagram of competition regimes using the transformed parametersη X and d X , defined in Equation (82). The green dot corresponds to the neutral coexistence point. The same conventions apply as in Figure 8. See Table 5 for the legend.

941
These relationships can be verified visually in Figure 9, which shows the competition 942 regimes in transformed parameter space. tions related to proliferation rates result in cell competition, and others do not [7]. In 956 this view, the inhibition of apoptosis can be regarded as a mutation that results in an 957 infinite tolerance, without affecting emission. Indeed, it has been shown in experiments that inhibiting apoptosis prevents cell competition [7,28]. are not viable on their own. The challenge in producing this outcome experimentally, 997 however, is that we would first need to identify an intrinsically nonviable mutant cell 998 type to assume the role of the loser. We therefore propose an alternative experiment 999 that could potentially simulate this behaviour with known cell types. The prevalent hypothesis is that cell competition is a mechanism for maintaining 1011 tissue health by eliminating unfit cells. However, what is meant by "fitness" in this context is not clear [44]. The classical definition of fitness is based on reproductive success and early experiments indeed linked reproductive fitness to cell competition, with winner cells having higher intrinsic proliferation rates than losers [45,46]. However,