Range Expansion with Mutation and Selection: Dynamical Phase Transition in a Two-Species Eden Model

The colonization of unoccupied territory by invading species, known as range expansion, is a spatially heterogeneous non-equilibrium growth process. We introduce a two-species Eden growth model to analyze the interplay between uni-directional (irreversible) mutations and selection at the expanding front. While the evolutionary dynamics leads to coalescence of both wild-type and mutant clusters, the non-homogeneous advance of the colony results in a rough front. We show that roughening and domain dynamics are strongly coupled, resulting in qualitatively altered bulk and front properties. For beneficial mutations the front is quickly taken over by mutants and growth proceeds Eden-like. In contrast, if mutants grow slower than wild-types, there is an antagonism between selection pressure against mutants and growth by merging of mutant domains with an ensuing absorbing state phase transition to an all-mutant front. We find that surface roughening has a marked effect on the critical properties of the absorbing state phase transition. While reference models, which keep the expanding front flat, exhibit directed percolation critical behavior, the exponents of the two-species Eden model strongly deviate from it. In turn, the mutation-selection process induces an increased surface roughness with exponents distinct from that of the classical Eden model.


Introduction
Spatial structure strongly influences the dynamics of evolving biological systems and often gives rise to qualitatively different outcomes as compared to well-mixed populations [1,2]. During the last years microbiological experiments have progressively been used to shed light on the dynamics of spatially extended systems [3,4,5,6,7,8,9], and spurred theoretical investigation [10,11,12,13,14,15,16]. Since microbial colonies typically grow from some initial seed, a particularly interesting question to ask is how both a colony's morphology and its internal composition are shaped by the growth rates of the different strains it is composed of, and by the interactions between these strains. Already the simplest scenario, the spreading of two selectively neutral strains or species shows intriguing phenomena which makes the evolutionary outcome quite distinct from well-mixed populations [17,18,19,20,21,22]. Experimental investigations of expanding E. coli and S. cerevisiae colonies containing two fluorescently labeled but selectively neutral strains have shown that the population differentiates along the growing front and thereby segregates into well-defined domains [18,22]. This is caused by demographic fluctuations: Since mainly cells at the leading front of the growing colony access nutrients and reproduce, the effective population size is small and neutral dynamics leads to local fixation of strains and thereby generates sectoring of the population [22].
In general, however, microbial communities are heterogeneous and composed of multiple strains, which may have different growth rates or show other kinds of distinct phenotypic features, cf. e.g. [23,16,24,25,26,27,28]. One well-known phenomenon, which is the focus of this work, are bursts of new sectors of mutants during the growth of bacterial colonies [29,3,30,31,22]. In a growing bacterial colony, which initially consists of one phenotype (the wild-type) only, mutations of this strain may appear during reproduction of individuals. If these mutations happen at the leading front and are beneficial, i.e., if the mutant strain has a larger growth rate than the wildtype strain, mutant regions along the front not only advance faster but also expand laterally. Consequently mutant sectors take over ever larger parts of the front in a quasideterministic fashion [32,33,22]. While deleterious mutations are, for the same reason, handicapped by selection, they may still form large clusters along the front, if they appear frequently enough. In this case mutant sectors are no longer spatially separated, but may coalesce. At a critical mutation probability, the mutants' selective disadvantage is effectively balanced and mutants may take over the front [22]. If back-mutations are prohibited, the front remains trapped in an all-mutant state. In the language of nonequilibrium statistical mechanics, the critical mutation rate marks a phase transition between an active state, for which the front is composed of wild-types and mutants, and an absorbing homogeneous state, composed of mutants only [34,35,36,37].
What makes this absorbing state phase transition in a growing bacterial colony interesting is the intricate interplay between the morphology of the growing front and the evolutionary dynamics of the colony. Depending, among others, on the particular type of bacterial strain, nutrient concentration and softness of the agar surface, bacterial colonies exhibit a kaleidoscope of possible morphologies [38,3,30,4,6,39,40]. The Eden model [41] has been devised to describe the growth of bacterial cultures with a compact morphology on which we focus in this work. A hallmark of this model is the generation of rough fronts with characteristic features closely resembling recent experimental observations [42]. The roughness of the front directly affects the trajectory of interfaces between domains of different strains [32]: the undulations of the front are imprinted in the meandering of the domain boundaries on all length scales, as has recently also been observed experimentally [18,22]. This surface-induced meandering speeds up coalescence of clusters, i.e., the roughness of the propagating front strongly affects the temporal evolution of the population's composition. This suggests that front roughening is highly relevant for the nature of the phase transition from a heterogeneous to a homogeneous population at the expanding front.
It is precisely this issue which we would like to address in this manuscript. To this end we study bacterial range expansion using a two-species Eden model, which incorporates surface roughness, selection and irreversible mutations. We intend to gain deeper insight into the interplay of these key features of the dynamics and their relative importance for the transition to the absorbing state. While our findings are mainly of general importance for a broader class of multi-species growth models, we also expect that real world range expansions of bacterial colonies are subject to this coupling and carry its signature in the evolving patterns. In the remainder of the introduction we give a concise overview of previous work on surface roughness and absorbing state phase transitions as relevant for this work. A more in-depth discussion and comparison with the results of our work is given in the final section.
Both, discrete numbers and roughness of the expanding front, are intrinsic to surface growth models [43,44,45], which mimic the stochastic advance of particles into empty space. One particular growth model, the Eden model [41] (of which there are three, slightly different, variants [46]) has been devised to describe the growth of bacterial cultures. Mesoscopically its evolution is captured by the Kardar-Parisi-Zhang (KPZ) equation [47]. The KPZ equation constitutes a robust universality class which incorporates many surface growth models, like e.g. ballistic deposition and solid-onsolid models. There have been a number of generalizations to multi-species growth models [48,32,49,50,51,52,53,54,55,56,57,58,59,60], notably one of the Eden model which incorporates selection [32]. The coupled influence of mutations and selection on kinetic surface roughening, which is one of the topics of this work, has not been analyzed in detail so far. Even when neglecting roughness, multi-species propagation is not treated easily in more than one dimension. The reason for this is the intricate interplay of creation, annihilation and merging of clusters, which contain only one kind of individuals, at the leading front of the colony.
In the case of irreversible mutations both analytical results and numerical simulations for range expansion with flat fronts (neglecting surface roughness) predict a transition to an all-mutant absorbing state, which emerges even for deleterious mutations at a critical mutation probability [22]. For such models with flat fronts [22], the dynamics closely resembles that of a contact process [61]. As a consequence, it belongs to a broader class of absorbing state phase transitions whose main representative is directed percolation (DP) [62], a dynamic version of percolation [63]. The DP universality class of phase transitions to absorbing states has been found to display an enormous robustness with respect to alterations of the microscopic update rules and is considered a paradigm of non-equilibrium statistical physics [34,35,36,37]. How roughness of the front may influence phase transitions to absorbing states has, to the best of our knowledge, previously not been addressed.
The paper is organized as follows. In Section 2 we introduce a two-species growth process on a 2d lattice to analyze the general properties of range expansions of asexually reproducing microorganisms. It explicitly includes irreversible mutations from wild-types to mutants and selection between the two strains. The evolving system's morphology intimately depends on the antagonistic effects of new mutant domains being created at the front and others loosing contact to it. For abundant, deleterious mutations a phase transition to an absorbing state exists, which changes both the evolutionary dynamics and surface roughening behavior of the system qualitatively. We discuss properties of both surface and bulk morphology and map out the phase diagram in Section 3. The system's critical behavior near the transition is affected by the front's roughness, since the temporal evolution of the system is restricted to the growing front. This alters the critical properties of the absorbing state phase transition. In Section 4.1 we determine its critical exponents, which are different from those of the DP class (which is most often found for flat systems [34,35,36,37]). In addition, the different birth rates of the two strains induce an enhanced width of the front near the phase transition. Close to the transition the roughness exponent of our model is severely enhanced compared to that of the Eden model [43,44,45]. In addition to KPZ behavior we identify and characterize a critical roughening regime in Section 4.2. We conclude with a discussion of our results and a comparison to related models which study the coupling between surface roughening and domain dynamics.

Eden Model with Mutations
Range expansion into hitherto unoccupied territory proceeds in a non-homogeneous manner on the length scale of individuals. Along the leading front local protrusions emerge randomly and subsequently expand, thereby creating a rough front and an overall forward movement. The main features of this growth process are well captured by the classical Eden model, which was developed to mimic the growth of microbial colonies [41,44,45]. While some multi-species extensions of the Eden model have been analyzed [48,32,49,50], surface growth experiments with competing microorganisms have only been performed recently [18,22]. They reveal intriguing, non-trivial patterns if the population is comprised of distinguishable sub-populations. Mutations and selection can alter the growth dynamics by giving some individuals a growth advantage or by introducing qualitatively different organisms.
In this work we employ a lattice gas model to analyze the influence of mutations and selection on range expansion at rough, fluctuating fronts. We model microbial range expansion with mutations as a cellular automaton on a 2d semi-infinite square lattice of extensions L × ∞ with periodic boundary conditions in the transverse direction (see Fig. 1). At a given time t, each site (i, j), with i ∈ {1, . . . L} and j ∈ N, is either empty Eden model with mutations (color online). An initially line-shaped bacterial colony of length L, consisting of only wild-types, shown in dark gray (red), grows into an empty half-space. Wild-types reproduce at rate 1, given that they have a vacant nearest-neighboring lattice site. Individual offsprings are either wild-types or mutants (light gray) with probabilities 1 − p and p, respectively. Mutants reproduce like wild-types but at birth rate b. Back-mutations are prohibited. In the L-direction (transverse) periodic boundary conditions apply.
or occupied with an individual, which in turn can be either a wild-type or a mutant. We identify empty sites, wild-types, and mutants with state variables s = 0, s = 1, and s = −1, respectively. The state of the system at a given time is, therefore, specified by the set of occupation numbers s i,j ∈ {−1, 0, 1}. The system evolves as individuals reproduce: empty sites with s = 0 can change their state if an individual on a nearestneighboring site reproduces. We assume that individuals do not die and hence any site with s = ±1 remains in its state indefinitely.
To implement random-sequential update of a configuration we apply a simplified version of Gillespie's stochastic algorithm [64]. Only individuals at the front, defined as the set of occupied sites with at least one empty nearest-neighbor site, can reproduce. Of these, an individual is randomly, but proportional to its birth rate, chosen to reproduce: Mutants reproduce with relative birth rate b ≥ 0, while the birth rate of wild-types is 1 and thereby sets the timescale. Let N wt and N mut denote the number of wild-types and mutants with free neighbors, respectively. The overall birth rate of the population (at which production events happen) is given by b tot := (N wt + bN mut ). To account for different birth rates, each wild-type individual with an empty neighbor is chosen to reproduce with probability 1/b tot , while each mutant individual with an empty neighbor is chosen with probability b/b tot . The new individual is placed on a random empty nearest-neighbor site of the individual which has been chosen to reproduce. During wildtype reproduction, a mutation may happen with probability p. Thus, if the reproducing individual is wild-type (s = 1), the offspring is a wild-type with probability 1 − p and a mutant with probability p. If the reproducing individual is a mutant (s = −1), the offspring is necessarily also a mutant. Note, that by these rules back-mutations and multiple mutations are prohibited, and that all mutants are identical and reproduce with the same birth rate b. Assuming exponentially distributed reproduction times, the expectation time until the next reproduction event is b −1 tot . Hence, we update time by t → t+b −1 tot whenever an individual reproduces. Since, after some initial transient period, the average front position moves at constant velocity, one may take the longitudinal coordinate j as a proxy for time t.
The model, as described above, is a generalization of version C of the Eden model as introduced by Jullien and Botet [46]. It is the biologically most realistic version, as it focusses on occupied sites, i.e., individuals, rather than on empty sites (version A) or bonds between adjacent occupied and empty sites (version B). In the limit of vanishing mutation rate our model reduces to the model of Saito and Müller-Krumbhaar [32]. While we consider the case of a homogeneous initial front with uni-directional mutations, they analyzed the temporal evolution of an initially heterogeneous front in the absence of mutations, see Section 5 for more details. If not stated otherwise, we use a line of wild-type particles as initial condition, i.e., s i,j = δ j,1 . Since diffusion is not included in the model, surface configurations become frozen in the bulk, as observed for patterns in range expansion experiments [18].

Phenomenology
We now turn to a phenomenological description of the morphology of the evolving colony, as obtained from stochastic simulations; a representative realization of the range expansion process is shown in Fig. 2a. Starting out from an initial line of wild-types the growth front moves forward and, as a result of the stochastic individual birth processes, some parts of the front expand more rapidly than others. This leads to front roughening which first appears on length scales of the lattice constant, but as time progresses, the typical size of protrusions and indentations of the front grows both longitudinally and laterally.
As the range expansion proceeds, mutation events occur where a wild-type individual gives birth to a mutant, whose direct descendants create a new mutant sector growing between two wild-type domains. Mutant and wild-type domains are separated by domain boundaries. Since the reproduction rates of mutant and wild-type individuals In the bulk, the presence of mutant clusters is discernible as a density of wild-type sites ρ(j) < 1 at longitudinal position j. Mutant clusters in the bulk are characterized by their typical longitudinal and transverse extensions, which are identified with the correlation lengths ξ and ξ ⊥ , respectively. The surface roughness of the population front is characterized by the width w, defined as the standard deviation of the front position. For this realization on a lattice of length L = 128, we used birth rate b = 0.9 and mutation probability p = 0.016. (b) Averaged observables. By averaging over many independent realizations of the growth process fluctuation effects are suppressed and mean observables can be defined. Shown here is the mean density ρ(j) as function of longitudinal position j, averaged over 500 realizations for the parameters used in (a). For large j, the density of wild-type sites settles to a stationary value ρ s , as creation and annihilation of mutant clusters (and boundaries) at the front equilibrate.
differ in general, the boundary between their respective growth sectors performs a biased random motion. While for beneficial mutations, with b > 1, sectors consisting of mutants broaden on average, they have a tendency to decrease in size for deleterious mutations where b < 1. Fluctuations in the trajectory of the boundary arise mainly for two reasons: On short length scales, they are due to the intrinsic stochasticity of the birth events. On larger scales, roughening of the front drives a super-diffusive meandering of the sector boundaries [48,32,18]: As the population locally always expands normal to the front, the roughness of the front is imposed on the trajectories of the sector boundaries.
While the domain boundaries move transversely through the system, they may encounter other boundaries, resulting in mutual annihilation and an ensuing merging of domains. There are two distinct types of coalescence events of domain boundaries, cf. Fig. 2a: Either boundaries of different mutant clusters meet such that they merge and form a larger mutant cluster, or two boundaries of the same cluster meet, in which case a mutant cluster loses contact to the growing front and is trapped by wild-types. Note that boundaries are always created and annihilated pairwise, and regions of wild-types and mutants alternate at the front.
The relative frequency of events creating and annihilating sector boundaries determines the ultimate fate of the expanding front. In wild-type dominated regions of the front, mutation events enhance phenotypic heterogeneity and create new boundaries. At the same time, merging of mutant clusters creates homogeneous mutant regions. If the selective advantage of wild-types is too small to trap mutant clusters, coalescence events promote growth of mutant clusters, leading to more uniform front populations. Since we do not allow for back-mutations, the expanding front may end up in an absorbing state where it is completely taken over by mutants. For a finite system, this will eventually always happen even for deleterious mutations. The main question to ask then is how the corresponding fixation time scales with the system size L.
Inspecting Fig. 2 one can discern several morphological features of the expanding population, which we discuss phenomenologically now, and analyze quantitatively in the following sections. We discriminate properties of the front and of the bulk of the population. An important bulk observable is the density of wild-type sites, at longitudinal position j. As indicated in Fig. 2b, the ensemble averaged density decays, for L → ∞, towards some stationary value ρ s , which serves as an order parameter. A value of ρ s = 0 corresponds to the absorbing state where all individuals at the front are mutants. Any finite value indicates phenotypic heterogeneity with both mutants and wild-types present at the front of the expanding population. Mutant clusters are, in general, anisotropic and one has to distinguish between their extension parallel and perpendicular to the preferred direction of the range expansion: The corresponding longitudinal and transverse correlation lengths are denoted by ξ and ξ ⊥ , respectively. The front of the population is characterized by its average speed and its roughness. A good measure for the latter is the width w, defined as the standard deviation of the front's position,  Here the wild-type is lost in the bulk and the absorbing state, with only mutants at the front, is reached very fast, as mutants merge into extended clusters. As opposed to this, the active phase (cf. panels E, I, and J) is characterized by long lasting heterogeneous fronts with both mutants and wild-types present. Here, the average spatial distance between distinct mutation events is larger than the extension of the created mutant clusters, which repeatedly lose contact to the front. A non-equilibrium phase transition separates the two phases, where low relative birth rates b and high mutation probabilities p balance and rich self-affine patterns evolve (cf. panels A, F and K). Going from panel J to panel B at constant birth rate b = 0.5, the system changes from a heterogeneous front of wild-types and mutants to a homogeneous all-mutant front. The transition is characterized by diverging length scales (mutant clusters), vanishing wild-type density, and enhanced width of the population front, as discussed in the main text.
from its average value, Here h(i, t) is the local position of the front, defined as the largest j for which s i,j = 0.

Phase Behavior: Active and Inactive Phase
The morphology of the expanding population, as discussed above, depends on the values of the mutation probability p and the relative reproduction rate b of the mutant individuals, cf. Fig. 3. One can clearly identify two distinct phases: In what we call the active phase the population front is composed of both wild-types and mutants in a heterogeneous mixture. Here mutants are continuously created by mutation events and the ensuing mutant sectors are subsequently lost again by a coalescence process where sector boundaries meet and the domain of mutants looses contact to the front. This typically happens in a parameter regime where mutations are deleterious and rare. In contrast, if mutations become more frequent and/or their reproduction rate becomes larger, mutant clusters have a significant probability to merge and/or to sweep through the system, and thereby completely take over the population front. This is termed the inactive phase since the absorbing state, with only mutants at the front, cannot be left anymore.
For all beneficial mutations, where b > 1, we are well in the inactive phase, independent of the probability p at which mutations appear. Here, a sector created from a single mutation has a finite opening angle and hence grows laterally on average, as observed in experiments [3,22]. Therefore, already a single mutant can sweep through the population and take over the population front, as the boundaries of the respective growth sector have, on average, a transverse velocity directed outwards [32,33,22]. This fixation process is further accelerated by the merging of neighboring mutant sectors, cf. Fig. 3

panel L.
For deleterious mutations, where b < 1, there is an antagonism between merging of different mutant clusters, creating large mutant domains, and closing of individual mutant clusters, leading to an enlarged fraction of wild-types at the front. Since mutations go from wild-type to mutant only, phenotypic heterogeneity is maintained as long as there are some wild-type individuals left at the front.
An isolated sector of mutants can survive only for a finite time interval due to stochastic effects. These enable mutants to "surf" population waves in expanding populations [65,19,66,67]: If mutations appear at the front of an expanding population they have a twofold advantage compared to mutations in spatially homogeneous settings. First, the front can be seen as a perpetual population bottleneck, where demographic fluctuations are enhanced, which in turn reduces the effect of selection. Second, offspring of mutations at the front can spread into unoccupied territory, where there is less competition. By this "founder effect" mutants form clusters at the front that can reach much higher frequencies and evade extinction much longer than would be expected in spatially homogeneous settings. However, in the long run an isolated sector tends to loose contact to the front as it is outgrown by the wild-type as a result of the lower birth rate of mutant individuals, cf. panels E, I, and J of Fig. 3. The phenomenology changes qualitatively when mutation events become more frequent. Then, nearby mutant sectors can merge and thereby counteract the loss of mutant sectors by the coalescence of sector boundaries [22], cf. panels B, C, and G, of Fig. 3. The founder effect promotes growth of individual mutations, but persistence of mutants is guaranteed by merging of domains. As a consequence of this antagonism one expects a phase boundary p c (b) between the passive and the active phase. Hallatschek et al. considered mutations appearing at a flat front and found that for deleterious mutations there exists a critical mutation rate where mutants take over. We here consolidate this interesting finding and incorporate the roughness of the front. We in detail analyze the properties of the transition from a heterogeneous to an all mutant front in Section 4.

Phase Diagram
Since our model explicitly excludes back-mutations, the absorbing state, for which the wild-type strain has lost contact with the population front and only mutant individuals are present, cannot be left. For finite systems, L < ∞, this absorbing state is eventually always reached, as even for deleterious mutations mutant clusters of arbitrary size can appear through rare fluctuations. As noted above, the key quantity to analyze is the average time for this to happen as a function of mutation probability p, relative birth rate b and system size L. This mean fixation time t f (p, b, L) is defined as the mean time t when the number of wild-type individuals with empty neighbors, N wt , becomes zero for the first time.
The mean fixation time generally diverges with growing system size, L → ∞. The results from our stochastic simulations, shown in Fig. 4, allow to distinguish three generic cases: (i) For large mutation probability p or beneficial mutations b ≥ 1 we find t f ∼ ln L. This is the same result as for well-mixed populations [68,69]. We take this asymptotic law as a hallmark for the inactive phase. Note that for isolated beneficial mutations (neglecting merging of mutant clusters) one finds that the extinction time scales linearly in the system size [22]. Hence, the logarithmic law arises from the merging events. (ii) In the active phase we have t f ∼ exp(cL), with some constant c. (iii) The phase boundary p c (b) between the active and the inactive phase is characterized by power law behavior of the mean extinction time: t f ∼ L z with a dynamical exponent z. This signature is well known from previous studies of phase transitions to absorbing states [37,36] (see Section 4). From our numerical data we estimate z = 1.05 ± 0.05.
Since the phase transition from the active phase to the inactive phase is accompanied by a qualitative change in the L-dependence of the mean fixation time t f from logarithmic to exponential behavior, we may use it to map out the phase diagram, see Fig. 5.
We find a phase transition line p c (b), which ends at the points (b, p) = (0, p * ) and (b, p) = (1, 0). For b = 0 mutants do not reproduce at all and each mutant is created by a reproducing wild-type individual in an independent event with probability p. For low p the clusters have a typical size of one site and are obstacles around which the wildtypes grow. For larger p more extended clusters are formed by mutation events that happen on neighboring sites. For the wild-types this corresponds to an Eden process (b) Same data as in (a), but rescaled with system size L and distance to the phase transition ∆, according to finite size scaling as detailed in Section 4.1. All data collapses onto a master curve as can be inferred from the semi-logarithmic plot. The inset is a double logarithmic plot of the same data, which depicts the characteristic scaling above and below criticality, t f = L z τ ± (L/∆ −ν ⊥ ). As a guide to the eye we included best fits (black lines) to these functions given by τ + (x) = 3.8 exp(10.5x)/(1 + x) and τ − (x) = 0.26 ln(1 + 10.5x)/x. on a site percolation network, where wild-type individuals are identified with occupied sites (present with probability 1 − p) and mutant individuals with empty sites (present with probability p). Thus (b, p) = (0, p * ) is a multi-critical point of the transition line, below which the wild-type strain, for L → ∞, keeps growing on the infinite percolating cluster. This implies that p * = 1 − p isp c ≈ 0.407, where p isp c is the percolation threshold of isotropic site percolation on a 2d square lattice [63].
For non-vanishing birth rates, b > 0, a newly formed mutant cluster can grow to some significant size. As a result clusters originating from more distant mutation events can merge before loosing contact to the front. The phase transition line marks the set of mutation rates p c (b), where, for a given birth rate of mutants b, the typical spatial distance between two mutation events becomes comparable to the typical extension of individual mutant clusters. Thus, as we move to larger and larger birth rates, the critical mutation probability p c declines.
For b → 1 mutations become neutral, i.e., selection pressure vanishes and the random motion of sector boundaries is not biased anymore. Individual clusters can grow unlimited by fluctuations, and consequently the smallest number of mutations suffices for the population to be in the inactive phase. If p = 0 no mutants are created and we recover the classical one-species Eden model composed of growing wild-types [41]. Thus, at (b, p) = (1, 0) the front is composed of only wild-type individuals and the absorbing state is never reached. Our numerical data suggests that p c approaches this point in a cusp-like singularity.
For b > 1 mutations are beneficial, which means mutations and selection are not antagonists anymore, and hence a phase transition is absent. Here, the mutants take over the system quasi-deterministically for all finite values of p.
To rationalize our findings for the phase diagram we consider a phenomenological mean field theory for the dynamics of the density of the wild-type strain at the front of the expanding population, n = N wt /(N wt + N mut ). On average, wild-type individuals are lost by mutation at a rate p, and gained or lost through natural selection depending on the relative growth rates of the wild-type and mutant individuals. The effect of natural selection is proportional to n(1 − n), which vanishes if the front is composed of one strain only. This leads to the following rate equatioṅ where s(b) is an effective selection strength. The functional form of s(b) can be determined phenomenologically: The rate equation (3) has the fixed points n 1 = 0 and n 2 = 1 − p/s(b). In the active phase, n 1 = 0 is unstable, while n 2 is stable. At the phase transition the two fixed points merge and interchange their stability in a transcritical bifurcation. Solving n 1 = n 2 gives the mean field transition line p mf c (b) = s(b). For b → 0 the problem reduces to isotropic percolation and hence s(0) = p * . Certainly the selection coefficient s(b) has to vanish if mutants reproduce as fast as wild-types, s(1) = 0, which gives us the other end of the transition line. Close to it the transition line can be approximated by a power law, p mf c (b → 1) ∼ (1 − b) µ , reflecting the cusp singularity we found in the simulations. The simplest ansatz for the mean field transition line, which fulfills these constraints is A fit to our numerical derived phase transition line gives µ = 1.41 ± 0.03. The phenomenological mean-field transition line, depicted in Fig. 5, is in very good agreement with our numerical results.

Critical Behavior
In this section, we in-depth investigate the properties of the absorbing state phase transition from the active into the inactive phase. Without loss of generality we focus on a fixed value for the birth rate of mutants, b = 0.5, for which we found the critical mutation probability to be p c (b = 0.5) = 0.159 ± 0.001. The qualitative behavior of all observables stays the same along the transition line. Close to the special points b = 0 and b = 1 crossover effects become more pronounced, which we do not examine here.

Bulk Properties
We first address bulk properties of the system, that is observables measured sufficiently far away from the rough growth front. Since we are dealing with a phase transition to an absorbing state, we expect to find four independent critical exponents [34,35,36,37]. A common choice of observables is: The stationary density ρ s := ρ(j → ∞) of active sites (wild-types); the survival probability P s (Survival in this context means the wildtype having contact to the growth front.) of a single active site (wild-type) in a front of inactive sites (mutants); the longitudinal correlation length ξ , and the transverse correlation length ξ ⊥ . For L → ∞, all these observables diverge like power laws in the control parameter ∆ := p c − p in the vicinity of the phase transition. The respective critical exponents are defined through The above observables are understood as ensemble averages. The stationary density ρ s and the survival probability P s are 0 in the inactive phase, ∆ < 0, since in this case any heterogeneous composition of the front is unstable. For finite systems, L < ∞, both the stationary density ρ s and the survival probability P s are identical to 0 for the inactive and the active phase. Though it may be extremely rare, it is always possible -even in the active phase -that large enough mutant clusters appear and finally lead to a front consisting of mutants only, which is the absorbing state. Hence, ρ s and P s are not particularly useful observables for L < ∞, and one instead has to examine the time-dependent density ρ(∆, t, L) and survival probability P (∆, t, L), and also time-dependent correlation lengths ξ (∆, t, L) and ξ ⊥ (∆, t, L).
We measured the critical exponents, defined in Eqs. (5a-5d), using different methods. Starting with a line of wild-types in the active phase we let the system evolve, disregarding realizations where the absorbing state was reached. This allows us to use stationary observables, instead of more involved dynamical ones. From each realization we extracted for each mutant cluster i) the longitudinal extension , ii) the transverse extension ⊥ and, iii) the number of cluster sites or mass m. In analogy to percolation theory [63], the correlations lengths are then calculated by where k runs over all observed clusters. The computed correlation lengths depend on both the distance to the phase transition ∆ and the system size L. However, close to the critical point all macroscopic observables are invariant under scaling transformations of the form [34,35,36,37] ∆ → c∆ , x → c −ν x , (7b) where c is some positive rescaling factor. This implies for the correlation lengths the following scaling forms, where the f # (x) are universal scaling functions with the asymptotic behavior Indeed, after rescaling, our simulation data for both correlation lengths ξ # (∆, L) collapse onto master curves for all L and ∆, if we take see Fig. 6a. Note that this means that close to the transition the longitudinal and transverse correlation length show the same scaling behavior, ξ ⊥ ∼ ξ , at least within the accuracy of our simulations. Simply put, a fluctuation of twice the size takes twice as long to decay, which explains that close to the phase transition fixation times are proportional to the system size, t f ∼ L, which we found in Section 3.3.
With the critical exponents ν ⊥ and ν at hand, rescaling of the fixation time t f (L, ∆) is easily achieved as well. Note that as a time-like observable t f scales like a longitudinal distance. Again employing the phenomenological scaling theory for absorbing state phase transitions, Eqs. (7a-7e), we obtain where z = ν /ν ⊥ is the dynamical critical exponent. The overall scaling function can be split up into two distinct parts, τ + and τ − , which exhibit characteristic logarithmic and exponential behavior for ∆ > 0 and ∆ < 0, respectively. Using this scaling, all data for the fixation times collapse, as shown in Fig. 4b.
In the inactive phase, the density ρ(∆, j, L) decays exponentially with j and the absorbing state, with ρ = 0, is reached fast. The decay length is proportional to ξ which can easily be found by fitting an exponential function to the measured density profile. Applying finite size scaling, we again find data collapse for critical exponents ν ⊥ = ν = 1.6 ± 0.1, see Fig. 6b.
Further evidence for these values of the critical exponents is obtained by calculating the (active state) correlation functions in the bulk. Averages are with respect to all lattice indices i, j and independent realizations. The correlation functions are expected to behave as where σ # = 2β/ν # . To obtain the correlation lengths we fitted this expression to our data with the fitting parameters ξ # . In plots of the correlation lengths, we again find good collapse of data for ν ⊥ = ν = 1.6 ± 0.2 (not shown).  Figure 6.
Rescaled correlation lengths from simulations started with a line of wild-types as initial condition for different lattice sizes. (a) Rescaled correlation lengths from cluster-size distributions in the active phase. Longitudinal correlation length ξ (upper plot) and transverse correlation length ξ ⊥ (lower plot) have been computed from mutant cluster masses m and cluster extensions and ⊥ (see text). Correlation lengths depend on both L and distance to the phase transition ∆. Taking ν ⊥ = ν = 1.6 ± 0.1 all data collapse to master curves under rescaling. (b) Rescaled longitudinal correlation length as found from decay of wild-type density in the inactive phase. The wild-type density ρ decays exponentially, with a decay length proportional to the longitudinal correlation length. Again all data collapse to a master curve for ν ⊥ = ν = 1.6 ± 0.1.
Next, we consider the stationary density ρ s (∆), which equals the fraction of wildtype sites for large j, that is sufficiently far away from the initial line for the density to relax to its stationary value. As shown in Fig. 7a, a double-logarithmic plot asymptotically gives power law behavior with the exponent Deviations from the asymptotic power law are well described by the finite-size scaling form where g(x) is a universal scaling function, cf. Fig. 7b. Two characteristic regimes can be distinguished corresponding to the scaling law ρ s ∼ ∆ β , for a given L and sufficiently far from criticality, and the scaling law, ρ s ∼ L −β/ν ⊥ , for a given distance ∆ from the critical point and system sizes smaller than the transverse correlation length, L ξ ⊥ . To measure β , the exponent associated with the survival probability when starting from a single seed, P s , different initial conditions must be used. Instead of a line composed of only wild-types, a single wild-type is placed in a line of mutants. A typical realization for these initial conditions is shown in Fig. 8a. From an ensemble  of realizations the survival probability P s (j) and the mean number of wild-types N (j), both as functions of distance j to the initial line, can be extracted. After transients have decayed and sufficiently far away from the front (where there might exist empty lattice sites), we find power laws for both quantities when close to the phase transition. For P s the exponent is given by −β /ν while for N it is denoted by θ. Fig. 8b shows results for P s (j), from where we obtain β /ν = −0.32 ± 0.2. Using our result for ν we find Similarly, from the simulation data for N (j) we obtain θ is related to the other critical exponents by the generalized hyperscaling relation [70] in d = 1 dimension Using our previous results for the critical exponents we find good agreement with the above result, Eq. (18), within the error margin. The critical exponents for the phase transition to an absorbing state are summarized in Table 1. For systems with infinitely many absorbing states, like ours, it is known that critical exponents, especially those related to single seed initial conditions (in our case β and θ), can vary significantly with the configuration away from the seed [71,72]. The correct exponents, which satisfy the hyperscaling relation of Eq. (19), are obtained if the configuration away from the seed is a typical configuration of the absorbing state. In our case this amounts to a rough front. To check for this dependence of the critical exponents, we initialized single wild-type seeds at arbitrary front positions of all-mutant Eden models grown to saturation. Data extraction works similar, but due to roughness and the non-uniqueness of the growth direction the initial seed is not necessarily the one with the smallest j-value, which complicates the evaluation slightly. These simulations take more time since independent saturated fronts have to be generated for every realization, which makes it harder to obtain data with high precision. However, we find β /ν = −0.33 ± 0.03 and θ = 0.31 ± 0.04 (data not shown), which indicates that our system is robust with respect to the initial configuration. We are therefore confident, that the critical exponents of the phase transition, see Table 1, are indeed correct up to the given precision.
To check if the DP universality class is recovered if roughness is neglected, we also considered a simplified model with synchronous update rules. Here, an individual at site (i, j) has two possible parents, at sites (i − 1, j − 1) and (i, j − 1). Depending on their states, we define the probabilities P (s i,j |s i−1,j−1 , s i,j−1 ) for the type of offspring s i,j , incorporating mutations and selection: If both parents are mutants, the offspring is inevitable mutant as well; if both parents are wild-types the offspring is wild-type if no mutation happens; if one parent is wild-type while the other is a mutant, the mutant reproduces with rate b/(1 + b) and the wild-type otherwise. In the latter case a mutation can happen, resulting in a mutant offspring. All events are summarized by the probabilities and P (−1|., .) = 1 − P (1|., .). It is easily seen that this corresponds to the Domany-Kinzel cellular automaton [73], with probabilities p 1 and p 2 . Indeed, for b = 0.5 we find the critical mutation probability for the flat model p flat c ≈ 0.077, corresponding to p 1 ≈ 0.62 and p 2 ≈ 0.93. This lies well on the transition line of the Domany-Kinzel automaton [74]. As one would expect, the critical exponents of the flat model are that of the DP universality class (data not shown). These conclusions are further consolidated by earlier studies for a two-species ballistic deposition model with kinetics that allows surface growth only at sites where one type of particles has an exposed position on top of the incidence or neighboring column [60]. These particular growth rules imply that the surface roughness does not affect the bulk dynamics such that it can be mapped onto a one-dimensional contact process [61] which is in the DP universality class [34,35,36,37]. In conclusion, the interplay between surface roughening and domain dynamics in our two-species Eden model is clearly responsible for the deviation of its critical behavior from that of the DP universality class.

Front Properties
We now turn to characterizing the front of the expanding population, i.e. to surface roughness properties of the model. As argued in Section 3.1, the presence of mutations changes the roughening behavior of the growth front qualitatively: Close to the phase transition large mutant clusters form extended parts of the leading front and since mutants and wild-types reproduce with different rates, the front's roughness is strongly enhanced. In the case of detrimental mutations, b < 1, mutant-dominated regions trail behind compared to the average front; see for example Fig. 3, panel F.
The width is the key observable in the analysis of kinetic surface roughening processes. These processes can be organized into universality classes, which are characterized by symmetry properties and the dimensionality of the process under consideration [43,44,45]. Each growth model's universality class has a unique set of two exponents: The growth exponent γ and roughness exponent α, defined through with the crossover time between these two regimes scaling as t × ∼ Lz, wherez = α/γ is the dynamic roughening exponent.  Figure 9. Width of the growth front as a function of time in the active phase (color online). The temporal evolution of the width w(t) of the growth front depends on both the system size L and the distance to the phase transition ∆. w initially grows Edenlike, γ = 1/3 (bold dashed line). Close to the phase transition (i.e. at small ∆) we find a crossover to a growth exponent γ = 0.73 ± 0.03 (bold straight line), as a result of different growth rates of mutants and wild-types. For intermediate times two scenarios can be distinguished, see main text for details: (i) Either the surface width saturates at a finite length-dependent value as a result of finite system size, or (ii) mutant clusters at the front saturate at a size smaller than L. In the latter case Eden growth is recovered on large length and time scales. As the width saturates we recover the Eden roughness exponent α = 1/2 but with a p-dependent amplitude of the saturation width. In the former case, enhanced roughening proceeds until saturation. The saturated width then scales like L α , with a different roughness exponent α = 0.91 ± 0.01. All data shown are for b = 0.5, but results hold for different birth rates b as well.
In the limit of vanishing mutation rate, p = 0, our model reduces to the Eden model, which falls into the well-known KPZ universality class [47], with exponents γ = 1/3 and α = 1/2. For p > 0, the time evolution of the width depends on whether the system is in the active phase or in the inactive phase. In the inactive phase we find a transiently enhanced width as compared to the Eden model and a subsequent return to it. We attribute this to mutants outcompeting wild-types rapidly, as is typical for the inactive phase, followed by Eden growth of the resulting mutant front (data not shown). As shortly discussed in Section 3, the active phase is characterized by a heterogeneous front. In finite systems, this is the case until a mutant cluster covers the whole system size L, which then takes the system into the absorbing state with subsequent return to Eden-like growth. Therefore, the roughness characteristics of the heterogeneous front in the active phase can only be analyzed if the absorbing state has not yet been entered. Hence, for the following analysis, realizations which reach the absorbing state are stopped, and simulation data are used only up to this point in time.
We now discuss the different regimes and scenarios of surface roughening in the active phase. Fig. 9 displays our numerical results for a range of mutation rates and two different systems sizes. Since we start with an initial condition where the lattice only contains wild-types, the initial surface roughening is that of a single species, i.e. Eden growth: w ∼ t γ with γ = 1/3. As time passes, mutants are created, leading to growing mutant clusters. Their extension in growth direction evolves like the correlation length ξ (∆, t, L), which initially grows proportional to t, as follows from scaling. Since these clusters grow slower than wild-type clusters this leads to additional surface roughening. As for small t both contributions to roughening are independent of ∆ and L, the width gets dominated by differential expansion velocities at a crossover time t ×,1 = O(1).
As a result of the ensuing strong coupling between domain dynamics and surface roughness, the interface width then grows more rapidly than in the Eden model. We observe a regime with an altered growth exponent γ = 0.73 ± 0.03, which becomes more extended upon approaching the phase transition to the inactive phase, ∆ → 0. There are two distinct scenarios by which this enhanced roughening regime may end: (i) Either the surface width saturates at a finite length-dependent value as a result of finite system size, or (ii) mutant clusters at the front saturate and their typical size reaches the asymptotic extension ξ ∞ ⊥ ∼ ∆ −ν ⊥ ≤ L. § The first case requires that one is sufficiently close to the critical point such that the cluster size ξ ⊥ does not saturate at ∆ −ν ⊥ ≤ L, but mutant clusters continue to grow during the whole roughening process. Then, our simulations show that there is a direct crossover from enhanced surface roughening, w ∼ t γ , to saturation, w ∼ L α with a roughness exponent α = 0.91 ± 0.01. The corresponding scaling form close to criticality reads w c (L, t) = L γ ŵ c (t/Lz ) (22) with the dynamic critical exponent for surface rougheningz = α /γ = 1.25 ± 0.05, cf. Fig. 10a. (Note that deviations from the critical scaling behavior are all due to initial transient Eden growth: all systems were started with an initial condition of wildtype individuals only.) A roughening exponent α close to unity indicates a very jagged surface. For models with large roughness exponents, especially for super-roughening with exponents larger than 1 (see, e.g., [75,76]), it is known that they may exhibit anomalous scaling [77,78], in the sense that the roughness exponent as determined from w(L, t → ∞) may only be an effective exponent. The actual exponent can be obtained from the structure factor or power spectrum [79,80,45] S(k, t) := h(k, t)h(−k, t) , (23) § A third scenario where the typical lateral extension of mutant clusters reaches the extension of the system ξ ⊥ ∼ L is also possible but excluded by our sampling method.
where h(k, t) is the Fourier transform of the height fluctuations. From Eq. (21) it follows that the structure factor of the Eden model scales as with the scaling function The rescaled structure factor of our model is shown as the inset of Fig. 10a. For the trivial case p = 0 we observe data collapse as expected for the Eden model. At criticality ∆ = 0 (p = 0.159), the scaling function for the power spectrum changes with time due to the growth dynamics of mutant clusters, i.e. the time dependence of the cluster size ξ ⊥ (t). For asymptotically large times, when a single mutant cluster spans the system, a broad regime emerges with a roughness exponent α = 0.92 ± 0.03 corroborating our results obtained from analyzing the interface width. This clearly shows that, for ∆ → 0, the Eden model with mutations and selection is no longer in the KPZ universality class, but exhibits different asymptotic roughening behavior. We attribute this behavior to the strong coupling between critical domain growth and surface roughening.
Moving further away from the critical point, typical domains may saturate at their asymptotic value ∆ −ν ⊥ ≤ L within a crossover time t ×,2 ∼ ∆ −ν . Then, on length scales larger than ξ ∞ ⊥ , roughening proceeds, but slower than before, and one recovers Eden-like roughening with w ∼ t γ . This regime is most pronounced if L ξ ∞ ⊥ , and finally saturates in a constant surface width with w ∼ L α . The crossover between critical roughening and return to Eden roughening is determined by a time scale inherent to the system which only depends on the distance to the critical point but not on the system size. We may, therefore, rescale time by the crossover time, t ×,2 = ∆ −ν , and rescale width w by the corresponding value w(t ×,2 ) = t γ ×,2 = ∆ −γ ν , cf. Fig. 10b. All curves collapse to a crossover scaling function, from which they only deviate due to finite size effects.

Conclusion and Outlook
A generalization of the Eden model to a two-species model, with uni-directional mutations and selection due to different reproduction rates, has interesting properties from the viewpoints of both dynamical phase transitions as well as kinetic surface roughening. We find that surface roughness has a marked effect on the critical properties of the absorbing state phase transition. While reference models [22], which keep the expanding front flat, exhibit DP critical behavior, the exponents of our generalized Eden model strongly deviate from it. In turn, the mutation-selection process in the bacterial colony induces a increased surface roughness with exponents distinct from that of the KPZ universality class. for a series of system sizes as indicated in the graph. Neglecting initial transients, all data collapses for α = 0.91 ± 0.01 andz = α /γ = 1.25 ± 0.05. Both growth exponent γ and roughening exponent α are markedly different from their KPZ counterparts. The inset shows the rescaled structure factor S(k, t) for system size L = 512. For the Eden model all data collapses, and from the black solid line, proportional to k −2 , we recover α = 1/2 from Eq. (25). For the critical system, ∆ = 0, we find an extended scaling regime with S ∼ k −2.85±0.05 (black dashed line) corresponding to α = 0.92 ± 0.03. (b) Crossover scaling. Neglecting initial transients, data can be rescaled with respect to the intrinsic crossover time t ×,2 ∼ ∆ −ν , where the mutant cluster extension ξ ⊥ saturates at ∆ −ν ⊥ . If the system size L is much smaller than ξ ∞ ⊥ , the front width grows to saturation with the new growth exponent γ = 0.73 ± 0.03 (black solid line), whereas for ξ ∞ ⊥ L the saturation width is reached after return to Eden-like roughening with γ = 1/3 (black dashed line). Saturation, due to finite size, happens where curves deviate from the common envelope and become horizontal.
For our model the critical exponents of the longitudinal and transverse correlation length, ν and ν ⊥ , are identical within the error margin. As the aspect ratio ξ /ξ ⊥ can still be different from unity, this does not make our system isotropic. However, it indicates that there is no longer a rigorous distinction between a time-like longitudinal and a space-like transverse direction as for DP. To some degree, this might have been expected from the microscopic update rules of the Eden model, which have no preferential direction of growth. Indeed, for vanishing birth rate of mutants, the system can be mapped to isotropic percolation, and it looses any preferential growth direction. One may also argue that front roughening in concert with locally normal growth of the interface leads to "mixing" of transversal and longitudinal directions, and thereby promotes the same scaling behavior of both directions. It remains to be seen how this result can be obtained from a renormalization group calculation.
The critical exponents summarized in Table 1 are significantly different from those of the DP universality class, which are observed for a broad range of systems exhibiting phase transitions to absorbing states [34,35,36,37]. The stability of the directed percolation universality class is reflected in the DP hypothesis, a conjecture formulated by Janssen [81] and Grassberger [82]. It asserts that models characterized by a positive one-component order parameter, which exhibit a phase transition to an unique absorbing state, belong to the DP universality class, provided that interactions are short-ranged and no additional symmetries or special kinds of disorder or fluctuations are present [34,35,36,37]. Strictly speaking, our model is not in contradiction with the DP hypothesis as it does not fulfill the condition of having a unique absorbing state. Instead, there are infinitely many configurations with only mutants at the front. However, even for models with infinitely many absorbing states one often finds DPclass universal behavior [71,72]. What distinguishes our model is that growth persists in the fluctuating absorbing state of an all mutant system. Most importantly, there is interplay between surface roughness and bulk properties: Even though we analyzed bulk properties of the system, the dynamics are restricted to the front and, therefore, convey imprints of the rough surface. This might also be a reason why it is difficult to find actual experimental systems that obey DP-like behavior: Any non-flat system may carry signatures of surface roughening in the numerical values of the critical exponents of the corresponding non-equilibrium phase transition.
The influence of surface roughness on domain boundaries, without the additional complication of mutations, selection and of merging clusters, is a rather complex problem on its own. The morphology of the boundary between two distinct but equally fast growing Eden clusters has previously been analyzed by Derrida and Dickman [48]. Numerically they found that the dynamics of the domain boundary depends on the local curvature of the front. It may perform a subdiffusive or superdiffusive random walk, and even move ballistically. For our model this implies enhanced complexity, because boundaries are preferentially created at protrusions, where wild-types are more prone to be, while merging of domain boundaries usually occurs at indentations of the front. Coarsening of domains in a two-species Eden model is in the focus of work by Saito and Müller-Krumbhaar [32], who also considered differential reproduction rates of the two species. As mutations are absent in their model, the faster growing species outgrows the slower one exponentially fast. In the neutral case they show that domain boundaries exhibit super-diffusive scaling and thereby explain the domain coarsening kinetics of their simulations. It would be interesting to find a suitable generalization of their analysis to include the case of irreversible mutations.
There are also growth models where surface roughness and particle configuration at the front mutually affect each other [49,50,56,57,58,59]. A theoretically wellstudied system is vapor deposition of binary films containing two different kind of molecules which have a tendency to phase separate [58,59]. Similar to the dynamics of the two-species Eden model, there is an interplay between surface roughening and phase ordering kinetics: the dynamics of the domain boundaries and the surface are coupled through the growth kinetics. One main difference between the models is that for binary films the non-equilibrium roughening process is coupled to the coarsening dynamics of a thermodynamic model (Ising model) exhibiting detailed balance, while in our model it is coupled to the far-from-equilibrium dynamics of an absorbing state phase transition. In the binary film particles of one type attach to domains containing particles of mainly the other type. If this is interpreted as "mutation", then mutations are bi-directional and symmetric in the binary film. In addition, the growth rules differ in some important aspects. While in our two-species Eden model, the two species have different growth rates, in binary film growth, particle attachment at domain boundaries and within domains differ [49,50,56,57,58,59]. A common finding of these studies on binary films is that phase ordering kinetics increases surface roughness on length scales comparable to domain size; this is phenomenologically similar to our findings but the critical exponents are quite different. The reverse coupling of surface roughness to phase ordering kinetics is subtle and depends on the details of the kinetic growth rules. If growth at domain boundaries is faster than within domains, the surface roughness imprints its scaling properties on the domain boundaries [59], as in our case. Else, it appears that the domain boundaries perform random walks or even show non-universal behavior [48,58,59].
The enhanced roughness, which is induced by differential front velocities, bears some similarities with pinning models [45]. There, inhomogeneities (e.g. obstructions) of the medium locally reduce the growth speed, similar to mutant clusters in our model. The presence of these heterogeneities changes the scaling of the interface, since it introduces a length scale ξ ⊥ , up to which one finds enhanced roughening. It is found that in the pinned phase, many pinning models can be mapped to DP, and the roughness of the surface is dominated by the DP critical exponents [83,84]. In contrast to pinning models, where the disorder is quenched, the mutant clusters in the Eden model with mutation are dynamically generated and strongly correlated with the growth dynamics. We regard the latter as the main reason why the ensuing roughness and growth exponents are not related to DP.
Summarizing, we here argue that for multi-species range expansion, surface roughness and domain dynamics interfere with each other which qualitatively changes both bulk and front properties. If absorbing front-states exist, phase transitions to absorbing states of a new universality class are possible. More research is needed to explore the coupling between surface roughening and evolutionary dynamics of range-expansion scenarios. It would be especially interesting to look at systems where the reproduction rate of an individual depends on the local composition of the front as might be the case in growing biofilms [24]. Such systems, where growth depends on a population's composition, have been analyzed in the context of one dimensional wave propagation [85], the dilemma of cooperation in spatial settings [86], games with cyclic dominance [87], and structured populations [88]. A discrete 2d growth model, which incorporates these effects, would be most interesting to analyze with the methods developed in this paper. Moreover, the universality class of our model should be tested for different lattices, and critical exponents should be determined for higher dimensional setups.