Laws of Spatially Structured Population Dynamics on a Lattice

: We consider spatial population dynamics on a lattice, following a type of a contact (birth– death) stochastic process. We show that simple mathematical approximations for the density of cells can be obtained in a variety of scenarios. In the case of a homogeneous cell population, we derive the cellular density for a two-dimensional (2D) spatial lattice with an arbitrary number of neighbors, including the von Neumann, Moore, and hexagonal lattice. We then turn our attention to evolutionary dynamics, where mutant cells of different properties can be generated. For disadvantageous mutants, we derive an approximation for the equilibrium density representing the selection–mutation balance. For neutral and advantageous mutants, we show that simple scaling (power) laws for the numbers of mutants in expanding populations hold in 2D and 3D, under both ﬂat (planar) and range population expansion. These models have relevance for studies in ecology and evolutionary biology, as well as biomedical applications including the dynamics of drug-resistant mutants in cancer and bacterial bioﬁlms.


Introduction
Before we begin describing some of our attempts to derive a number of mathematical laws for biological population dynamics, one of the authors (N.L.K.) would like to express her eternal gratitude to Michael Tribelsky, who was her teacher in the Physics Department, Moscow State University, at the beginning of the 1990s. Without his advice and guidance, N.L.K. would not be what she is now. Tribelsky's relentless optimism (often disguised) has taught her to overcome whatever difficulties life has posed. N.L.K. always follows Tribelsky's principle that problems (in life, as well as in science) must be addressed at the same pace, and not faster, than they are thrown at you ("Проблемы нужно решать по мере их поступления"). Tribelsky's "special course", delivered to the theoreticians from the Low Temperature Division, opened up a universe of diverse phenomena that are, at the same time, fascinating and amenable to understanding. He helped instill the sense of wonder at the beauty of the surrounding world, and this has stayed with N.L.K. no matter what subject matter she happened to focus on, from language evolution to virus dynamics.
Much has been said about the existence of laws in biology. Physicists in particular often feel that some beautiful (and, hopefully, simple) formal relations must exist to help us navigate the complexities of life. In this paper, we look for universal biological laws in the area of population dynamics, which is an area relevant for studies in ecology and evolutionary biology. A rich literature exists about the evolutionary dynamics of mutant spread and invasion, investigating aspects such as fixation probabilities and fixation times of different kinds of mutants in constant populations [1,2], as well as mutant dynamics in growing populations [3][4][5], where mathematical approaches were motivated by the famous Luria-Delbruck experiments [6]. Traditionally, much of this work has focused on wellmixed (non-spatial), homogeneous populations. In a wide array of biological scenarios, however, cells and organisms evolve in more complex settings, such as in spatially structured habitats (bacterial biofilms, cells in tissues and tumors) and within heterogeneous population structures (such as stem and more differentiated cells in tissues [7] or bacteria with different degrees of specialization in biofilms [8,9]). In the last 15 years, theoretical work has extended our understanding of mutant emergence in spatially structured populations [10][11][12][13][14][15][16][17][18]. An excess of mutational jackpot events was observed in spatial compared to well-mixed systems; such events result from mutations arising at the surface of expanding, spatially structured populations, surfing at the edge of range expansions, and appearing as mutant "bubbles" or "slices". These jackpot events have been implicated in the finding that the average number of neutral mutants when the total population reaches a given threshold size is significantly larger in spatial compared to non-spatial settings [18]. Further work by our group [19] showed that the evolutionary dynamics of mutants in spatially structured, expanding populations are characterized by additional complexities.
In general, spatial population dynamics are significantly more complex compared to the mass-action dynamics, and analytical approximations are not easily derived. In this paper, we describe two examples where our group has been successful in obtaining analytical approximations for laws of spatial dynamics: (a) steady-state density in space; (b) scaling laws for the number of mutants in expanding populations.
(a) Analytical approximations for steady-state density. In [20], we considered the phenomenon of range expansion, the process in biology by which a species spreads to new areas. We derived a numerical methodology that allowed for efficient computations of the number of individuals as the species expands its range in space. In addition, we found approximations for the steady-state density of populations, or the core density of expanding colonies, in several cases, such as the von Neumann and Moore neighborhoods on a square lattice and for a honeycomb neighborhood on a hexagonal lattice; see Table 1. For the same death-to-birth ratio, the grid is more packed under the Moore neighborhood, because of the availability of more neighbors per site. As a consequence, the equilibrium density corresponding to the Moore neighborhood is higher and closer to that of massaction. A general formula that provides an approximation for the steady-state density as a function of the number of neighbors is also given ( Table 1). In this paper, we provide a novel mathematical derivation of these laws, which generalizes to any neighborhood size (Section 3). Table 1. Steady-state densities for different types of lattices, as well as the mass-action system and a generalized, N b -neighbor system; see [20]. Here, δ = D/L, the division-to-death rate ratio. Red dots denote the focal cell and the blue ones stand for its neighbors.

Lattice
Geometry Density von Neumann Laws of disadvantageous, neutral, and advantageous mutant spread in different geometries. In [19], we studied the evolutionary dynamics of disadvantageous, neutral, and advantageous mutants in spatially expanding populations, where the growing front was characterized by different symmetries in two dimensions (2D) and 3D; see Table 2. In particular, while disadvantageous mutants grow linearly with N (and differences among different cases are more subtle), neutral mutants grow as a decreasing power of N as we go from a 2D flat front (which is essentially a 1D growth), to a 2D and 3D range. There are always fewer neutral or advantageous mutants in the mass-action (exponentially growing) case compared to any spatially restricted growth. Here (Sections 4 and 5), we provide a derivation of these laws, which follows [19]. Table 2. Mutant scaling laws for expanding populations; see [19]. Here, N is the population size, u is the mutation rate, and α > 1. See text for details. In the case of disadvantageous mutants, the number of mutants always scales with the first power of N. In addition to this scaling law, for the case of disadvantageous mutants, we also derive an approximate expression for the number of mutants in quasi-steady-state, given that the number of wild-type cells is N:

2D
see Section 4.1 for details and Section 2 for definitions (this expression is valid for a 2D von Neumann grid, and can be generalized to other cases). We show that the proportion of mutants at selection-mutation balance is higher for spatially distributed systems compared to well-mixed systems at equilibrium.

General Setup and Agent-Based Modeling of Population Dynamics
In order to describe the spatial growth and turnover of cells, as well as the population dynamics of species, we use a continuous-time Markov process, which is a generalization of the usual birth-death process and a type of a contact process. We assume that individuals of one or two different types exist, which we call "wild-type" and "mutant" individuals (to reflect the versatility of these models, we use both "individuals" and "cells" to refer to the biological agents under consideration). At the core of the description is a lattice that specifies possible locations of cells. This can be viewed as a geometric network, where each node is connected with its neighbors. The state space consists of different locations of the cells of the two types on the lattice. Cells are characterized by division and death rate parameters, which are denoted by L w and D w for the division and death rates of the wild-type cells and L m and D m for division and death rates of mutants. In the case where only one type of cell is considered, the notation is simply L and D for the two rates.
Each of the nodes can either be empty or contain one cell of either type. During an infinitesimally small time increment, ∆t, a given wild-type cell can attempt a division with probability L w and death with probability D w . If division is attempted, then an offspring is placed in one of the neighboring locations, chosen randomly and uniformly; division only happens if the chosen node is empty. The offspring cell is wild-type with probability 1 − u and mutant with probability u, where 0 ≤ u < 1 is the mutation rate. Mutant cells divide and die according to similar rules, except in the models considered here, we always assume that an offspring of a mutant cell is a mutant cell.
A number of biological phenomena can be studied by slight modifications of this process. To study the growth laws in the absence [20] or presence [19] of mutants, we assume an infinitely large grid. To study the quasi-steady-state (which describes, e.g., the turnover of cells in homeostasis), we make the grid finite and impose relevant boundary conditions.
For numerical explorations of these processes, agent-based modeling (ABM) is often used; see, e.g., [21][22][23]. For example, consider a two-dimensional ABM on a square grid with individuals of two types. A spot on the grid can be empty or can contain a cell, which is either wild-type or mutant. At each time step, the grid is randomly sampled N times, where N is the total number of cells currently in the system. If the sampled cell is wild-type, the cell attempts division (described below) with a probability proportional to L w or dies with a probability proportional to D w . When reproduction is attempted, a target spot is chosen randomly among the N b nearest neighbors of the cell. A neighborhood may contain, e.g., N b = 4 cells (the von Neumann neighborhood) or N b = 8 neighbors (the Moore neighborhood). If that spot is empty, the offspring cell is placed there. If it is already filled, the division event is aborted. (This modeling choice represents the assumption that the probability of divisions is reduced under more crowded conditions. This is similar to the logistic growth term often used in deterministic models.) The offspring cell is assigned wild-type with probability 1 − u, and it is a mutant with probability u. If the sampled spot contains a mutant cell, the same processes occur. Attempted division occurs with a probability ∝ L m , and the cell dies with a probability ∝ D m . Initial and boundary conditions are determined by the specific geometric setting investigated. For 2D spatial simulations, an n × n square or an n × w rectangular domain could be considered. At the boundaries of the domain, a spot is assumed to have fewer neighbors, i.e., more division events will fail. The process starts with a single wild-type cell, placed, e.g., in the center of the grid. Simulations always stop before the boundary of the grid is reached.
In what follows below, we derive some approximations of important observables from these types of dynamics, which have a clear biological meaning. We discuss both single-type populations and the co-dynamics of wild-type and mutant individuals.

Analytical Approximations for Steady-State Density (A Single Type)
In this Section, we consider populations consisting of only a single cell type with division and death rates L and D < L and no mutations (u = 0).

Equilibrium Density in Space: An Analytical Method
The process of range expansion (colonization) is one of the basic types of biological dynamics, whereby a species grows and spreads outwards, occupying new territories. Spatial modeling of this process is naturally implemented as a stochastic ABM of the type described above, with individuals (in this case, of only a single type) occupying nodes on a rectangular grid, births and deaths occurring probabilistically, and individuals only reproducing onto un-occupied neighboring spots. This approach is known to be computationally expensive. In [20], we derived a set of efficient computational tools, which were shown to be in good agreement with the underlying stochastic process of spatial expansion. As part of the method development, we were able to obtain approximate expressions for the quasi-steady-state (core) density of the individuals for different types of grids. In [20], we provided the density formula for the contact process: where δ = D/L and N b is the number of neighbors. Equation (2) was derived for several cases of N b , but no general derivation that would work for a given N b was supplied. Here, we present such a derivation, starting with N b = 2 and generalizing to any N b .
In 1D, with only N b = 2 neighbors, Equation (2) gives which, although a slight overestimation of the density, provides a good approximation. To derive this formula, consider a 1D ABM of the type described above. At the equilibrium, denote by p 1 the probability that a cell is located next to an existing cell and by p 0 the probability that a cell is located next to an existing empty spot. These two quantities can be estimated numerically. A 1D realization of the population at a fixed moment of time can be viewed as a Markov chain with states {0, 1}, where state 0 corresponds to the absence of a cell at a given spot and state 1 denotes the presence of a cell at a given spot. The transition probability matrix is given by where P 11 is the probability that there is a cell on the right of a cell, P 12 is the probability that the spot on the right of a cell is empty, P 21 is the probability that there is a cell on the right of an empty spot, and P 22 is the probability that an empty spot is on the right of an empty spot. The steady-state probability distribution of this process is (ρ, 1 − ρ), where ρ is the probability that a given spot contains a cell. This can be found as the eigenvector corresponding to the unit eigenvalue and is given by This quantity (given the numerically calculated p 0 and p 1 ) is a very good approximation of the actual (numerical) density. Another way to derive this connection between ρ and p 0 , p 1 is as follows: where the left-hand side is the probability to have a cell at a given point, the first term on the right assumes that there is a cell to the left of that point (ρ) and, then, the given point contains a cell with probability p 1 , and the second term on the right assumes that there is no cell to the left of that point (1 − ρ) and, then the given point contains a cell with probability p 0 . Solving this for ρ gives expression (4). Now, suppose that a cell is located at a given location. Then, a cell is located to the right of it (which we call the focal location) with probability p 1 . We have, after a single update of the contact process: The first term on the right assumes that there was no cell at the focal location (1 − p 1 ), but a cell on its left was chosen (probability 1/N) that divided to its right (L/2) or that a cell on its right exists (probability p 0 ), was chosen (1/N), and divided to its left (L/2). The second term assumes that there was a cell at the focal location (p 1 ) and it did not die (1 − D/N). The equality follows from the assumption of having an equilibrium.
Similarly, we can assume that there is no cell at a given location, then to its right (the focal location), the probability to have a cell, p 0 , satisfies: The first term on the right assumes that there was no cell at the focal location (1 − p 0 ), but a cell on its right exists (probability p 0 ), was chosen (1/N), and divided to its left (L/2).
The second term assumes that there was a cell at the focal location (p 0 ) and it did not die (1 − D/N).
Solving Equations (6) and (7) for p 0 and p 1 , we obtain: Substituting this into Equation (4), we obtain expression (3). This analysis easily generalizes to other systems (including higher dimensionalities), where the number of neighbors is given by N b . Instead of Equations (6) and (7), we have Solving this system, we obtain Substitution into Equation (4) (which holds in these systems because Equation (5) remains valid in these systems) yields Equation (2). Figure 1 plots the quasi-equilibrium density approximated by Equation (2) as a function of D/L for two finite values of N b corresponding to the von Neumann and Moore neighborhoods. Figure 1a compares them with the well-known mass-action result, showing that higher connectivity of the underlying network corresponds to a higher density of the cells. Figure 1b and Figure 1c compare the results with numerical simulations, demonstrating that Equation (2) provides a very good approximation. We notice that the approximation works better at higher densities, compared to lower ones. In the derivation, we assumed that the state of a given position on the lattice only depends on its immediate neighbor, and not on other, more distant, sites. This assumption becomes less valid at low densities, because macroscopic structures with strong correlations over distances > 1 form. Therefore, the approximation becomes worse as the death rate approaches the division rate.

Disadvantageous Mutants and Selection-Mutation Balance in Spatial Models
In the remainder of the paper, we consider two types of individuals, wild-type and mutants. This Section deals with disadvantageous mutants (defined below), while Section 5 focuses on neutral and advantageous mutants. We start with an Ordinary Differential Equation (ODE) formulation.

A Basic ODE Formulation
Let us denote the wild-type population as x(t) and the mutant population as y(t). Suppose that mutations happen at the rate u, and D w < L w . Then, the competition dynamics of cells in a mass-action system can be formulated as follows: where K is the carrying capacity. This system illustrates differences in the behavior of advantageous and disadvantageous mutants. If , that is the mutant excludes the wild-type and takes over. In this case, we say that the mutant is advantageous. If the inequality above is reversed, the mutants are disadvantageous and the so-called selectionmutation balance is reached, where mutants and wild-type cells coexist. This happens because the (stronger) wild-type cells produce mutants at a nonzero rate u, but the latter cannot take over and remains at a fraction of the population. For simplicity, let us assume that the mutation rate is low: Then, we can say that the mutants are disadvantageous if and then, the equilibrium solution is given bȳ where we neglected terms of the order u/γ, while the number of mutants is where we neglected terms of the order (u/γ) 2 . Below, we calculate the equilibrium densities of disadvantageous mutants and (advantageous) wild-type cells in a spatially distributed system at steady-state. This will also correspond to the densities in the core of an expanding system away from the advancing front.

A Spatial Description: Equations for the Densities
We restrict our description to a 2D square grid, with the von Neumann neighborhood (that is, each location has four nearest neighbors); the methodology is generalizable to the Moore neighborhood (eight neighbors). We use a method similar to that of [20]. Two random variables describe the state of the stochastic system at each spatial location, x: ρ x describes wild-type cells, such that and η x describes mutant cells, such that Note that ρ x and η x cannot be equal to one simultaneously; an empty spot corresponds to ρ x = η x = 0. We assume that wild-type cells have division and death rates L w and D w and mutant cells have division and death rates L m and D m . Wild-type cells mutate with probability u, and no back mutations are considered.
Denote the expectation of ρ x and η x by where we assumed that the expected values do not depend on spatial location, since we are interested in spatially homogeneous equilibrium solutions. We havė where the product x is empty, and the summation goes over all the neighbors of point x, which reproduce into location x at rates L w /N b and L m /N b if they are wild-type or mutant, respectively. The superscript in the notation ρ (k) x , η (k) x refers to the quantity at a location, k, neighboring the focal location, x.
Let us consider the von Neumann neighborhood (N b = 4). In the right-hand side of Equation (15), the terms are the summation having the form and in Equation (16), there are also terms of the form In the expressions above, we have ρ x ρ x is zero at location x (k) , and the three types of dyads are defined as follows: is the probability to have two wild-type cells at two neighboring locations; is the probability to have a wild-type cell and a mutant at two neighboring locations; is the probability to have two mutant cells at two neighboring locations. Figure 2a illustrates these three configurations. In terms of these correlations, Equations (15) and (16) can be rewritten aṡ The correlations for the three dyads that appear in these equations require their own equations to close the system. Let us derive an equation for W. We havė where we assume that one of the points in the dyad contains a wild-type cell (term ρ (k) x ), while the other point is empty (term (1 − ρ x )(1 − η x )), and that one of its neighbors (location x (j) ) contains a wild-type cell, which reproduces faithfully into point x at rate L w (1 − u)/N b . Note that either of the two points could be empty, which results in the multiplier of two in the first term on the right-hand side. Similarly, either of the dyad's locations can experience cell death, resulting in the negative rate 2D w . In order to calculate the average, we need to consider terms Note that here and below, the operation of averaging makes the expression independent of the actual location x. Further, the superscripts (k) and (j) do not refer to any specific neighbor of x, but to any neighbor of x; in particular, location x (j) may be the same as or different than location x (k) . In the case when the two locations are different, the correlation (21) is presented in Figure 2b, on the left. In the equations forṀ andİ, the following expressions appear in addition to Equation (21): These correlations are shown in Figure 2b, center and right. Therefore, denoting by a and b either ρ or η, we evaluate the average of the form which corresponds to a dyad with one of the locations (location x) empty and the other (location x (k) ), containing type "a", while a neighbor of x (location x (j) ) contains type "b". First, let us assume that location x (j) is different from location x (k) . Under von Neumann neighborhoods, this implies that x (j) and x (k) are not each other's neighbors because on a square grid, there could not be a non-degenerate triangle with diameter 1 or less. expression (22) is equal to The expression in the denominator is calculated as follows: Depending on the types at location x, the expressions in the numerator of Equation (23) can be of two types: and they are calculated in (17) and (18), respectively.
Next, we assume that location x (j) is the same as x (k) . Then, if types "a" and "b" in expression (22) are different, then we obtain x η (k) x = 0. If the types are the same, then we obtain expression (17) or (18). To summarize, expressions of type (22) are given as follows: The equation for W is then given bẏ Similarly, the other two equations can be derived: The closed system of equations for ρ, η, W, I, and M is given by Equations (19), (20), and (24)-(26).

Selection-Mutation Balance Solution
Solving these equations in the steady-state exactly is difficult, but if the mutation rate u 1, we can find the approximate solution. We start by setting u = 0 and obtaining the steady-state solution. Apart from the trivial solution and a negative solution, there are two symmetric solutions where only one type survives (competitive exclusion). We use the one where the wild-type excludes mutants: where the superscript corresponds to the zeroth order in the expansion in terms of small u. Note that, as expected, the expression for ρ (0) coincides with the one given in Table 1 for von Neumann grid (N b = 4). We then look for the first correction by substituting inserting in the system of five equations, keeping only the first order of expansion in u, and solving for ρ (1) , . . . , M (1) . We obtain This is the equilibrium solution corresponding to mutation-selection balance in the presence of spatial interactions. This approach is valid as long as the wild-type is advantageous (inequality (12)). In the opposite scenario, this solution is unstable, and the system converges to the mutants excluding the wild-type.
Under selection-mutation balance, of interest is the equilibrium proportion of mutants in the system given by It is instructive to compare this quantity with the equilibrium proportion of mutants in a mass-action system, ν eq m−a =ȳ/x, Equation (14): In other words, the relative contents of mutants (in proportion to the wild-types) are higher in spatially distributed systems compared to mass-action systems at equilibrium. This result resembles our recent analysis presented in [24], where we showed that the mean number of disadvantageous mutants is higher in fragmented populations compared to well-mixed populations. A comparison with non-equilibrium, growing, well-mixed populations is presented in [19] and not discussed here.

Applications and Comparison with Computations
The result in Equation (29) is applicable to two relevant scenarios. • Quasi-equilibrium density in finite populations. If simulations are continued until a finite grid is filled, the population reaches a quasi-equilibrium state where wild-type and disadvantageous mutant cells coexist. For a 2D square grid under a von Neumann neighborhood, the density of the mutants is approximated by η in Equation (28), while the density of wild-types is given by ρ (0) , Equation (27). The total numbers of mutant and wild-type cells are obtained by multiplying the quantities η and ρ (0) by the total number of grid points, respectively. Note that this scenario is not interesting in the case of advantageous or neutral mutants, as the entire population will eventually consist of mutant cells. • Number of mutants in spatially expanding populations. Simulating a growing population (on a large grid where the boundaries are not reached), we can ask how the number of disadvantageous mutants scales with the total number of cells, N. Since the core of the expanding colony is in quasi-equilibrium, Equation (29) shows that the number of mutants grows as the first power of N. Results for the scaling laws for neutral and advantageous mutants are discussed in the next Section.
The expected number of mutants predicted theoretically was compared with the results of numerical simulations. This was done in the following way. At size N, the number of mutants (in the von Neumann case) is predicted to be Nuν eq vN ; see Equation (29). Solving the equation Nuν eq vN =const, we can obtain the pairs (L m , D m ) of mutant kinetic rates corresponding to a predicted given number of mutants in a system of size N. Figure 3a shows the predicted number of mutant as a contour plot. The closer to the "neutrality" line (see inequality (12)), the larger the predicted number of mutants. The solution of equation Nuν eq vN = const is shown in Figure 3b and for five points from the solution set, the predicted number of mutants (given by Nuν eq vN = 10) is compared with the numerically obtained mean (plotted together with the standard deviation, Figure 3c. We can see that for larger mutant division rates, the deviation from the theory becomes significant. Figure 3d shows a histogram of the numbers of mutants for the parameters corresponding to the fifth point. One can see that the distribution has a long tail and a very large standard deviation. This is the consequence of macroscopic structures ("slices") that cannot be handled by the present, local, method.

Laws of Neutral and Advantageous Mutant Spread in Different Geometries
In the previous Section, we showed that (under the assumption of small mutation rates) the number of disadvantageous mutants in a spatially growing population scales with the total number of cells, as ∝ Nu. If cells are advantageous or disadvantageous, they obey different scaling laws [19]. In this Section we present simple derivations for those laws in different geometries.

Derivation of the Growth Laws
Below, we consider several scenarios that differ by the dimensionality of the grid (2D or 3D) and by the direction of growth (planar growth vs. range expansion; see Table 2.

Two-Dimensional Flat Front
Let us first assume that the death rate of cells is equal to zero. Consider cells growing along the surface of a cylinder of width W. This represents a one-directional growth process, where during each generation, we assume that W new cells appear, and the the total population is given by N = LW, where L represents the number of layers. The value of L is proportional to the number of generations, and thus to the physical time, t: The following calculation estimates the growth law of mutants. Every time a new layer (of width W) is added, the mean number of new mutations is given by Wu. Suppose that mutants are neutral. Then, each such mutation will give rise to an array of daughter mutant cells of width 1; see Figure 4. The length of this array is given by L − i, where i is the layer at which the mutation occurred. Therefore, the total expected number of neutral mutants is a cylinder of length L given by where we assumed L 1. Note that in this derivation, we assumed that the number of mutants is small compared to the total population, and individual mutant clones do not interact. In a more precise calculation, the number of wild-type cells in each layer is smaller than W because of the existence of mutants, and thus, the rate of new mutant production is smaller than Wu. We, however, assume that uLW 1, such that the number of mutants is relatively small.
Note that the number of neutral mutants decreases with W; see Figure 5; the largest number of mutants is achieved in the case of W = 1, a one-dimensional expanding array of cells.  Next, let us consider advantageous mutants. In this case, each new mutant gives rise to a triangular clone. In the first layer, the width of the clone is 1; in the next layer it is 1 + s; in the kth layer, it is 1 + (k − 1)s, where parameter s ≥ 0 measures the advantage of the mutants (with s = 0 corresponding to neutral mutants). Therefore, we have where for the approximation, we assumed that Ls 1. Furthermore, for this simple calculation to be valid, we need to assume that the wedges created by mutants do not come close to the cylinder's width, W, that is Ls W. In particular, Equation (31) can be valid for small values of W > 1, but only for mutants that are neutral for practical purposes (s 1). Note that when N is fixed, the total number of cell divisions that the system has undergone is also fixed. The number of mutants however is vastly different depending on the spatial configuration. It is the highest for W = 1 (one row of cells) and decreases drastically with the width of the cylinder. This is consistent with the notion that spatial restrictions result in a heightened number of mutants, the 1D space (W = 1) being the most spatially restrictive system. The reason for this is that in 1D, a mutant, once created, blocks the whole range of expansion and prevents wild-type cells from reproducing. The wider the front, the weaker this effect. Further, we note that in the special case where W = 1, mutant advantage does not play a role, and the number of advantageous, neutral, and even disadvantageous mutants is given by the same formula, Equation (30).
In the derivations above, a zero death rate of cells was assumed. This means that the colony spreads as a solid mass, where all of the spots in the core are occupied. Including a nonzero death rate does not change the geometric argument presented here, because the only difference now is that the expanding population is "porous", such that the same number of cells occupies a larger number of spots on the grid. Therefore, adding a nonzero death rate does not alter the scaling laws derived here and in the other cases. For this reason, we present calculations assuming a zero death rate. Numerical calculations in Section 5.2 confirm that the scaling laws remain the same in the presence of nonzero death rates.

Two-Dimensional: Circular Range Expansion
Next, we turn to the dynamics of neutral mutants on a circle. Let us suppose that the radius of the circle is R and N = πR 2 . The size of the colony increases via surface growth with N ∝ t 2 and R ∝ t.
As the range expansion proceeds, the circular layer of radius r will on average give rise to 2πru new mutations. Each mutation will result in a wedge expanding outwards. If the new mutation occurred in the layer with radius r, the number of mutating cells in layer r is 1. The number of mutants in the next layer is given by r+1 r , because under the assumption of mutant neutrality, the fraction of mutants in each new layer of radius j > r (with surface 2πj) should stay constant and equal to 1 2πr . For layer j, the number of mutants is then given by j/r. This gives rise to the following calculation: (the approximation is valid for R 1). For advantageous mutants in a growing 2D circle, the fraction of mutants will grow with each layer: where we assumed Rs 1. For this approximation to be valid, the mutant wedges should not exceed the circumference of the colony. Strictly speaking, this results in the condition Rs << 2πR, that is s 1. For larger values of s, the events where the mutant covers the whole surface of the colony are no longer negligible.

Three-Dimensional Flat Front
In a 3D space, let us first consider a solid cylinder of constant radius R 0 , where initially, the cells are situated as a layer at the bottom of the cylinder and proceed to grow by adding layers of size πR 2 0 . Each generation contributes πR 2 0 u new mutants, and as the colony grows to length L (and volume 2πR 2 0 L), we have in the neutral case: which is similar to the 2D flat front expansion. If the mutants are advantageous, then their number will increase from layer to layer, giving rise to conical wedges. This gives rise to the following calculation: where Ls 1 for the approximation, and the approach is valid as long as the wedge radius is smaller than that of the cylinder, sL R 0 .

Three-Dimensional Range Expansion
Next, we consider a 3D expanding sphere. For a sphere of radius R, we have N = 4/3πR 3 , and the surface is given by 4πR 2 . The size of the colony increases via 3D surface growth with N ∝ t 3 . Each spherical layer of radius r will on average give rise to 4πr 2 u new mutations. Each mutation will result in a conical wedge expanding outwards. If the new mutation occurred in layer with radius r, the number of mutating cells in layer r is 1. The number of mutants in a layer of radius j > r is given by (j/r) 2 , because under the assumption of mutant neutrality, the fraction of mutants in each new layer should stay constant (and equal to 1 4πr 2 ). Therefore, we write: (the approximation is again valid for R 1). If the mutant in a growing 3D sphere is advantageous, the fraction in each layer will increase according to the fitness advantage s and stretch from layer to layer in the same way as for the neutral mutants. We therefore have As before, the approximation holds if Rs 1. The method assumes that the mutant colony's size in each layer does not come close to the surface area, which amounts to the inequality s 1.

Exponential (Non-Spatial, Mass-Action) Growth
Finally, for exponentially growing populations, similar formulas could be derived. In particular, for neutral mutants, we have and for advantageous mutants with advantage α (which is the ratio of the net growth rate of mutants and the net growth rate of wild-type cells), we have see [25], Equation (14c); see also Equation (13) there for a more general formula.

Comparison with Numerical Simulations
We ran numerical simulations to check the results derived here. Figures 6 and 7, which are described in detail below, illustrate the accuracy and applicability of the approximations derived above. They contain colored lines, which are the results of the numerical simulations, and black "guides for the eye", which represent power-law functions characterized by the powers predicted by the theory. The figures show that these scaling laws hold over large intervals of N, the population size.   Table 3. The rest of the parameters are u = 5 × 10 −5 , W = 100 (except G, where W = 1000). Figure 6 shows results for the case of 2D systems, with flat front expansion presented in Figure 6a and range expansion in Figure 6b. Plotted are the average numbers of mutants (the vertical axis) measured at different population sizes, N (the horizontal axis). Different curves correspond to different parameter values, summarized in Table 3.
To simulate flat front expansion (Figure 6a), ABM simulations were performed in cylindrical geometry, with an n × W rectangular domain of width W. We started with an array of W wild-type cells at the left boundary of the domain and imposed periodic boundary conditions in the transversal direction. In each simulation, the cell population was allowed to grow to a size N, and the number of mutant cells at this size was recorded. Such simulations were performed repeatedly, and the average number of mutants was calculated. Simulation runs, in which the total cell population went extinct due to stochastic effects, were ignored.  Table 4.  Figure 6a represent neutral mutants in the absence (A) and in the presence (B) of cell death. The black solid lines are guides for the eye with slope 2 in the log-log plot, representing the quadratic scaling law (see Table 2). Curves (C-G) represent advantageous mutants, and the dashed lines are guides for the eye with slope 3 in the log-log plot, representing the cubic scaling law (see Table 2). The different cases of advantageous mutants include systems with and without cell death, cases where mutants are advantageous by divisions (that is, have a larger division rate and the same death rate, compared with wild-type cells), and cases where mutants are advantageous by deaths (that is, have a smaller death rate and the same division rate, compared with wild-type cells). The cubic scaling law holds at least for part of the N values for all these cases; see [19] for more details.

Curves (A) and (B) in
To simulate range expansion (Figure 6b), we performed simulation on a square grid, starting with a single cell in the middle and letting the population expand outwards. Simulations resulting in population extinction were discarded, and all simulations were stopped before the boundary of the grid was reached. In Figure 6b, curves (A) and (B) again represent neutral mutants without and with cell death (with solid guides to the eye having slope 3/2). The rest of the curves again explore advantageous mutants under different assumptions, with dashed guides for the eye having slope 2. In all cases, the scaling laws in Table 2 are confirmed.
Finally, Figure 7 demonstrates numerical results for 3D expanding colonies (with parameters listed in Table 4). The different cases considered include neutral and advantageous mutants experiencing 3D flat and range expansion, with the advantage realized through differences in division and death rates.

Conclusions
In this paper, we reviewed some recent results on the behavior of populations in lattice models. Both homogeneous populations (that is, populations consisting of a single type of individuals) and evolving populations (wild-type and mutants) were considered, and some simple laws derived approximate population densities in different dimensionalities, geometries of growth, and lattice types. The methodologies developed here can be extended to other cases, for example Equation (1) of mutant density was derived for a 2D von Neumann grid, but the methodology can be generalized, e.g., to the Moore lattice and also to 3D systems.
Our results have further practical applications, for example for understanding the dynamics of drug-resistant mutants in solid tumors [26] or the laws of evolution in bacterial biofilms [27].
The laws described in this paper can be viewed as a way to reduce a complex (spatial, stochastic) process to a small number of important observables, such as the mean equilibrium population density or the expected number of mutants. These quantities are shown to obey some simple rules, which relate these observables to the microscopic (kinetic) parameters of the cellular turnover and to system geometry. Having these simple rules can be very useful, if one, for example, fits a stochastic agent-based modeling (ABM) to a set of data: some parameters can be extracted by applying these laws and solving for the unknown quantities.
Aspects of the theory presented here can be tested in biological experiments. Even though lattice models present a certain simplification of reality, the scaling laws of mutant growth ( Table 2) are quite versatile, as they have been shown (numerically) to hold in models with different neighborhood structures and over a large variety of assumptions on cells' kinetic parameters. While some of these laws were tested previously (see, e.g., [18]), others, such as 3D flat and 3D range expansion laws, have not yet been confirmed experimentally.
Having simple laws that connect the macroscopic state with microscopic variables can be useful in interpreting experimental results. The laws in Table 1 describe the equilibrium density of populations of cells undergoing a turnover. By comparing the observed population density in experiments with the theoretical predictions, one could simply deduce the division-to-death ratio of cells growing in a given system. This can provide a useful quantitative measure of cells' state. In particular, the differences in cellular density can inform one of the changes in kinetic parameters that arise as a consequence of changes in the microenvironment. Examples of relevant experiments include investigations of the effect of various drug treatments on the cellular turnover.
Finally, our techniques can be extended to describing more complex systems. One potential application is deriving spatial laws of the microbial colony dynamics, such as bacteria growing in a microfluidic trap [28][29][30]. Deriving rules of global dynamics of microbial communities from local interaction rules is a task that is similar in spirit to the efforts that were described in the present paper.