Competition between cell types under cell cycle regulation with apoptosis

Competition between different cell types plays a crucial role in bacterial ecology, developmental biology, and tumor growth. We study how cells compete in 2D by using both a mean ﬁeld theory and particle-based simulations. We employ a mechanical model that incorporates a stylized form of cell cycle regulation to control cell division events. This is extended to treat multiple cell types with different rates of programed cell death (apoptosis) and characteristic cell-cycle control pressures. Analytic predictions for the invasion speed and the coexistence line are in agreement with simulations. Synchronization in the cell division/apoptosis events can emerge, leading to oscillations in pressure and cell-cycle activity. We also study the invasion or elimination of small (pre)cancerous colonies. We show how a Laplace pressure at the colony interface, here controlled by differential cell adhesion, shifts the coexistence line; there are cell types that can invade when starting from a large colony but will be eliminated if the colony is small.

Cell competition occurs in a variety of contexts. These range from essentially collaborative, e.g., mutual organization of tissues during development or wound healing, to competitive, e.g., different bacterial species challenging each other for space within some environment. Cell competition not only maintains tissue homeostasis, but can be seen as the origin of tumor development when cancer cells outcompete healthy host cells [1][2][3][4][5][6]. Recent experimental and simulation studies have shed light on how stresses impact cell growth and proliferation [7,8] and how homeostatic pressure differences influence the competition outcome [9,10]. Previous studies [11][12][13][14] have not incorporated any sort of internal biochemical/genetic regulatory network, e.g., including the cell cycle and its associated signaling and apoptotic pathways, insofar as these are not already considered part of the cell cycle. In what follows, we consider such a generalized form of the cell cycle. While pressure-dependent apoptosis is a well known phenomenon, it is important to understand that there are very few mechanisms by which pressure can directly affect the cell, apart from at extreme pressure when the cell mechanically ruptures. Mechanical rupture is a direct physical response at the level of a cell, but such direct physical effects are not the focus of the present work. Under less extreme conditions pressure does not directly influence a cell's decisions to grow, divide, or die. Rather these decisions are made by a complex genetic regulatory machine that we refer to as the * m.s.turner@warwick.ac.uk (extended) cell cycle. The cell does not possess a physical pressure sensor, directly coupled to growth and death. Rather it indirectly senses mechanical and/or osmotic stresses at the (supra)molecular level, e.g., via an effect on membrane focal adhesion complexes or cortical cytoskeletal activity; the precise mechanism is unknown. The state of these molecules is then signaled to the cell cycle via biochemical (not physical) signaling pathways. This information is integrated over some response time (analogous to τ r in the model we use below). The cell then makes a biochemical decision, based on this signaling. Note that this is not at all the same picture as a direct pressure sensor; such a thing does not exist in a cell in any meaningful sense. It is our belief that a stylized version of the cell cycle regulatory mechanism must therefore be included if we are ever to have any prospect of better integrating future biological experiments (e.g., the effect of genetic manipulations) into a hybrid biochemical-physical model that deals with the stresses induced in the physical tissue. Indeed such a model also seems necessary to create a conceptual bridge between cell biologists, who consider the cell a biochemical "computer" and physicists, who build models based on the Newtonian response of generalized fluids.
Competition outcomes in vivo depend both on how pressure affects growth and division and on the rate of apoptosis. These properties vary between cell lines or following mutations [15][16][17][18]. In this study we ask how pressure and apoptosis together control the competition output. Our primary question is to ask how the role of pressure-dependent apoptosis is tensioned against division rate in determining competition outcomes. We conclude that for "sharp" apoptotic response, in which cells that have been stressed for a long time are much more likely to die, the characteristic pressure underlying the division rate is qualitatively the more important.
Tumor development begins with one single mutated cell but not all such cells will end in tumor formation; many are eliminated due to some competition with host cells. Only those that proliferate lead to tumors [19][20][21][22][23]. With this in mind we will also study how initial conditions impact the competition outcome.
Model. We model cell growth and division using a mechanochemical model that incorporates a stylized version of cell-cycle regulation, similar to that previously used to study the unrestricted growth of colonies made up of a single cell type [24]. We first outline how this model can be generalized to multiple cell types and then how it can be extended to include apoptosis, before using it to studying cell competition. Our stylized model can certainly be extended to include more biochemically complex responses. Nonetheless, we see it as about the simplest possible version of a model that combines a (stylized) cell cycle with a physical description. The actual parameter values that control it are likely very different in different systems, as discussed below. Our model is distinct from others in the literature as it directly controls the cell volume, and does not involve density explicitly (although they are related as ρ = 1/v). It also correctly captures the feature that cell decisions, like division and apoptosis, are controlled by the cell cycle, an internal biochemical/genetic regulatory network. Such cell decisions are not directly controlled by physical variables, like pressure or density, although these do affect the cell cycle [25,26]. The cell cycle is an orbit in a high dimensional biochemical/genetic space. Here we simplify by projecting this onto a stylized two dimensional control space parametrized by (polar) variables r, the cell cycle activity, and θ an irreversible angular phase, assumed to advance at a rate proportional to the activitẏ In our simulations ω is assigned, post division, to each daughter cell and is drawn from a normal distribution with mean ω = 2π/τ div and standard deviation τ div /10, with τ div being the nominal division time of an unstressed cell at r = 1. The activity is assumed to obeẏ involving a cell cycle signal integration time τ r and a characteristic reference pressure p r , above which the cells eventually become quiescent and stop dividing. The cell volume change between division events is then a combination of growth, proportional to activity r with a dimensionless growth rate c, and reversion to a minimal quiescent volume V q , with a characteristic recovery time τ v Cell division occurs at θ = 2nπ with n ∈ Z and is a process that replaces one cell with two daughters, each with half the volume 1×V → 2×V/2 [24]. In what follows we nondimensionalize all times (rates) and lengths in units τ div = 1 and V q = 1. Cells are assumed to undergo programed cell death at rate K (r), controlled by the cell cycle [27][28][29][30]. We assume a simple form for the apoptosis rate for each cell with k (0) the death rate at r = 0. We choose a sharp apoptotic response near r = 0 by taking γ = 8 in what follows so that quiescent cells die more frequently than those that are more actively dividing [28,29,[31][32][33][34].
Our particle-based simulations in 2D use Eqs. (1) and (4) to characterize each cell, carrying its own value of r, θ and V . These equations are parametrized differently for each cell type: cells of type i have reference pressure p ri and apoptotic rate constant k (0) i . In general they could also be assigned different growth rate c i , quiescent volume (area in 2D) V qi , and cell cycle times τ vi and τ ri , although in what follows we set these last four parameters to be the same for all cell types, for simplicity.
These parameters can vary widely between different organisms and cell types. In microbial populations, such as E. coli, cells are known to have a strong resistance to pressure, with division only stalling during extended, extreme pressure conditions. Such conditions can occur in long channels, reminiscent of so-called "mother machine" devices [35]. In the context of our model this means one would expect these systems to have a high reference pressure p r and a long division recovery time τ r . They also tend to become quiescent, rather than apoptotic, and so would have a very small value of k (0) . In contrast, MDCK cells are known to be rather sensitive to local pressure and/or density and quickly switch to quiescent and/or apoptotic response [36]. Animal epithelial cells also seem to be rather sensitive to small pressure changes as evidenced by the growth profile in spheroids [33]. Both cases correspond to a smaller value of p r , shorter τ r , and larger k (0) . In the context of the transition to cancer it is generally agreed that lack of control of cell cycle regulation combined with apoptotic misfunction enables the invasion of tumor colonies [37]. We can model this by either imposing a higher p r or lower apoptotic rate k (0) for cancer cells, compared with wild-type cells. When there is intratumor heterogeneity, several tumor colonies can coexist [38]; although for simplicity, we consider only the case of two cell types in what follows.
Cells are physically modeled as disks with area V , having a strong repulsive interaction with all cells and a weaker, shortranged attraction with cells of the same type, as described in [24]. These serve to maintain confluence while minimizing cell overlap. Cells moving on a substrate experience a friction force, opposing their velocity, with friction constant ζ = 0.01, non-dimensionalized in units of the smaller reference pressure. We neglect active cell motility for simplicity. Our choice of sliding friction and a lack of motility are natural assumptions for nonmotile E. coli. Epithelial tissues can have strong friction with the substrate due to cell adhesions and active motile forces exerted by the cells are important in determining their motion, particularly near colony edges. While we neglect motility in the present work we have previously shown that our model qualitatively reproduces several key properties of motile colonies such as their cell area distributions, their cell division times, and the crossover from exponential growth to a growth regime in which the tissue boundary moves at constant speed [24]. Because motile stresses are not a focus in this work we neglect them. In what follows we assume that the B cells have the smaller reference pressure, without loss of generality. The dynamics is then overdamped and is numerically determined by using force balance to solve for 033156-2 the cell velocities. Cell division events create two overlapping disks, oriented along an axis with random direction, with the same combined area as their mother. Results from a typical simulation are shown in Figs. 1(a)-1(d); see also Supplemental Material movies [39]. Here synchronization in cell division/apoptosis events leads to oscillations in the pressure, as well as both the number and activity of the proliferating (A) species, with a timescale of a few division times. This synchronization appears quite robust as it persists across a broad regime of parameters, see Figs. 1(e)-1(g) and S1 to S9. We find that the amplitude of the oscillations increases with the cell cycle times τ r and τ v . The frequency of the oscillations increases with the death rate of the cell type but is independent of the reference pressure. This synchronization is reminiscent of developmental processes [40][41][42][43][44], although signaling mechanisms are likely to be more complex in this context. Oscillations like these cannot be recovered in literature models that do not contain a faithful representation of the cell-cycle oscillator as there is no oscillator to entrain (phase lock) in such models. We see this as an important feature of our work. Figure 2 shows the results from a suite of simulations performed under different conditions. The line where A and B cells coexist indefinitely is seen to be in good agreement with the continuum analysis developed below. See also Figs. S1 to S9 and SM movies 1-6.
We also study the invasion or elimination of small (pre)cancerous colonies, see Figs. 2(d)-2(f). Here a Laplace pressure arises due to a surface tension at the colony interface, here controled by differential cell adhesion; we set the attraction between A and B cells to zero. This shifts the coexistence line. A noteworthy result is that there are cell types that can invade when starting from a large colony but will be eliminated if the colony is small, i.e., they lie between the corresponding coexistence lines. As the initial colony size is decreased it becomes more likely that stochastic colonies go extinct purely due to chance, leading to less precise agreement with the noiseless continuum theory. These results are broadly consistent with what is known from experimental studies [6]: small tumors often do not invade and there is an established role in stresses from the tumor "microenvironment".
Continuum theory. With a view to better understanding the simulation results we compare with a mean field theory. Here we construct v to be the division-adjusted mean volume of a cell. In the mean field this can be writteṅ The term −r ln 2 reflects halving in cell volume due to cell division with rate r. An unstressed cell with r = 1 has division-adjusted mean volume v = 1 + τ v (c − ln 2). In what follows we choose v = 5 and τ r = τ v = 1. The continuity equation involves the local cell velocity u and must be adjusted to account for apoptosis at rate K (r) Herev/v accounts for changes in the mean cell volume, r ln 2 is the cell growth rate, and K (r) is the death rate. Note that cell division is volume preserving and so does not appear in Eq. (6) explicitly. We assume cells move on a substrate with a friction coefficient ζ , determining the gradient of pressure p assumed to be of Darcy-like form This constitutive relationship is appropriate for quasi-2D systems in which cells reside on a frictional substrate. For tissues in 3D, if an external frame of reference is absent, a covariant description may be required. In this case the leading order stress gradients would depend on second spatial derivatives of the local velocity field [45]. However, it remains an open question whether a Darcy-like term still dominates in cell colonies in 3D, e.g., due to the role of a semi-permanent extra-cellular matrix. The question of whether frictional interactions come from some sort of contact friction/lubrication or from direct adhesion is immaterial, at least at the level of Newtonian theories (like ours) that involve only a friction coefficient (or viscosity) although it is possible that non-Newtonian viscoelastic effects may also play a role.
In the presence of the flow field u the laboratory-frame cell activity r and volume v are given by the advected forms of Eqs. (2) and (5). In what follows we will study two cell colonies in physical contact with u i , r i , v i , and p i describing each of the two cell types i ∈ {A, B}. In order to determine the system dynamics in general one must simultaneously solve for all these variables using two sets of equations of the form of Eqs. (2) and (5) to (7) with continuity of pressure and velocity as boundary conditions. Channel flow. Following the geometry of most of our simulations we consider cells confined in a periodic 2D channel, aligned along the x direction. The interface between the two cell types is parallel to the y direction. We seek solutions in which there is a steady-state front of dividing cells moving with constant speed s. We transform according to z = x − st so that the boundary between A (z < 0) and B (z > 0) cells lies at z = 0 and the equations can be cast in 1D by symmetry. We use a prime ( ) to denote ∂ ∂z and write the comoving, advected versions of Eqs. (2) and (5) to (7) as Solving these equations and matching the boundary values in a general setting is numerically tedious. However, we antici-pate that cell types will often be quite closely matched, e.g., mutant strains of the same bacteria or cancer/healthy tissue where tumor growth rates are usually extremely slow. In such cases the boundary will be moving slowly and spatial variations in the pressure and other variables will be small. This motivates a perturbative approach in which we will expand the variables about their values deep in the bulk (on either side of an interface). Denoting small variables with a δ-prefix we have s = δs, u i = δu i , v i =v i + δv i , p i =p i + δ p i , and r i =r i + δr i with thev,p, andr values representing the bulk values at z → −∞ (the A phase) or z → ∞ (the B phase). At zeroth order the cells are stationary with homeostatic activitȳ r i , volumev i , and pressurep i that are completely determined by the death rate k (0) i and reference pressure p ri according to At first order, from which one can show with a recovery length q −1 i defined by The boundary conditions at the interface z = 0 then determine the problem and hence the interface velocity We compare Eq. (23) with the interface speed obtained from simulations in Fig. 3 and find good agreement close to the coexistence line. The discrepancy between the theory and simulations on moving away from the coexistence line are due to a combination of factors, including (i) failure of the perturbative analysis, (ii) the emergence of interface roughness in the simulations, and (iii) the emergence of synchronization (global oscillations) not captured in the theory that maps discrete division and death events into a continuous process. The coexistence line in Fig. 2 3. The mean interface speed s, positive when A is invading, is nonzero off the coexistence line and is in best agreement with the mean field theory Eq. (23) (solid lines) close to coexistence. In panels (a)-(c) p rA /p rB is varied at constant is varied at constant p rA /p rB = 1. For these cuts on the phase plane the coexistence line lies at p rA /p rB = 1 and k (0) A /k (0) B = 1. Here the periodic simulation box size is 200×30.
can invade another unless its death rate k (0) is more than two orders of magnitude greater. The significance of this for biology is in how the reference pressure seems to naturally emerge as a strong control mechanism for cell invasion, even though the fact that cancer cells typically have severely impaired apoptotic response is commonly understood as the primary reason for their invasiveness.
In this work we develop a model for competing cell colonies that employs a mechanical model combined with a stylized form of cell cycle control, incorporating apoptosis.
We study the competition between two cell colonies with different characteristic cell cycle control pressures and apoptosis rates. The outcome depends on both apoptotic rate constants and the ratio of the two characteristic cell cycle pressures. Using both simulations and a mean field theory we identify the coexistence line and show how the characteristic pressure can naturally emerge as the stronger control parameter. We further predict the invasion speed, away from coexistence. We also show how interfacial tension between different cell types can suppress the invasion of small cell colonies (tumours). Finally, we observe that the birth/death events of invading cells can partially synchronize, leading to the emergence of global oscillations in cell number, pressure, and cell cycle activity. Such a feature is inaccessible at the qualitative level in models that do not encode an explicit oscillator, and is reminiscent of similar synchronisation in developmental biology. We hope that these predictions motivate experimental studies on timescales of many division events not typically the focus of current experiments.