An Emerging Allee Effect Is Critical for Tumor Initiation and Persistence

Tumor cells develop different strategies to cope with changing microenvironmental conditions. A prominent example is the adaptive phenotypic switching between cell migration and proliferation. While it has been shown that the migration-proliferation plasticity influences tumor spread, it remains unclear how this particular phenotypic plasticity affects overall tumor growth, in particular initiation and persistence. To address this problem, we formulate and study a mathematical model of spatio-temporal tumor dynamics which incorporates the microenvironmental influence through a local cell density dependence. Our analysis reveals that two dynamic regimes can be distinguished. If cell motility is allowed to increase with local cell density, any tumor cell population will persist in time, irrespective of its initial size. On the contrary, if cell motility is assumed to decrease with respect to local cell density, any tumor population below a certain size threshold will eventually extinguish, a fact usually termed as Allee effect in ecology. These results suggest that strategies aimed at modulating migration are worth to be explored as alternatives to those mainly focused at keeping tumor proliferation under control.

r ∈ L node on the lattice σ ∈ {m, r} cell phenotype in the model (moving m and resting r) K ≥ 2d capacity of each lattice node (2d velocity channels and K − 2d resting channels) n(r, k) ∈ {0, . . . , K} cell number at a node r and time k (r, k) = n(r, k)/K cell density at a node r and time k (0 ≤ (r, k) ≤ 1) ρ(r, k) = n(r, k) /K mean cell density at a node r and time k (0 ≤ ρ(r, k) ≤ 1) k ∈ N automaton time step r d probability that a cell dies r b probability that a resting cell undergoes mitosis r s probability that a moving cell changes to resting 1 − r s probability that a resting cell changes to moving κ ∈ R intensity of the phenotypic switch θ ∈ (0, 1) critical cell density value at which the probabilities to switch from resting to moving and vice versa are equal P(·) probability T b average cell cycle time  1 LGCA model description Our LGCA is defined on a 2-dimensional square lattice L ⊂ Z 2 . To every lattice node r ∈ L, velocity channels (r, c i ), i = 1, . . . , b, are associated. The parameter b defines the number of nearest neighbors on the lattice. Here, we choose b = 4 and c i ∈ {(1, 0), (0, 1), (−1, 0), (0, −1)}.
In addition, a fixed number of rest channels, (r, c i ), i = b + 1, . . . , K with c i = {(0, 0)} is introduced. The LGCA imposes an exclusion principle on channel occupation, i.e. at any time at most one cell is allowed in each channel at every lattice node. Thus K defines the maximal capacity of cells per node. We represent the healthy cells by empty channels and model explicitly two tumor cell phenotypes, denoted by σ ∈ {m, r}: moving (m) and resting (r) tumor cells. The occupation numbers at time k, η i (r, k), i = 1, . . . , K, are random Boolean variables that indicate the presence (η i = 1) or absence (η i = 0) of a tumor cell in channel (r, c i ). The local configuration of cells at a node r at time k is described by a vector η(r, k) = (η 1 (r, k), . . . , η K (r, k)) ∈ {0, 1} K . The total number of tumor cells at a node r and time k is defined by where n m (r, k) and n r (r, k) are the number of moving and resting cells at a node r at time k, respectively. The dynamics of our LGCA arise from the repeated application of three operators: propagation (P), reorientation (O) and cell reactions (R). The propagation and reorientation operators define cell movement, whereas the cell reaction operator changes the local number of cells on a node. The composition R • O • P of the three operators is applied independently at every node r of the lattice and at each time k to evaluate the configuration at time k + 1.
The cell reactions operator comprises cell death and cell proliferation where the latter is affected by the migration-proliferation dichotomy. In detail, we consider three stochastic processes: phenotypic switching (S), cell death (D) and proliferation(B). For the sake of simplicity, the interaction between cells is restricted to the node itself. Moreover, we consider interaction rules which are independent of the cells' distribution to channels and solely depend on the total number of cells n(r, k) within the node. During one automaton time step, the number of cells on a node r changes by subsequent application of a switch (S), death (D) and proliferation (B) step, n = (n m , n r ) The precise update rules are described below. Since we assume that individual cells decide to switch, die or divide independently from each other, the corresponding transition probability for the cell number per node follows a binomial distribution.
• Phenotypic switch (S): Cells can either rest or move. The phenotypic switch depends on the local cell density 1 = n/K. Phenotypes are changed with probabilities r s ( ) and 1 − r s ( ) that denote the probabilities of a moving cell to become resting and vice versa. Moving and resting populations can exchange cells if there is enough free space on the lattice. Two steps are independently performed: (1) Define M 1 = min(n m , (K − b) − n r ) which is the potential number of cells that can switch from the moving to the resting phenotype, taking into account the effect of local volume exclusion. The transition probability that there are j 1 successful events 'moving→resting' is modeled by (2) Define M 2 = min(n r , b − n m ) which is the potential number of cells that can switch from the resting to the moving phenotype, taking into account the effect of local volume exclusion. The transition probability that there are j 2 successful events 'resting→moving' is given by Finally, the difference between successful switches 'moving→resting' and 'resting→moving' is used to update the number of moving and resting cells so that (n m , n r ) → (n m − (j 1 − j 2 ), n r + (j 1 − j 2 )) = (n m , n r ). Note that the switching process does not change the total number of cells but the ratio between moving and resting cells.
• Cell death (D): Each tumor cell is assumed to die with probability r d , independently from the other cells. Thus, the probability that j 1 moving and j 2 resting cells die, starting with n cells is P D ((n m , n r ) → (n m − j 1 , n r − j 2 )) = n j r j d (1 − r d ) n −j if j 1 ≤ n m and j 2 ≤ n r , j = j 1 + j 2 , 0 otherwise. (4) • Cell proliferation (B): Birth of cells is dictated by the migration-proliferation dichotomy, which means that cells are allowed to proliferate only when they rest, i.e. when they are positioned on a rest channel. In addition, there must be empty rest channels to place the offspring into. Therefore, M = min(n r , (K − b) − n r ) defines the potential number of resting cells that are allowed to proliferate starting from n cells after the switch and death processes. The transition probability that there are j new cells after the proliferation step is then given by Altogether changes occurring in the cell number of a given node between two consecutive times k and k + 1 can be encoded in a transition matrix P with where the entries of the transition matrix of the switch process P S is determined by (2) and (3), and the transition matrices P D and P B of the death and proliferation processes, respectively, are defined by (4) and (5).
Subsequent to the cell reactions step, a reorientation step takes place where cells are randomly redistributed to the channels, providing a new local cell configuration at each node. Because of the migration-proliferation dichotomy, only the cells of migratory phenotype move. Therefore, moving cells are distributed to the velocity channels, while resting cells are distributed to the resting channels, resulting in a configuration η R•O (r).
The final update step is the application of a propagation operator. Cells in the velocity channels, i.e. moving cells, are transported one lattice unit in directions determined by their velocities. The movement of cells in the propagation step is purely deterministic. The spatiotemporal automaton dynamics is therefore completely described by the microdynamical equation where the change in the occupation numbers is −1 if a channel looses a cell, 0 if nothing is changed and 1 if a channel gains a cell, respectively.

Scaling of the LGCA
We have chosen a nondimensional scaling with unit lattice spacing and unit time scale k ∈ N in our LGCA simulation. Thus, cells residing at velocity channels move one lattice unit per unit time step. Our dimensionless simulations can easily be rescaled such that the temporal and spatial scales fit to specific applications. In our case, cell death and proliferation probabilities are related to the real-dimensional average life time T d of a cell and the average cell cycle time T b , respectively, by where τ is a sufficiently small real-dimensional time step length. This scaling is chosen because of the following arguments, detailed for the case of cell death. Let X 1 , X 2 , . . . be independent Bernoulli trials with possible outcomes "cell dies" (X i = 1) and "cell survives" (X i = 0). The corresponding probability distribution is given by P(X i = 1) = r d . Let N be the number of trials up to the first time a cell dies, i.e. N := min{n|X 1 = . . . X n−1 = 0, X n = 1}. Then N is geometrically distributed such that P( If τ is the real-dimensional time step length, the life time of a cell is given by τ N . Hence, the average cell life time in the model is T = E[τ N ] = τ /r d which must be related to the real-dimensional average life time T d . Similar arguments hold for the real macroscopic time taken for a cell to proliferate (cell cycle time) and for a cell to change its phenotype. Please note that if the LGCA parameters r s , r d and r b depend on time τ , the transition matrix of the cell reactions becomes time-dependent: Besides scaling the LGCA time intervals, also the lattice spacing can be adjusted. If the lattice spacing is , the mean-square displacement of a cell per time step τ in the LGCA is proportional to 2 /τ . Thus, by scaling x = r ∈ R, the relationship between the LGCA motility and measured movement of cells is characterized by the diffusion constant D = 2 /τ .

Equilibrium state under fast switching assumption
When deriving the mean-field equation, we assume that the switch dynamics is much faster than the cell number changes due to proliferation and death, that is r s , 1 − r s r b , r d . Further, we consider the low density regime where the carrying-capacity effects can be neglected when considering the switching. Then, by time scale separation, the fractions p σ := n σ /n, σ ∈ {r, m} of resting and moving cells, respectively, are in detailed balance with respect to the switch dynamics. This amounts to Hence p r = r s or, equivalently, n r = r s n.
4 Mean-field approximation of the LGCA dynamics The LGCA model is composed of a cell movement and a cell reaction process. In the following, we separate the two processes and derive mean-field approximations for each one. We show that the cell movement in the automaton model can be approximated by a degenerate diffusion term, while cell reaction in the automaton leads to a nonlinear reaction term in the resulting partial differential equation description.

Scaling limit of the cell reaction process
In the following, we are interested in the time evolution of the cell density and neglect any spatial dependence between nodes. We assume that the system is in equilibrium with respect to the switching (fast switching assumption, see section 3). We consider the Markov process (N k ) k≥0 which describes the number of cells at a node r ∈ L in the discrete phase space {0, . . . , K}. The transition matrix of this process, derived from (6) and (10), is given by and Here, the minimum function min( . This approximation underestimates the minimum function but provides an intuitive explanation as it describes the fraction of cells that are able to proliferate times the fraction of space that can be filled. Please note that, for simplicity, we use a notation where the dependence of the switch function r s from the node density is dropped.
The probability to be in state j after one time step is W (i, j)P k (i) and in matrix notation P k+1 = W P k .
The transition of a microscopic process to a macroscopic description requires a temporal scaling relation between the macroscopic and microscopic variables (see section 2 for details). We assume a small parameter τ > 0 that scales the time variable t = kτ such that where now the transition matrix W depends on the real-dimensional time step size τ . Rearranging terms gives with intensity matrix Q := lim τ ↓0 where Q D and Q B are the intensity matrices for the death and birth processes, respectively. A direct calculation gives and similarly From the master equation (13), the temporal evolution of the mean cell density ρ t = N t /K is obtained by Assuming that E[N 2 t ] = E[N t ] 2 , which holds if the variance of N t is small, for instance if K is large, we finally arrive at a macroscopic description of the LGCA birth-death process where R d = 1/T d and R b = 1/T b are the real death and proliferation rates, respectively, see section 2.

Scaling limit of the cell migration process
In the following, we consider the automaton cell migration process without cell reactions, under the fast switching assumption, see section 3. In our LGCA, the probability of a cell to jump from one node r to a neighboring node r + c i is equal to the probability of a cell residing in the respective velocity channel. In the isotropic case, cells within velocity channels are redistributed with probability 1/b. To connect the discrete mechanism with a continuum model we consider the expected cell density of a node r at time k, defined by ρ(r, k) = (r, k) , in the scaling limit when time step length and lattice spacing go to 0. The change in average occupancy of a moving cell at site r during the next time step is given by where the terms can be interpreted as the probability that a cell at site r + c i , i = 1, . . . , b, is moving to site r and the probability that a moving cell at site r does not attempt to leave the site. In the following, for simplicity, we will analyze a one-dimensional situation where motility events only take place in the horizontal direction. This approach can easily be extended to higher dimensions, because of the model isotropy.
In one dimension, equation (20) reduces to where the terms correspond to the probability that a cell at site r − 1 is moving and moves to the right, the probability that a a cell at site r + 1 is moving and moves to the left and the probability that a moving cell at site r does not attempt to leave the site. Approximating r s (ρ) around ρ 1 such that r s (ρ) = r s (0) + r s (0)ρ + r s (0)ρ 2 + O(ρ 3 ), equation (21) becomes To convert the finite difference equation (22) into a continuous macroscopic partial differential equation (PDE), we identify the average density ρ(r, k) by its continuous counterpart ρ(x, t), where x = r ∈ R and t = kτ ∈ R + , with , τ ∈ R + , which gives Expanding all terms in (23) in truncated Taylor series in powers of and τ up to second order gives To proceed, we consider the diffusive limit, i.e. → 0, τ → 0 and lim ,τ →0 .
Finally, the jump process of our LGCA can be approximated by a degenerate diffusion equation 5 Kernel density estimation of the averaged cell density The kernel density estimation, also known as Parzen-Rosenblatt density estimation, is a well known nonparametric way of estimating the probability density function underlying a finite set of observations [1,2]. We use a built-in Mathematica function to estimate the kernel density of the averaged cell density. 6 Stability analysis of cell reaction mean-field equation In the following, we analyze the stability behavior of the macroscopic growth equation (19) derived in section 4.1. Thus, let us consider the kinetic equation where with positive constants r b and r d and Obviously, F (0) = 0 and F (1) = −r d < 0. The first derivative of F with respect to is where r s ( ) = 1 2 κ sech 2 (κ( − θ)). For the trivial fixed point ρ 1 = 0 of (27), it is In order to identify fixed points of (27) and to determine their switch behavior, we will consider two conditions tanh(κθ) > 1 − 2 and Lemma 1. The zero ρ 1 = 0 of F ( ) is a stable equilibrium of (27) if condition (31) holds. Under condition (32), ρ 1 = 0 is an unstable equilibrium of (27).
For r b > r d , we will consider the two cases κ < 0 (repulsive case) and κ > 0 (attractive case).
Remark 3. If, values of ρ exist for which r s (ρ)(1 − ρ) > r d r b , then (27) shows bistable behavior and population extinction can be observed for small densities if (27) holds.

Movie Legends
Movie S1. Population extinction in the LGCA go-or-grow model. Initially, in a radius of 10 from the center one rest and one velocity channel is randomly occupied (threshold 0.25).