Genetic drift and selection in many-allele range expansions

We experimentally and numerically investigate the evolutionary dynamics of four competing strains of E. coli with differing expansion velocities in radially expanding colonies. We compare experimental measurements of the average fraction, correlation functions between strains, and the relative rates of genetic domain wall annihilations and coalescences to simulations modeling the population as a one-dimensional ring of annihilating and coalescing random walkers with deterministic biases due to selection. The simulations reveal that the evolutionary dynamics can be collapsed onto master curves governed by three essential parameters: (1) an expansion length beyond which selection dominates over genetic drift; (2) a characteristic angular correlation describing the size of genetic domains; and (3) a dimensionless constant quantifying the interplay between a colony’s curvature at the frontier and its selection length scale. We measure these parameters with a new technique that precisely measures small selective differences between spatially competing strains and show that our simulations accurately predict the dynamics without additional fitting. Our results suggest that the random walk model can act as a useful predictive tool for describing the evolutionary dynamics of range expansions composed of an arbitrary number of genotypes with different fitnesses.

Introduction A competition between stochastic and deterministic effects underlies evolution. In a wellmixed system such as a shaken culture of the yeast microorganism Saccharomyces cerevisiae, stochastic competition between individuals, mutations, and selection dictate the dynamics of the population [1]. In spatially structured environments, active or passive dispersal of individuals also plays an important role. The local "well-mixed" dynamics must be coupled to the motion of individuals, leading to strikingly different evolutionary dynamics, even in the absence of selection [2][3][4][5][6][7]. A model laboratory system that can be used to explore the coupling between local "wellmixed" effects and spatial deterministic and stochastic dynamics is a microbial range expansion [8], in which a population expands into an unoccupied region of a hard agar Petri dish. Non-motile microbes expand outwards from their initial position due to a combination of growth coupled with random pushing by neighboring cells and leave behind a record of their genetic competition as they cannot move and cease reproducing once the population becomes too dense [8]. A frozen genetic pattern of four competing strains of E. coli marked by different fluorescent colors can be seen in Fig 1. Spatial structure is present in the frozen genetic patterns because the microbes at the expanding frontier produce daughter cells of the same color that migrate only a small fraction of the front circumference within a generation. Hallatschek et al. [8] identified the key role of genetic drift in producing these sectored patterns; the small population size at the front of an expanding population [9,10] enhances number fluctuations (i.e. genetic drift), eventually leading to the local fixation of one strain past a critical expansion radius R 0 . The decrease in genetic diversity as the small number of individuals at the frontier expands is referred to as the "Founder effect" [11].
Outside of the laboratory, range expansions occur naturally during the spread of invasive species such as the bank vole in Ireland [12] or the cane toad in Australia [13], and played a role in the evolutionary history of humans when migrating out of Africa [14]. In these natural expansions, populations may have many competing genotypes, or alleles, each instilling a different fitness. Even if a population is originally clonal, mutations may create new alleles that compete with one another to proliferate, a phenomenon known as clonal interference [15].
An allele's fitness is often determined by its corresponding expansion velocity. Faster expanding individuals will colonize more territory and will block slower strains from expanding, resulting in the increased abundance of 'faster' alleles at the frontier [13,16,17]. If the curvature of a microbial colony can be neglected and its front is sufficiently smooth, it has been shown both theoretically and experimentally that the domain wall of a faster expanding strain will displace a slower expanding strain at a constant rate per length expanded after an initial transient, resulting in a characteristic triangular shape [17] as shown on the right side of Fig 1. If the curvature of the expansion is not negligible, the sector boundaries will trace logarithmic spirals [17].
Even in the most simple scenario when de-novo mutations and mutualistic or antagonistic interactions are ignored, the dynamics of many competing alleles with varying fitnesses at the front of a range expansion have neither been quantified theoretically nor explored in laboratory experiments. Prior laboratory experiments focused on the dynamics of a single sector of a more fit strain (representing a competing alelle) of yeast sweeping through a less fit strain [17] in regimes where stochastic wandering of genetic boundaries was not expected to be important. Recent experimental work studied how fast a single more fit strain swept through a less fit strain in a range expansion and compared the dynamics to the same strains in a well mixed test tube [9].
In this paper, we experimentally and numerically investigate the dynamics of four competing strains (alleles) of E. coli with varying selective advantages initially distributed randomly at the front of a radial range expansion. The eCFP (blue) and eYFP-labeled (yellow) strains expanded the fastest, followed by the non-fluorescent (black) strain, and finally the mCherrylabeled (red) strain. The differences in expansion speeds are reflected in Fig 1 as follows: the yellow/blue bulges at the front of the expansion are larger than the black bulges which are larger than the red bulges. The significant random undulations at the frontier, however, significantly mask the selection-induced bulges. Left: A four-color E. coli range expansion. Four strains of E. coli differing only by a heritable fluorescent marker were inoculated on an agar plate in a well-mixed droplet and expanded outwards, leaving behind a "frozen record" of their expansion. Virtually all growth occurred at the edge of the colony. The markers instilled different expansion velocities: our eCFP (blue) and eYFP (yellow) strains expanded the fastest, followed by our black strain, and finally our mCherry (red) strain. As a result of the differing expansion velocities, the yellow/blue bulges at the frontier are larger than the black bulges which are larger than the red bulges, although the significant stochastic undulations at the front mask their size. The microbes segregate into one color locally at a critical expansion radius R 0 due to extreme genetic drift at the frontier [8]. After segregated domains form, genetic domain walls diffuse and collide with neighboring walls in an "annihilation" or "coalescence" event indicated by an A or C, respectively. Right: Illustration of the relevant parameters used to model range expansions. Here, a faster expanding, more fit yellow strain with expansion velocity u Y is sweeping through a less fit red strain with expansion velocity u R , in a regime where the curvature of the colony can be neglected. The length expanded by the colony is L = R − R 0 . We characterize domain wall motion per differential length expanded dL and the wall's differential displacement perpendicular to the expansion direction dx. v YR w is a dimensionless speed, characterizing the yellow-red (YR) domain wall's average expansion dx per length expanded dL, i.e., v YR w ¼ dx=dL: D w is the domain walls' diffusion coefficient per length expanded; it controls how randomly the domain walls move. We treat the dynamics of our four strains as a one-dimensional line of annihilating and coalescing random walkers using the parameters R 0 , D w , and v ij w , where ij represents all possible domain wall types. https://doi.org/10.1371/journal.pcbi.1005866.g001 As is evident from Fig 1, the size and location of a monoclonal sector can be described by the locations of its boundaries. When two boundaries collide, they either annihilate if the neighbors to the left and right of the collision are the same or coalesce if the neighbors are different, as illustrated by the A and C respectively on the left side of Fig 1. We therefore describe our expansions as a one-dimensional line of annihilating and coalescing random walkers, a description that has been used extensively in previous work (see Ref. [2] for a review). To account for the radial geometry of our colonies, we allow the frontier to inflate, corresponding to the increasing perimeter of the colony as its radius increases. Past the radius R 0 where genetic domains originally form, we describe the random motion of genetic domains by a diffusion constant per length expanded D w (see Fig 1) [18]. If dx characterizes the displacement of a domain wall perpendicular to the expansion direction and dL is the distance the colony has expanded (the radius that the colony has grown) as illustrated on the right side of Fig 1, where we neglect the circumferential curvature in this small region, we define the diffusion constant per length expanded as 2D w = dVar(x)/dL where Var(x) hx 2 i − hxi 2 is the variance and the brackets indicate an average over many domain walls. Note that D w has dimensions of length. Similarly, differences in expansion velocities between neighboring strains will lead to the deterministic displacement of domain walls per length expanded as the faster expanding strain will reach the contested point on the front before a slower growing strain as mentioned above [17]; we characterize this deterministic motion by a dimensionless "wall velocity," [18] v ij w ¼ dhxi=dL, where i is the strain to the left of the domain wall and j is the strain to the right. Note that v ij w ¼ À v ji w . The dynamics of an arbitrary number of neutral competing strains in an expansion (i.e. v ij w ¼ 0 for all domain walls) is well understood as the dynamics can be described as a onedimensional system of annihilating and coalescing random walkers [19][20][21] which is equivalent to a one-dimensional q-state Potts model [22,23] governed by zero-temperature Glauber dynamics [24] or a q-opinion Voter model [25,26]. Many theoretical predictions and analyses of this system exist; of particular relevance to this paper are the relative annihilation and coalescence rates per collision as q is varied [27][28][29] and the calculation of spatial correlation functions [28]. To map standard linear results onto an inflating ring (i.e. including R 0 in the models), one can use a conformal time transformation [30][31][32]. Fewer results are available in the presence of selection, i.e. when domain walls have deterministic biases (nonzero v ij w ) [33]. Analytical results are rare because the moment hierarchy of this model does not close [2] as discussed in S1 Appendix.
In this paper, we measure and predict three quantities relevant to the evolutionary dynamics of our four competing strains of E. coli in radial range expansions: the average fraction of each of our four strains, the two-point correlation functions between our strains, and the relative annihilation and coalescence probabilities per domain wall collision (see Fig 1), a quantity that has received theoretical attention [27][28][29] but has neither been explored experimentally nor investigated in the presence of selection. We measure these three quantities using an image analysis toolkit (available on GitHub, complete with examples of how to use it [34]) that extends experimental techniques for two-color (two-allele) range expansions [8,9,17,18,35,36] to an arbitrary number of competing strains. We next use an efficient radial simulation (also on GitHub [34]) of annihilating and coalescing random walkers with deterministic wall velocities to determine what sets the scale of the dynamics and to synthesize our experimental and theoretical results. We show that three key combinations of R 0 , D w and v ij w control the dynamics of our four strains. We conclude with suggestions for future studies. The details of our experimental, theoretical, and simulation methods are given in the last section.

Experimental results
We begin by reporting our measurements of the average fraction of each strain, the two-point correlation functions between strains, and the relative rates of annihilations and coalescences as a function of length expanded for our four competing strains of E. coli. As discussed in the Materials and Methods, we found that our eCFP and eYFP strains had the fastest expansion velocities followed by the black strain and finally the mCherry strain (see Table 1). We expected that our experimental measurements would reflect this hierarchy of speeds; faster expanding strains should have a larger fitness than slower expanding ones. To illustrate the presence of selection, we used neutral theory (discussed in detail in S1 Appendix) as a null expectation; selection caused deviations from the neutral predictions. To calibrate neutral theory to our experiments we fit R 0 and D w , two model parameters illustrated in Fig 1, following the procedures discussed in the Materials and Methods. The fit values of R 0 and D w can be seen in Table 2. In later sections, we show how to predict the average fraction, two-point correlation functions, and relative rates of annihilation and coalescences using our random-walk model and simulation. Table 1. The expansion velocity u i and each strain's selective advantage relative to the mCherry strain s iR = u i /u R − 1 measured over the course of seven days for radii greater than R 0 (the radius where distinguishable domain walls formed) averaged over three independent experiments conducted on separate sets of agar plates. s wm iR ¼ g i =g R À 1 was the fitness of each strain relative to mCherry in liquid culture with respect to their basal growth rates g i and g R . The radial expansion velocity fitness s iR did not match the well-mixed liquid-culture fitness s wm iR . However, every strain in liquid culture still grew faster than mCherry. Interestingly, the black strain grew faster than the eCFP and eYFP strains in liquid culture while on agar, the eCFP and eYFP strains expanded faster than the black strain. See the Materials and methods for additional information.

Strain
Average Velocity u i (mm/day)  Table 2. Parameters used in the annihilating and coalescing random-walk model. We experimentally measured R 0 , D w , and v ij w using the procedures outlined in the Materials and methods so that we could compare experimental results with our model's predictions.

Parameter
Value Description Average fractions. The average fraction of strain i at a length expanded of L = R − R 0 is defined as where f i (ϕ, L) is the local fraction of strain i at angle ϕ and length L (i.e. at a pixel specified in polar coordinates by ϕ and L). The angular brackets represent an average over many range expansions and f i is normalized such that ∑ i f i (ϕ, L) = 1 for each location in the colony as discussed in the Image Analysis section. In the neutral case, the average fraction of each strain should equal their inoculated fractions and should be independent of length expanded. Selection forces the average fractions of less fit strains to decrease. We measured the average fraction versus radial length expanded in two separate sets of experiments where we inoculated different fractions of our eYFP, eCFP, and mCherry strains. In one experiment, we inoculated the eYFP, eCFP, and mCherry strains with equal initial fractions of 33% while in the other we inoculated 80% of the mCherry strain and 10% each of the eCFP and eYFP strains. We conducted 20 replicates in each case and calculated the average fraction of each strain using our image analysis package. Fig 2 displays the trajectories of the 20 expansions and the mean trajectory (the average fraction) as ternary composition diagrams for both sets of initial conditions [37].
In both sets of experiments, we observed a systematic drift away from the mCherry vertex as a function of radius as illustrated by the mean trajectories shown as insets. We witnessed two cases where the 10% initial inoculant of the eCFP strain became extinct, represented by the pinning of trajectories to the absorbing boundary connecting the eYFP and mCherry vertex, a consequence of the strong genetic drift at the frontiers of our E. coli range expansions. These measurements indicate that the mCherry strain was less fit than the eCFP and eYFP strains, consistent with the order of the radial expansion velocities. Average fraction of each genotype as a function of length expanded for 20 radial expansions each when equal fractions of eCFP, eYFP, and mCherry were inoculated (left) and when 10% eCFP, 10% eYFP, and 80% mCherry were inoculated (right). The red dot indicates the composition at the radius R 0 = 3.50 mm where distinct domain walls form and the blue dot indicates the composition at the end of the experiment. The red dots are dispersed about the initial inoculated fractions due to the stochastic dynamics at the early stages of the range expansions when R < R 0 . The highly stochastic trajectories illustrate the importance of genetic drift at the frontier in the E. coli range expansions. The smaller ternary diagrams display the average fraction over all expansions vs. length expanded for each set of experiments. For both initial conditions, we see a small systematic drift away from the mCherry vertex indicating that the mCherry strain has a lower fitness, in agreement with the independent radial expansion velocities of each strain (see Table 1). Note that two replicates on the right resulted in the complete extinction of eCFP due to strong spatial diffusion, indicated by the trajectories pinned on the absorbing line connecting the eYFP and mCherry vertices. Two-point correlation functions. Next, we measured the two-point correlation functions given by where f i (ϕ, L) is again the local fraction of strain i at angle ϕ and expansion length L. F ij gives us the probability that strain i is located at an angular distance of ϕ away from strain j at a length expanded L. Note that F ij = F ji and F ij (ϕ) = F ij (−ϕ). Although the average fraction is constant in the neutral case, the two-point correlation functions broaden due to the coarsening of genetic domains [2]. Neutral q-color Voter models analytically predict the form of the twopoint correlation functions [2] as seen in equation (S1.3) in S1 Appendix. Deviations from neutral predictions are caused by selection. Analytical results describing these deviations are not available for reasons discussed in the S1 Appendix (the hierarchy of moments does not close); numerical simulations must be used to calculate the precise shape of the correlation functions as seen in second half of our Results section. Regardless, selectioninduced deviations can be understood in the limit of both large and small angular separations. For large angular separations, spatial correlations will be negligible; the two-point correlation functions will consequently factorize and plateau at the value F ij = F i F j where F i is the average fraction at length L from above. Therefore, in neutrality, the two-point correlation functions F ij should plateau at F ij ¼ F 0 i F 0 j , the product of the initial fractions inoculated of strains i and j (in neutrality, F i does not change). Selection can thus be identified by comparing the experimentally measured plateau value to the neutral prediction value. Furthermore, in the limit of zero angular separation, it is known that @ ϕ F ij measures the density of ij domain walls [2] (where i 6 ¼ j). In general, if strain i is less fit than the other strains, it will have fewer domain walls, decreasing the domain-wall density and thus the slope near ϕ = 0.
We measured the correlation functions between each pair of strains in three sets of experiments where we inoculated equal well-mixed fractions of the eCFP, eYFP, and black strains, then eCFP, eYFP, and mCherry, and then finally all four strains. We conducted 20 replicates of each experiment, measured all two-point correlation functions at the final radius of R = 10 mm corresponding to a length expanded of L = R − R 0 = 6.5 mm, and averaged the results. In Fig 3, we plotted the neutral correlation function prediction and compared it to the experimentally measured correlation functions.
The two-point correlation functions in the experiment between eCFP, eYFP, and the black strains (first column of Fig 3) are consistent with the order of radial expansion velocities (see Table 1). The correlation between the eCFP and eYFP strains plateaued at a higher value than the neutral prediction while the correlation between eCFP and black plateaued at a lower value, indicating that the eCFP and eYFP strains were more fit. The self-correlation for the black strain, F BB , also plateaued at a value below eCFP, eYFP, and the neutral prediction, further indicating that it had a smaller fitness. The self-correlation data was more noisy than the correlation between strains, however; we consistently found that correlations between strains were better at detecting fitness differences than self-correlations.
In contrast, combining eCFP, eYFP, and mCherry in one set of experiments and all four strains in another revealed that mCherry had a larger fitness defect. Correlation functions including mCherry always plateaued at a significantly smaller value than correlation functions excluding it. Furthermore, off-diagonal (bottom-row of Fig 3) correlation functions involving the mCherry strain had a smaller slope at zero angular separation, indicating that less mCherry domain walls were present and that the mCherry strain was less fit than the others. The twopoint correlation functions were thus consistent with the black strain having a small selective disadvantage relative to eCFP and eYFP and the mCherry strain having a larger disadvantage relative to all others.
Annihilation asymmetry. The last quantity we measured was the relative rate of annihilations and coalescences per domain wall collision; examples of annihilations and coalescences can be seen on the left side of Fig 1. Many theoretical results exist describing the neutral dynamics of annihilations and coalescences and they are summarized in S1 Appendix. To succinctly quantify the difference between the annihilation and coalescence probabilities per wall collision, we define the "annihilation asymmetry" ΔP(L) = P A (L) − P C (L) as the difference in probability of obtaining an annihilation versus a coalescence per collision at a distance expanded of L. If q neutral colors are inoculated in equal fractions, it can be shown that Note that in neutrality, the annihilation asymmetry ΔP is independent of the length expanded L; it depends only on the number of strains q inoculated in equal fractions. In the presence of selection, however, less fit strains should be squeezed out as the length expanded L increases, forcing q and thus ΔP to change.

mm) in three sets of experiments where we inoculated 20 replicates with equal fractions of our eCFP, eYFP, and black strains (left), then eCFP, eYFP, and mCherry (center), and finally all four strains (right).
The shaded regions in these plots indicate standard errors of the mean. Using the measured diffusion coefficient D w and initial radius where domain walls form R 0 (see Table 2), we also plot the theoretical neutral two-point correlation functions (black dashed line; see eq. (S1.3)). The colors of each plotted correlation function were chosen to correspond to their composite strain colors; for example, two-point correlation correlation functions associated with mCherry were red or were blended with red. The subscripts correspond to the color of each strain: C = eCFP, Y = eYFP, R = mCherry, and B = Black. As judged by the magnitude of the deviation from neutral predictions, the black strain has a small selective disadvantage relative to eCFP and eYFP and the mCherry strain has an even greater disadvantage, in agreement with the independent radial expansion velocities of each strain (see Table 1). https://doi.org/10.1371/journal.pcbi.1005866.g003 Genetic drift and selection in many-allele range expansions To gain insight into the behavior of ΔP, for the case of q neutral colors in equal proportions, we have lim q!1 ΔP(q) = −1 (only coalescences), ΔP(q = 3) = 0 (equal numbers of annihilations and coalescences), and ΔP(q = 2) = 1 (only annihilations). The quantity ΔP thus provides a simple way to characterize the annihilation/coalescence difference in a single curve that varies smoothly between −1 and 1 as 2 q < 1. In S1 Appendix we develop and discuss the case when strains are inoculated in non-equal proportions (see supplementary equations (S1.8)-(S1.10)); in that scenario, it is useful to define a "fractional q" by inverting eq (3) to read q = (3 + ΔP)/(1 + ΔP) (i.e. a fractional q can be evaluated for a given ΔP).
To experimentally quantify the annihilation asymmetry, we examined the average cumulative difference in annihilations and coalescences vs. the average cumulative number of domain wall collisions as colonies expanded; ΔP is given by the slope of this quantity and can be seen in Regardless of which strains were inoculated and their selective differences, our results were consistent with the neutral theory prediction in eq (3) for q = 2, q = 3, and q = 4 as judged by the overlap of the black dashed line with the shaded standard error of the mean in each case. ΔP appeared to be constant as a function of length. We also tested an initial condition where we inoculated strains in unequal proportions: we inoculated 10% of eCFP and eYFP and 80% of mCherry. This experiment again matched the neutral prediction of ΔP % 0.51 (and correspondingly q % 2.33) within error. Evidentally, as discussed in more detail below, certain observables like the average fraction and two-point correlation functions show stronger signatures of selection than others like the annihilation asymmetry.

Simulation results
In this section, we introduce three key combinations of our random walk model's input parameters R 0 , D w , and v ij w (see Fig 1) that control the evolutionary dynamics of our four competing E. coli strains. Using simulation, we show that we can utilize these key combinations to collapse the simulated evolutionary dynamics (focusing on the experimental quantities we measured above: the average fraction, two-point correlation function, and annihilation asymmetry) of an arbitrary number of competing strains in a range expansion.
Key parameters. What key combinations of the variables seen on the right side of Fig 1 govern the evolutionary dynamics of our competing strains? Our goal is to describe the dynamics as a function of length expanded by our colonies L = R − R 0 with R 0 the initial radius where domain walls form, the domain wall diffusion coefficient per length expanded D w (units of length), and all wall velocities per length expanded v ij w (dimensionless). The two-point correlation functions must include an additional independent variable: the angular distance ϕ between strains.
Investigating the width of a single sector of a more fit allele sweeping through a less fit allele, as illustrated on the right of Fig 1, reveals important parameter combinations (see S1 Appendix for additional details). In a linear expansion, the deterministic, selection-induced growth of a sector of genotype i sweeping through a less fit genotype j will scale as v ij w L while its diffusive growth will scale as ffiffiffiffiffiffiffiffi ffi D w L p . At short lengths expanded, diffusion will thus dominate deterministic growth, and at larger lengths selection will dominate diffusion. A crossover expansion length L ij s [2,9,32] beyond which selection dominates follows by equating the deterministic and diffusive growth, The factor of 2 in front of v ij w and 4 in front of D w arises because we are monitoring the distance between two domain walls (i.e. a sector); similar arguments can be applied to describe the motion of individual walls. It is worth noting that the chirality of sector boundaries reported in E. coli range expansions [18,36] would result in a wall velocity pointing in the same direction (left or right) for every domain wall. We can ignore this constant bias in our models because sectors will still expand at the same rate despite an additional superposition of all domain walls moving in a specific direction. To avoid complications arising from chirality in this paper, we focus on quantifying the growth of sectors, i.e. the distance between two domain walls, as opposed to tracking the motion of an individual domain wall whenever possible. L ij s is the characteristic length that the colony must expand in order for selection to dominate over diffusion for strain i sweeping through strain j and acts as the first key parameter.
Upon repeating this argument for domains on a radially inflating ring (see S1 Appendix), we identify L ij I [32,38] as the inflationary analog of L ij s : the expansion length beyond which selection dominates over diffusion, and find κ ij is a dimensionless prefactor that can be thought of as an "inflationary selective advantage" The slope of this plot gives the annihilation asymmetry ΔP. The shaded regions represent the standard error of the mean between many experiments. We use the notation C = eCFP, Y = eYFP, B = black, and R = mCherry. Despite the presence of selection, ΔP was consistent with the standard neutral theory prediction of eq (3) for q = 2, q = 3, and q = 4 (equal initial fractions of q strains), as judged by the overlap of the black dashed lines with the shaded areas in every case. We also explored an initial condition where we inoculated unequal fractions of three strains; we inoculated 10% of both eCFP and eYFP and 80% of mCherry. Our experiments agreed with the prediction of ΔP % 0.51, or an effective q % 2.33, from the neutral theory developed in supplementary equations (S1.8)-(S1.10).
controlling the expansion length at which selection dominates over diffusion and is given by Fig 5 illustrates the importance of κ ij ; it displays the ratio of the inflationary to the linear selection length scale L ij I =L ij s as a function of κ ij from the numerical solution of eq (5). We find that the ratio of the length scales has the asymptotic behavior Thus, if κ ij ) 1, inflation can be ignored (relative to selection and genetic drift), and the inflating selection length scale approaches the linear selection length scale. In contrast, if κ ij ( 1, the inflationary selection length will be many times larger than the linear selection length scale [32]. As κ ij becomes smaller, inflation and genetic drift dominate over selection for a larger length expanded. κ ij is the second key parameter describing the dynamics of our system. Note that in contrast to a linear expansion which just features competition between genetic drift and selection (captured by the quantity L ij s ), a radial inflation has three, separate, competing effects: genetic drift, selection, and inflation. κ ij quantifies the strength of selection relative to inflation and diffusion.
The third and final key parameter is the characteristic angular correlation length between selectively neutral genotypes. This parameter arises naturally when In contrast, if κ ( 1, the inflationary selection length scale L I will be many times larger than the linear selection length scale L s , indicating that selection will be weak compared to inflation and diffusion (but will ultimately dominate at very large lengths expanded). The three black points correspond to measurements of the κ ij that govern the dynamics of our competing strains; N stands for the two selectively neutral strains (eCFP and eYFP), B for black, and R for mCherry (red). See the Predicting experimental results with simulation section for more details. https://doi.org/10.1371/journal.pcbi.1005866.g005 Genetic drift and selection in many-allele range expansions analytically calculating the neutral two-point correlation functions from the Voter model (see eq. (S1.3)). The parameter also has an intuitive description. When moving into polar coordinates, the angular diffusion coefficient D ϕ is related to the standard linear domain wall diffusion coefficient by D ϕ = D w /R 2 . The characteristic scale for the radius is R 0 ; the angular diffusive growth of domains should consequently scale as . Note that this characteristic angular length does not depend on the total number of strains; it describes the diffusive coarsening of a single strain sector propagating through one or more other strains.
We have now identified the three key parameters that govern the evolutionary dynamics of our competing strains. L ij s ¼ D w =ðv ij w Þ 2 is the length that a linear expansion must grow in order for selection to dominate over diffusion for strain i sweeping through strain j, ( 1) or inflation (κ ij ) 1) may be neglected relative to other effects in radially inflating expansions, and sets the characteristic angular correlation length between selectively neutral genotypes. These key parameters are listed in Table 2.
Collapsing the evolutionary dynamics with the key parameters. We used simulations of annihilating and coalescing random walkers constrained to lie on the edge of an inflating ring with deterministic biases due to selection (see the Simulation methods section for additional details) to investigate the effect of the parameters R 0 , D w , and the set of all v ij w on the evolutionary dynamics of our competing strains. As we varied R 0 , D w , and v ij w , we calculated the average fraction of each strain, the two-point correlation functions between strains, and the relative rate of annihilations and coalescences per domain wall collision (the quantities we measured experimentally). We also investigated the role of the three key combinations of parameters L ij s , κ ij , and ϕ c for both linear and radial expansions.
We first simulated q = 3 competing strains where two neutral strains swept through a third less fit strain with wall velocity v w , similar to our experiments with two neutral strains (eCFP and eYFP) and the less fit mCherry strain. The three strains were numerically inoculated in equal proportions. Note that in this simulation, there was only one non-zero v w and consequently one L s ¼ D w =v 2 w and one k ¼ We varied v w from 10 −3 v w 10 −1 and N 0 from 10 2 N 0 10 5 (altering R 0 = N 0 a/(2π)) and computed the average fraction F of the less fit strain and the annihilation asymmetry ΔP. We found that both F and ΔP from simulations with identical κ, despite different values of R 0 and v w , collapsed if L, the length traveled was rescaled by L s as seen in We now consider the collapsed curves F(L/L s , κ) and ΔP(L/L s , κ) as a function of the parameter κ as seen in Fig 6. κ had a pronounced effect on both quantities. For κ ≳ 5 the dynamics of F and ΔP approached the dynamics of a linear expansion at all L/L s , illustrated by the bright pink line on the left and the bright pink dots on the right of Fig 6; the more fit strain swept so quickly through the less fit strain that the colony's radial expansion could be ignored. As κ decreased, the less fit strain was squeezed out more slowly due to the inflation of the frontier, resulting in slower transitions from q = 3 to q = 2 colors and consequently slower transitions from ΔP = 0 to ΔP = 1. For κ ( 1, ΔP barely shifted from 0 over the course of the simulation. Interestingly, ΔP peaked at a finite L/L s for small κ; it is not clear what causes this effect, but it may be related to the transition from linear to inflation-dominated dynamics as L increases. Additional simulations revealed that for expansions composed of many strains with different fitnesses (multiple v ij w ) and consequently various κ ij , all of our observables (F, ΔP, and F ij ) could again be collapsed onto a master curve by rescaling L by any one of the selection length scales (i.e. L=L ij s ) and by rescaling ϕ by ϕ c ; the set of κ ij specified the master curve. An example of a simulation with that exhibits collapsed dynamics for three κ ij can be seen in Supplementary S4 Fig. To summarize the results of this section, we found that we could collapse the average fraction F, annihilation asymmetry ΔP, and the two-point correlation functions F ij by DPðL; D w ; R 0 ; fv ij w gÞ ! DPðL=L ij s ; fk ij gÞ ð9Þ where the brackets indicate a set of variables parameterized by i and j (i.e. fv ij w g represents the set of all ij wall velocities). As long as L was rescaled by any selection length scale L ij s and ϕ was rescaled by the characteristic angular correlation length ϕ c , the set of {κ ij } completely dictated the evolutionary dynamics.

Predicting experimental results with simulations
A major goal of this paper is to test if the annihilating and coalescing random-walk model can predict the experimental evolutionary dynamics of our four competing strains (alleles) with different fitnesses (radial expansion velocities). To the best of our knowledge, analytical results for the random-walk model are unavailable (as discussed in S1 Appendix); we consequently Genetic drift and selection in many-allele range expansions used our simulations to predict the dynamics. In this section we quantify the three key parameter combinations for our experimental expansions and then use them to predict the evolutionary dynamics of all four of our competing E. coli strains in an independent experiment.
In the last section, we found that our simulation dynamics could be collapsed onto master curves for a fixed set of k ij ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi R 0 =L ij s q by rescaling the length expanded L by any single . These simulated master curves were invariant to the alteration of simulation parameters provided that the set of κ ij remained the same. This insight allowed us to develop a novel method of characterizing the experimental dynamics. Namely, we could experimentally determine L ij s , κ ij , and ϕ c , collapse the experimental data the same way as the simulations (i.e. L=L ij s , ϕ/ϕ c ), and compare the two to predict the dynamics of many competing alleles in a range expansion. As discussed below, this technique ultimately allowed for accurate predictions of the evolutionary dynamics of the four competing strains and, surprisingly, allowed us to make much more precise measurements of selective differences between strains.
As mentioned above in the Experimental Results section, using the procedures outlined in the Materials and Methods, we had previously determined R 0 = 3.50 ± 0.05 mm and D w = 0.100 ± 0.005 mm ( Table 2). In order to fit L ij we needed to measure v ij w . By tracking the growth of a more fit sector sweeping through a less fit strain (see the Materials and methods), we found that each strain swept through mCherry with a wall velocity of v iR w ¼ 0:06 AE 0:02 (as seen in Table 2); we could not detect the wall velocity of the eYFP and eCFP sweeping through the black strain.
In principle, the measured values of R 0 , D w , and v ij w should have allowed us to totally calibrate the three key parameter combinations. For example, followed from the measurement of L ij s using the known value of R 0 . Unfortunately, the final parameter L ij s ¼ D w =ðv ij w Þ 2 was more difficult to calibrate. Using v iR w ¼ 0:06 AE 0:02, we found that L iR s ¼ 30 AE 20 mm; the error on this value was too large for it to be predictive in our simulations. Furthermore, as we were unable to accurately measure the wall velocity of the eCFP/eYFP strains sweeping through the black strain, we could not calculate the corresponding selection length scale. We therefore needed a new technique to determine L ij s . As our eCFP and eYFP strains were neutral within error, we treated our system as composed of one neutral (N) eCFP/eYFP strain, a red (R) mCherry strain, and a black (B) strain (q = 3 colors). As the eCFP/eYFP expanded faster than the black followed by the mCherry strain, we needed to determine the values of L NR s , L NB s , and L BR s . To fit L ij s more precisely than that from our direct measurement of wall velocity, we competed pairwise combinations of strains in range expansions (i.e. the eCFP/eYFP strain vs. mCherry) and calculated the two-point correlation functions F ij (L, ϕ) at the maximum length expanded of L = 6.5 mm. As there were only two competing strains, there was only one L s . To fit the value of L s , we began by rescaling the experimental length expanded L by L s and ϕ by ϕ c ( Table 2) and calculated the resulting k ¼ ffiffiffiffiffiffiffiffiffiffiffi R 0 =L s p . Note that L s simultaneously rescales the length expanded L by L s and sets the value of k ¼ ffiffiffiffiffiffiffiffiffiffiffi R 0 =L s p , changing the shape of the collapsed correlation function. We then ran a simulation at the set value of κ (the chosen simulated values of L s and ϕ c did not matter due to the collapse) and then compared the collapsed experimental dynamics to our simulation. Fig 7 illustrates the fitting procedure by displaying the experimentally rescaled two-point correlation function F NR (the solid red line) at a length expanded of L = 6.5 mm between our eCFP/eYFP strain (N) and our mCherry strain (R)

(inoculated at fractions of 2/3 and 1/3 respectively) and simulated universal correlation functions corresponding to different values of L s (dashed lines).
To determine the best-fitting value of L s , we calculated the sum of the squared displacements weighted by the inverse of the experimental standard error squared between experiment and simulation. The best-fitting L s was determined by finding the value which minimized the weighted sum of squares. To estimate the error in our fit, we assigned each potential value of L s a probability proportional to the inverse of the weighted sum of squares, normalized the probability distribution, and set the error in our fit of L s to the confidence intervals of the probability distribution.
Our fit values of L ij s and κ ij using this technique are listed in Table 3; the values of κ ij are also plotted in Fig 5. Although this technique was about a factor of 5 more precise than using the measured wall velocities v iR w to determine L ij s , the upper bounds of the 95% confidence intervals were still very large as seen in Table 3; the potential values of L s had a very large tail.
To test that the resulting L s and κ could accurately predict the experimental dynamics at all L and not just the L where the correlation functions were fit, we plotted the experimental average fraction and correlation functions (solid lines, Fig 8) as we varied L and compared their values to those predicted by simulation (dashed lines, Fig 8). Fig 8 uses the same set of experimental data as that from Fig 7. The simulation using the fit parameters always closely tracked the experimental values at all L, suggesting that our fitting technique was robust and could be used to describe the dynamics of our strains.  Genetic drift and selection in many-allele range expansions Having determined L ij s and κ ij from pairwise competitions between strains, we tested if we could predict the average fraction, the two-point correlation functions, and the annihilation asymmetry when the four E. coli strains were grown together (treating the eYFP and eCFP strains as neutral, so q = 3) in an independent experiment. We inoculated the four strains in equal proportions. Fig 9 shows experimental measurements of the average fractions and twopoint correlation functions (solid lines) together with simulated predictions (dashed lines) that used the independently fit values; no additional fitting parameters were used. The predicted average fractions and correlation functions closely tracked the dynamics for L ≳ 3 mm. We   Genetic drift and selection in many-allele range expansions attribute the deviations for L ≲ 3 mm to image analysis artifacts resulting from the presence of the black strain (see the Image Analysis section in the Materials and methods). At the largest length expanded of L = 6.5 mm where artifacts were minimal, the experiments matched the predictions within error. All average correlation functions at this length expanded were successfully predicted by the simulations; we only display F NR for simplicity.
In addition to predicting the average fractions and correlation functions, the simulation with our fit L ij s and κ ij predicted that the annihilation asymmetry would deviate only slightly from neutrality (at most a change of 0.1) over the length expanded by our colonies in every experiment, agreeing with our findings (Fig 4). This can be readily observed in Fig 6 which displays a simulation of two neutral strains (our eCFP and eYFP strains) and a less fit strain (our mCherry strain) inoculated in equal fractions. If we rescale the maximum distance expanded by our colonies of L max = 6.5 mm by the smallest selection length (this results in the largest possible change of ΔP) of L NR s ¼ 13 mm, L max /L s * 0.5 and ΔP increases from 0 (neutrality) to at most 0.1. This small deviation from neutrality is within the uncertainty of our experimental measurement of ΔP. Evidently, certain quantities, like the average fraction and correlation functions, show signs of selection before others (in this case, the annihilation asymmetry).
The quantitative agreement between our model and our experiments suggests that the onedimensional annihilating-coalescing random walk model can indeed be used to predict the dynamics of many competing strains with different fitnesses in a range expansion.

Discussion
We investigated the evolutionary dynamics of four competing strains of E. coli in radial range expansions with differing selective advantages. We measured the average fraction F i of each strain, the two-point correlation functions F ij between strains, and the annihilation asymmetry ΔP with our image analysis toolkit [34]. Our simulations, which model the expansions as a one-dimensional line of random walkers subject to deterministic drift on an inflating ring, showed us that these three quantities could be collapsed onto universal curves for a fixed set when the length expanded by the colony L was rescaled by any of the lengths L ij s ¼ D w =ðv ij w Þ 2 and the angular distance between strains ϕ was rescaled by 0 c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8D w =R 0 p . To test if the random walk model could predict experimental dynamics, we independently calculated experimental values of L ij s , κ ij , and ϕ c and compared the dynamics between the two. The simulations accurately predicted the dynamics of the average fraction, correlation functions, and annihilation asymmetry when all four of our strains were present with no additional fitting parameters. The annihilation asymmetry ΔP is a quantity unique to range expansions with three or more strains and, to the best of our knowledge, has not been studied previously. Our results illustrate the importance of considering domain wall annihilation and coalescence when more than two strains are present and suggest that the annihilating-coalescing randomwalk model can act as a useful predictive tool when describing the evolutionary dynamics of range expansions composed of an arbitrary number of competing alleles with different fitnesses.
Along the way, we introduced a new technique that compared universal simulated correlation functions to experimental correlations to fit L ij s . The resulting values of L ij s were about a factor of 5 more precise than directly evaluating L ij s ¼ D w =ðv ij w Þ 2 with the wall velocities extracted from the growth of sectors. Given our fit L ij s , we evaluated v ij w using L ij s ¼ D w =ðv ij w Þ 2 and the known D w ; we compare these values to those extracted from single sectors in Table 4. The wall velocities from both measurements agreed within error, but the wall velocities obtained from our correlation method were at least a factor of two more precise than tracking single sectors. Although the correlation technique dramatically increased the precision in our evaluation of L ij s , the resulting precision increase for v ij w was less pronounced as v ij w / 1= ffiffiffiffi ffi L ij s p . Nevertheless, it is clear that the correlation technique can be used to precisely extract small differences in fitness between spatially competing strains.
Our work illustrates that the annihilating-coalescing random walk model can predict the experimental dynamics of an arbitrary number of competing alleles with different fitnesses in microbial range expansions. It is possible that this model could predict the dynamics of range expansions occurring outside of the laboratory, especially if the expanding organisms' underlying motion did not completely smear out the population's spatial structure; the organismal motion could potentially be accounted for by increasing the domain wall diffusion coefficient D w . To predict the dynamics of expansions, however, the annihilating-coalescing walk model relies on a key set of parameters: the set of L ij s , the set of κ ij , and ϕ c . We found that the set of L ij s could not be predicted from the independent radial expansion velocities of our strains; standard techniques [17] using the relative ratio of expansion velocities to predict v ij w , and thus L ij s , yielded inconsistent results (see S2 Appendix where we quantify the discrepancy and postulate why it occurred). As the set of L ij s is so fundamental to the evolutionary dynamics of range expansions, future work should investigate why relative radial expansion velocities could not be used to accurately predict v ij w and thus L ij s and whether this phenomenon is specific to E. coli range expansions or our specific strains. It would also be interesting to incorporate the reported super-diffusive motion of domain walls [8,9] into our simplified simulations and theoretical analysis. The random walk model's ability to successfully predict the evolutionary dynamics of our experiments suggests that annihilating and coalescing genetic domain walls subject to diffusion and selection-induced displacement provide a useful conceptual framework from which to understand range expansion dynamics.

Materials and methods Strains
We used four E. coli strains (labelled BW001, BW002, BW003, and BW012) with a DH5α background and plasmids whose sequences coded for spectrally distinguishable fluorescent proteins. The unique colors were obtained by using the plasmid vector pTrc99a [39] and the open reading frame for the respective fluorescent proteins. Strains BW001, BW002, and BW003 expressed eCFP (cyan/blue), Venus YFP (yellow), and mCherry (red) respectively, and were identical to the E. coli strains eWM282, eWM284, and eWM40 used in Ref. [40]. Note that these three strains were isogenic and differed only by the open reading frames corresponding to their respective fluorescent proteins. The final strain, BW012, was a mutated descendant of strain BW002 (yellow) that fluoresced at a decreased intensity, appearing black, while retaining its ampicillin resistance from the pTrc99a vector. Throughout this work, no additional mutations were introduced or observed. We therefore consider that these four strains Genetic drift and selection in many-allele range expansions correspond to four different alleles. Throughout the paper, we refer to the strains as eCFP, eYFP, mCherry, and black.

Experimental setup
To prepare saturated cultures, strains were inoculated in 10mL of 2xYT media and were shaken for approximately 16 hours at 37˚C. After vortexing each saturated culture and obtaining their concentration via optical density (OD-600) measurements, appropriate volumes (e.g., 1:1:1 mixtures of three strains) were added to an Eppendorf tube with a final volume of 1mL. The Eppendorf tube was then vortexed to uniformly mix the strains. A volume of 2 μL was taken from the vortexed tube and placed on center of a 100 mm diameter Petri dish containing 35 mL of lysogeny broth (LB), ampicillin at a concentration of 100 μg/mL, and 1.25% w/v bacto-agar. The carrier fluid in the resulting circular drop evaporated within 2-3 minutes, depositing a circular "homeland" of well-mixed bacteria onto the plate. After inoculation, plates were stored for 8 days upside down (to avoid condensation) in a Rubbermaid 7J77 box at 37˚C with a beaker filled with water; the water acted as a humidifier and prevented the plates from drying out. The plates were occasionally removed from the box and imaged (at roughly 24 hour intervals) using the brightfield channel to determine the radius of the colony as a function of time. On the eighth day, the plates were imaged in both fluorescent and brightfield channels. The number of replicate plates used are stated next to the respective experimental results. If we noticed that a mutation had occurred during an expansion (mutations usually presented themselves as unexpected large bulges at the front of a colony or as distortions in fluorescent intensity), we discounted the colony.

Image acquisition and analysis
We imaged our range expansions with a Zeiss SteREO Lumar.V12 stereoscope in four channels: eCFP, eYFP, mCherry (fluorescent channels), and brightfield. In order to analyze a colony with a maximum radius of approximately 10 mm using a single image, we stitched four images together with an overlap of 20% using AxioVision 4.8.2, the software accompanying the microscope. We blended the overlapping areas of the images to lessen the impact of background inhomogeneities. An example of a stitched image can be seen on the left side of Fig 10. Stitching introduced small artifacts such as vertical lines near the center of our expansions; we verified that these did not affect our results.  [34]. Images were acquired for four overlapping quadrants and stitched together to obtain a single image with a large field of view. Overlapping regions were blended to minimize inhomogeneities. To obtain the binary masks, pixels with fluorescence above background noise were marked as "on." A visual comparison of the raw data and the masks confirm that our binary masks accurately reflect the location and shape of individual sectors. https://doi.org/10.1371/journal.pcbi.1005866.g010 Genetic drift and selection in many-allele range expansions To extract the local fraction of each strain per pixel, we first created binary masks for each fluorescence channel indicating if the corresponding E. coli strain was present. We utilized the "Enhance Local Contrast" (CLAHE) algorithm [41] in Fiji [42], an open-source image analysis platform, to help correct for inhomogeneities in background illumination. After applying the CLAHE algorithm, a combination of automatic thresholding and manual tracing yielded a binary mask of each channel, an example of which is shown in Fig 10; the image on the left is an overlay of an experimental range expansion's fluorescent channels and the image on the right is the overlay of the corresponding binary masks. A small amount of manual tracing was required near the edges of our colonies because our fluorescent lamp provided uneven illumination; resulting dark regions could barely be identified above background noise. As we mainly used manual tracing near the edge of the colonies where the monoclonal sectors were well defined, we found that our procedure was very reproducible. To alleviate this problem, future work could utilize brighter strains or a more advanced imaging setup.
We mapped the binary images to the local fraction of each E. coli strain in the following way: if N binary masks (corresponding to N colors) were "on" at a pixel, the local fraction of their corresponding channels was assigned to be 1/N. Although this assignment produces inaccuracies (i.e., if one strain occupied 90% of a pixel and the other occupied 10%, our algorithm would register both as 50%), domain boundaries were the only areas besides the homeland and the early stages of the range expansions where multiple strains were colocalized. The black strain was defined to be present at pixels reached by the range expansion in which no other strains were present. Although this definition introduced errors at radii close to the homeland with significant color overlap, the error became negligible at large radii as quantified in Supplementary S5 Fig. Once we determined the fraction of each strain at each pixel, we were able to extract quantities such as the total fraction of each strain in the colony and spatial correlations between strains at a given expansion radius.
The mask in Fig 10 highlights that sector boundaries can be used to determine local strain abundance. Although it is possible to extract the position of every domain wall from each strains' local fraction, it is challenging to actually track a single wall due to collisions between walls. To address this problem, we created a binary mask of the edges in our images and labelled the edges of each domain. Annihilations and coalescences were counted manually within Fiji [42]; automated measures were not accurate enough.
It is worth pointing out that in this paper, we ignore the three-dimensional structure of our colonies and describe them by our two-dimensional images taken with the stereoscope. We justify this approximation because the initial diameter of our colonies is at least a factor of 10 larger than their height (less than 1 mm as judged by a ruler), so they are effectively twodimensional, and because the strain composition of our colonies does not vary with height inside the colony. We confirmed that strain composition does not vary with height by using a confocal microscope to probe the internal structure and also by taking a pipette tip, scratching it through a sector, growing the cells touched by the tip in overnight culture, and verifying that plated single colonies from the culture were the same color as the sector.

Measuring radial expansion velocities u i
We used the average expansion velocity of each strain for radii R > R 0 as a proxy for selective advantage, similar to previous work [17,35]. In three independent sets of experiments using different batches of agar plates (the main source of variability in our experiments), we measured the diameter of 12 expansions of each strain approximately every 24 hours following the protocol for range expansions with two or more strains. To account for biological variance, sets of four of the 12 colonies were created from independent single colonies; no statistical difference was seen between biological replicates. The diameters were determined by manually fitting a circle to a brightfield image of the expansion three times and averaging the measured diameters. Fig 11 shows the average radius increasing with time for each strain from one of our experiments. In every experiment, the eCFP and eYFP strains had the fastest expansion velocities (the respective datapoints overlap in Fig 11), followed by the black strain, and then finally the mCherry strain. The expansion velocity slowly decreased as a function of time; we attribute this to nutrient depletion in the plates.
The radial expansion velocity of each strain was obtained by using linear regression to fit the radius versus time for radii greater than R 0 . We calculated the average radial expansion velocity between the three sets of plates and reported its error as the standard error of the mean; see Table 1. Additionally, we quantified the dimensionless selective advantage of each strain relative to the slowest growing mCherry strain following [17] via s iR = u i /u R − 1 where the R indicates the mCherry strain (red) in each experiment. The selective advantages were consistent, within error, when we calculated the velocities u i and u R over different time intervals. We averaged s iR across our three experiments and reported its error as the standard error of the mean as seen in Table 1.
The eCFP and eYFP strains had an average selective advantage of 9%, similar to the experiments of Weber et al. [35] which found, despite the fact that they used different E. coli strains and plasmids, that the expression of mCherry decreased the expansion velocity of their strains by approximately 15% in certain "fast growth" environmental conditions. Our black strain had an approximately 6% enhancement over the mCherry strain. Differences in radial expansion velocities of this magnitude have been used to study yeast S. cerevisiae and E. coli range The error bars (comparable to symbol size at early times) are the standard errors of the mean calculated from 12 replicate expansions for each strain. The eYFP and eCFP strains had the fastest expansion velocities (data points overlap in the plot) followed by black and then mCherry. R 0 is the radius at which expansions with competing strains typically demix into one color locally; R 0 is approximately 1.75 times the initial inoculant radius of 2 mm (see Fig 1).
expansions in the past [9,17]. To investigate the source of this fitness defect, we took the plasmids from our original strains, inserted them into a different set of clonal DH5α cells, and inoculated the new eCFP, eYFP, and mCherry strains in equal proportions in a range expansion. We saw that the average mCherry fraction decreased by 10% at a radius expanded of R = 10 mm, matching the results of Fig 2, suggesting that the presence of the plasmids was responsible for the fitness defect.
From Table 1, it is clear that the variance in s iR was large between different sets of agar plates. Although s iR varied significantly, the order of expansion velocities between the strains was consistent; the eCFP and eYFP strains always expanded faster than the black strain which expanded faster than the red strain. Importantly and in stark contrast to s iR , the demixing radius R 0 , wall velocities v ij w , and diffusion coefficient D w were very consistent between sets of plates (measured below), resulting in consistent evolutionary dynamics between our competing strains.

Comparing well-mixed fitness to fitness from expansion velocities
To test if the radial expansion velocity differences were related to the basal growth rates of our strains in liquid-culture, we competed all of our strains against mCherry in 10mL flasks of 2xYT growth media. We created three independent replicates of each pairwise competition (9 tubes in total) in the flasks by inoculating 60% of mCherry and 40% of the other strain in midlog phase. We passaged saturated E. coli samples into new media every 24 hours; we determined the population composition of each flask using a BD LSR Fortessa FACS machine when the cells were passaged. We competed the strains for 72 hours, corresponding to approximately 45 generations (doubling times). We found that every strain grew faster than mCherry in every replicate in liquid culture as judged by a decrease in fraction of mCherry over time. Following previous work [43], we used the decrease in mCherry fraction to determine the dimensionless "well-mixed" selective advantage s wm iR ¼ g i =g R À 1 of the strain competing against it, where g i is the growth rate of strain i and g R is the growth rate of mCherry. We list the measured values of s wm iR in Table 1. The radial expansion velocity fitness s iR did not agree with the liquid-culture fitness s wm iR within error, in contrast to previous experiments with the yeast Saccharomyces cerevisiae [17]. However, every strain in liquid culture still grew faster than mCherry. Interestingly, in wellmixed culture, the black strain had the largest growth rate followed by eCFP and eYFP (they had the same growth rate) and then mCherry, disagreeing with the order of radial expansion velocities (where black expanded slower than eCFP and eYFP). As the growth rate differences were small, it is possible that additional factors allowed the eCFP and eYFP strains to expand faster than black on agar. Our E. coli switched from log to stationary phase in the 24 hour cycle; the changing environment may have resulted in a different order of fitnesses compared to the agar plates as well. Future work should investigate how the eCFP and eYFP strains expanded faster than the black strain despite a smaller basal growth rate and should also investigate how such small growth rate differences in liquid culture resulted in such large differences in radial expansion velocities on solid agar.
Measuring the local fixation radius R 0 When calibrating our model to experiment, the precise value of R 0 did not matter as long as each strain's local fraction could be accurately measured at that radius. Therefore, to maximize the length over which we could quantify range expansion growth, we defined the local fixation radius R 0 as the minimum radius where our image analysis package became accurate. For R < R 0 , our package predicted equal fractions of each strain due to the overlap of each channel in the homeland (see Fig 10). Therefore, to determine R 0 , we inoculated radial expansions with three strains in unequal proportions; we used 10% of two strains and 80% of another. The minimum radius where the fractions agreed with their inoculated values was R 0 = 3.50 ± 0.05 mm as seen in Supplementary S6 Fig. We found that this value of R 0 worked for all colonies.

Measuring the domain wall diffusion coefficient D w
Past work has found that E. coli colony domain walls fluctuate diffusively in certain conditions [18] and super-diffusively in others [8]. In our expansions, the domain walls appeared to fluctuate super-diffusively (as judged by tracking the position of domain walls and determining their variance vs. length expanded), but we were able to successfully fit the evolutionary dynamics using a diffusive theory. Creating a super-diffusive theory to describe the evolutionary dynamics of our system is beyond the scope of this paper. To obtain an effective diffusion constant D w and to test if the diffusive approximation adequately described our experimental dynamics, we fit the neutral Voter model's prediction of heterozygosity. The heterozygosity is the probability that two points separated by an angle of ϕ at a length expanded of L = R − R 0 are occupied by different strains and is thus a measure of spatial genetic diversity. The neutral Voter model's prediction of heterozygosity can be given in terms of the two-point correlations used in the main text or can be explicitly written as (see S1 Appendix) p is a characteristic angular correlation length (one of the key combinations of model parameters from the main text). For q colors inoculated in equal fractions, H 0 = 1 − 1/q.
We fit H(ϕ, L) to our experimentally measured heterozygosity of two neutral strains (eCFP and eYFP) on three independent sets of agar plates each with 14 range expansions. We averaged the heterozygosity at each L as can be seen in Fig 12 (error bars were omitted for readability; the same figure with error bars can be found as Supplementary S7 Fig). As we had previously measured R 0 = 3.50 ± 0.05 mm, and H 0 = 1/2 for two neutral strains inoculated at equal fractions, D w is the single free parameter in eq (11). We consequently fit D w at each L with non-linear least-squares, averaged the D w from the three independent experiments, and found D w = 0.100 ± 0.005 mm; the reported error is the standard error of the mean between the experiments. The value of the diffusion constant is on the same order of magnitude as that from previous work [18]. Fig 12 shows the Voter model's fit (dashed lines) together with the experimental heterozygosity (solid lines) for one set of plates using our values of D w and R 0 . The fit closely matches the experimental heterozygosity suggesting that a diffusive description of E. coli domain motion is justified. We use this value of D w for all strains. In principle, D w may depend on ij, the particular domain wall type. However, we checked that the measured value of D w did not vary for our all ij (all strain) combinations by examining the variance in domain wall position versus length expanded; the variances agreed within error and were thus consistent with a constant D w . The two-point correlation functions in the main-text were well fit by a constant D w as well. Unlike the Voter model and our simulations, the experimental heterozygosity at zero separation H(L, ϕ = 0) fails to vanish due to overlap between strains at domain boundaries; this effect is less pronounced at large radii because the effective angular width of boundaries decreased. The discrepancy between the theoretical and experimental heterozygosity is larger at small lengths expanded because the overlap between strains is larger; our image analysis is consequently less accurate.

Measuring the domain wall velocities v ij w
We used image analysis to directly quantify v ij w from the angular growth of more-fit sectors. Characteristic single sectors of each strain sweeping through the mCherry strain can be seen on the left side of Fig 13. In radial expansions, more fit strains should, on average, sweep  [34]. We tracked the angular growth of sectors sweeping through a less fit strain, ϕ − ϕ 0 , as a function of ln(R/R 0 ) to obtain v ij w . Right: 40 traces of each strain sweeping through mCherry from one set of agar plates. The translucent lines are the individual traces, the solid lines are the mean angular growth hϕ − ϕ 0 i, and the shaded area is the standard error of the mean. The slope of the mean angular growth is 2v ij w . https://doi.org/10.1371/journal.pcbi.1005866.g013 Genetic drift and selection in many-allele range expansions logarithmic spirals through less fit strains at large lengths expanded, as verified in yeast expansions [17]. It can be shown that the average angular width of a sector of strain i sweeping through strain j is given by (see S1 Appendix for more details) where ϕ is the angular width at radius R and ϕ 0 is the initial angular width of the domain at R 0 . 2v ij w can thus be extracted from the slope of a linear regression fit of hϕ − ϕ 0 i vs. ln(R/R 0 ) as seen on the right side of Fig 13. By tracking domain walls directly, we found that more fit strains (eCFP, eYFP, black) swept through the less fit mCherry strain with a wall velocity of v iR w ¼ 0:06 AE 0:02. We could not accurately measure the wall velocity of the eCFP and eYFP strains sweeping through the black strain. The wall velocity was significantly smaller than expected from the basal independent expansion velocities of our strains (Table 1); potential explanations for this phenomenon are discussed in S2 Appendix. The magnitude of the velocities were consistent between experiments (using 40 single sectors on three sets of plates) but were too imprecise to be predictive in our models.

Simulation methods
Lattice simulations of range expansions, especially radial ones, can suffer from artifacts arising from the preferred directions of the lattice. It is possible to use an amorphous Bennett model lattice [44] to mitigate some of these effects [32]. Instead, we developed a simple off-lattice method that treats the domain walls as annihilating and coalescing random walkers moving along the edge of an inflating ring. The basic idea of the simulation is illustrated in Fig 14. We incorporate both the random, diffusive motion of the domain walls as well as deterministic movement due to selection. The radial expansion procedure is most easily understood by first considering a linear range expansion simulation for which the simulation steps are as follows: 1. Create a line of N 0 microbes of width a at the linear frontier. Assign each microbe one of the q potential alleles.
2. Identify genetic domain walls by locating neighbors with different alleles; assign type ij to each wall where i and j are the strains to the left and right respectively. Assign a relative "growth rate" r ij to each wall characterizing the bias in the probability that strain i divides into new territory before strain j. Two such domain walls are shown in a radial expansion in Fig 14. 3. Choose a wall at random and move it a distance a (the width of the cells) to the left or right; this represents the competition to reproduce and occupy new space at the population frontier. We use periodic boundary conditions for the domain wall positions along the line, meaning that the domain walls live on a circle, as shown for the radial case in Fig 14. (a) Jump to the right with a probability of P r ¼ 1 2 1 þ r ij À Á or to the left with probability Note that domain walls separating neutral strains (r ij = 0) will jump to the left or right with equal probability and that r ij 1. Our simulation's diffusion coefficient per length expanded, characterizing the random motion of the domain walls, can be shown to be D w = a/2 when r ij is small while its wall velocity per length expanded, characterizing the deterministic displacement of domain walls due to selection, can be shown to be v ij w ¼ r ij 1. Our algorithm, thus far simulating only linear expansions, can easily be extended to simulate radial geometries. To incorporate the radially inflating perimeter, we note that a domain wall at a radius R will jump an angular distance of δϕ = a/R, as shown in Fig 14. As the radius of our experimental expansions increases approximately linearly with generation time, we describe its radius as R = R 0 + at. We thus account for inflation by using a time-varying angular Fig 14. A schematic of the simulation procedure for a radial expansion. The initial population is a circle of cells of radius R 0 = N 0 a/2π, where N 0 is the initial number of cells and a is a cell width. During each time step (generation), the expansion advances a distance a; the radius consequently grows according to R(t) = R 0 + at where t is the time in generations. The dashed circle shows the population after one generation time. Each domain wall position is tracked on the inflating ring (solid lines). At each time step, domain walls (two shown) hop to the left or right with probability P l and P r , respectively, with an angular jump length δϕ a/R(t), and the position is updated (dashed lines). After each domain wall movement, the time in generations is incremented by 1/N where N is the number of domain walls present. For a linear simulation, the radius is simply not inflated in time, i.e. R(t) = R 0 . jump length of If there are N 0 individuals at the frontier, R 0 is given by R 0 = N 0 a/(2π). This modification of the domain wall step size δϕ is the only difference between the radial and linear cases! In contrast to algorithms that follow the position and state of every organism at the front of a colony, our algorithm only tracks the positions of domain walls and is consequently much faster per generation as the sectors coarsen, allowing for simulations of larger colonies. All strains were inoculated in equal fractions except for the experiment with 10% of eCFP, 10% of eYFP and 80% of mCherry. The annihilation and coalescence rates (the slope of the respective curves) decrease as radius increases as there are less domain walls due to previous collisions and also because inflation decreases the probability of two walls colliding per length expanded. As the number of colors increases, coalescences occur more often than annihilations. (TIF) We ran four simulations where we varied v w (the velocity that the two other more fit strains swept through the less fit strain) and R 0 such that L s ¼ D w =v 2 w changed but k ¼ ffiffiffiffiffiffiffiffiffiffiffi R 0 =L s p was fixed. Each simulation has a different symbol in the plot. We found that F ij could be collapsed at the same L/L s as long as κ remained fixed (we arbitrarily set it to κ = 0.4) and as long as the angular variable ϕ was rescaled by 0 c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8D w =R 0 p . If ϕ c approached the system size ϕ c % 2π, F ij could not be collapsed onto the above curves due to finite size effects. Note that even though we only show F RR , all correlation functions F ij could be collapsed using this procedure. To estimate the image analysis artifacts introduced by using a non-fluorescent, black strain we performed an experiment with three fluorescent strains (eCFP, eYFP, and mCherry in equal initial proportions) and analyzed the data twice: once where we included all three fluorescent channels and once where we excluded the eCFP channel and treated it as if it were a black strain. We compared the black-substituted average fractions F i (the dashed lines) to the real fractions as a function of radius (the solid lines). At a small radius relative to R 0 = 3.5 mm, the error from introducing a black strain was large; this is likely because we defined black as the absence of any other channels and channels typically had large overlaps close to the homeland. At large radius, the error from introducing a black strain was negligible. (TIF) S6 Fig. Determining R 0 . To fit the radius R 0 where our image analysis package became accurate, we inoculated 80% of mCherry, 10% of eCFP, and 10% of eYFP in 10 range expansions and tabulated the average fraction of each strain. The inoculated fractions are illustrated by dashed lines. As seen in the plot, at a radius of approximately R 0 = 3.50 ± 0.05 mm the measured average fractions were closest to the inoculated fractions. Our image analysis package inaccurately predicted fractions in the homeland because of significant overlap between the strains. (TIF) We simulated a single fit sector sweeping through a less fit strain. It is expected that the fit strain sector dynamics satisfy hϕ − ϕ 0 i = 2v w ln(R/R 0 ) and Var(ϕ) = 4D w (1/R 0 − 1/R), as seen in S1 Appendix. To test that our simulation appropriately reproduced this behavior, we quantified the average angular growth hϕ − ϕ 0 i and angular variance Var(ϕ) as we varied the simulation parameters N 0 (initial number of cells), r (selective advantage of the fitter strain), and d (distance the colony expanded each generation). The cell width, a, was kept constant. These parameters relate to the sector dynamics via D w = a 2 /(2d), v ij w ¼ ar ij =d, and R 0 = (N 0 a)/(2π). We confirmed that both the average angular growth hϕ − ϕ 0 i and angular variance Var(ϕ) had the correct functional form and dependence on the microscopic parameters (the dashed black line). In the main text, we used d = a for simplicity. (TIF)