Bacterial range expansions on a growing front: Roughness, Fixation, and Directed Percolation

TL;DR

This paper models how bacterial populations expanding on a surface exhibit rough growth fronts and fixation, using coupled stochastic equations and renormalization group analysis to understand their critical behavior.

Contribution

It introduces a novel coupled stochastic PDE framework linking surface roughness with directed percolation in bacterial range expansions.

Findings

Roughness influences the critical behavior of bacterial fixation.

Scaling laws at the DP transition are derived through renormalization group analysis.

The model predicts specific growth patterns and genealogical structures.

Abstract

Directed Percolation (DP) is a classic model for nonequilibrium phase transitions into a single absorbing state (fixation). It has been extensively studied by analytical and numerical techniques in diverse contexts. Recently, DP has appeared as a generic model for the evolutionary/ecological dynamics of competing bacterial populations. Range expansion -- the stochastic reproduction of bacteria competing for space to be occupied by their progeny -- leads to a fluctuating and rough growth front, which is known from experiment and simulation to affect the underlying critical behavior of the DP transition. In this work, we employ symmetry arguments to construct a pair of non-linear stochastic partial differential equations describing the co-evolution of surface roughness with the composition field of DP. Macroscopic manifestations (phenomenology) of these equations on growth patterns and genealogical tracks of range expansion are discussed; followed by a renormalization group analysis of possible scaling behaviors at the DP transition.

Figures

Tables

Table 1. Table 1: Comparison of critical exponents for coupled DP and KPZ equations obtained from RG analysis in d = 4 − ϵ 𝑑 4 italic-ϵ d=4-\epsilon with uncoupled DP in d = 4 − ϵ 𝑑 4 italic-ϵ d=4-\epsilon and d = 1 𝑑 1 d=1   Henkel  et al. . We also include the numerically determined exponents for the two-species Eden model in d = 1 𝑑 1 d=1 from Ref.  Kuhr  et al. ( 2011 ) ; noting that the comparison to the trends from ϵ italic-ϵ \epsilon -expansion is not particularly favorable.

Exponent	$z_f$	${\nu }_{⟂}$	$\beta$	${\beta }^{\prime }$	$z_h$	$\chi$
DP/KPZ ($d=4-ϵ$)	$2+\frac{ϵ}{6}$	$\frac{1}{2}-\frac{5ϵ}{24}$	$1-\frac{ϵ}{3}$	$1-\frac{4ϵ}{3}$	$2$	$-ϵ/6$
DP ($d=4-ϵ$)	$2-\frac{ϵ}{12}$	$\frac{1}{2}+\frac{ϵ}{16}$	$1-\frac{ϵ}{6}$	$1-\frac{ϵ}{6}$
DP ($d=1$)	$1.58$	$1.10$	$0.28$	$0.28$
Two-Species	$1.05$	$1.6$	$0.50$	$0.51$	$1.25$	$0.91$
Eden Model ($d=1$)

Equations

$${\partial }_{t}f=K{\nabla }^{2}f+{\tau }^{\prime }f(1-f)+\sqrt{2D_{f}f(1-f)}\eta .$$

$${\partial }_{t}f=K{\nabla }^{2}f+\tau f-\frac{g}{2}{f}^2+\sqrt{2D_{f}f}\eta .$$

$${\partial }_{t}h=v_0+\nu {\nabla }^{2}h+\frac{\lambda }{2}(\nabla h{)}^2+\cdots +\sqrt{2D_h}\xi .$$

$${\partial }_{t}h$$

$${\partial }_{t}f$$

$$+{\beta }_1\nabla f\cdot \nabla h+{\beta }_{2}f{\nabla }^{2}h+\frac{c}{2}f(\nabla h{)}^2.$$

$${\partial }_{t}h$$

$${\partial }_{t}f$$

$$+{\beta }_1\nabla f\cdot \nabla h+{\beta }_{2}f{\nabla }^{2}h+\frac{c}{2}f(\nabla h{)}^2.$$

$${\tau }_{\mathrm{eff}}=\tau +{\beta }_2{\nabla }^{2}h+\frac{c}{2}(\nabla h{)}^2.$$

$${\partial }_t{N}_i=\nu {\nabla }^2{N}_i+{\alpha }_i{N}_i+\sum _j{M}_{ij}{N}_j+\sqrt{2D{N}_i}{\eta }_i,$$

$${\partial }_{t}N=\nu {\nabla }^{2}N+\overline{\alpha }N+\sqrt{2DN}\eta ,$$

$${\partial }_t{f}_i$$

$$+\sum _j{M}_{ij}{f}_j+2D{f}_i(1-f_i)/N+\sqrt{2D/N}{\eta }_i^{\prime },$$

$${\partial }_{t}h=\frac{2\nu }{\lambda }\overline{\alpha }-\frac{2\nu D}{\lambda N}+\nu {\nabla }^{2}h+\frac{\lambda }{2}(\nabla h{)}^2+\frac{2\nu }{\lambda }\sqrt{\frac{2D}{N}}\eta .$$

$${\partial }_{t}h=\nu {\nabla }^{2}h+\frac{\lambda }{2}(\nabla h{)}^2+\frac{2\nu }{\lambda }\overline{\alpha },$$

$${\partial }_t{f}_i=\nu {\nabla }^2{f}_i+\lambda \nabla f_i\cdot \nabla h+({\alpha }_i-\overline{\alpha })f_i+\sum _j{M}_{ij}{f}_j;$$

$$N_i(\mathbf{x},t)=\int \frac{d^d{\mathbf{x}}^{\prime }}{(4\pi \nu t{)}^{d/2}}$$

$$\times$$

$$h(\mathbf{x},t)=\frac{2\nu }{\lambda }\ln N\approx \underset{{\mathbf{x}}^{\prime }}{\max }\left[h({\mathbf{x}}^{\prime },0)-\frac{|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^2}{2\lambda t}\right],$$

$$\frac{1}{(4\pi \nu t{)}^{d/2}}\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^2}{4\nu t}\right){\left[e^{Gt}\right]}_{ij}⟹$$

$${\int }_{({\mathbf{x}}^{\prime },0)}^{(\mathbf{x},t)}\mathcal{D}[\mathbf{x}(t^{\prime })]\exp \left[-\frac{1}{4\nu }{\int }_0^{t}d{t}^{\prime }{\left(\frac{d\mathbf{x}(t^{\prime })}{d{t}^{\prime }}\right)}^2\right]$$

$$\times \prod _{m=1}^n{e}^{{\alpha }_{i(m)}(t_m-t_{m-1})}\prod _{m=0}^n{M}_{i(m),i(m+1)}.$$

$$h(x,t)$$

$$⟨[h(\mathbf{x},t)-h({\mathbf{x}}^{\prime },t^{\prime }){]}^2⟩=|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^{2\chi }g\left(\frac{|t-t^{\prime }|}{|\mathbf{x}-{\mathbf{x}}^{\prime }{|}^{z_h}}\right),$$

$$\frac{D\mathbf{v}}{Dt}\equiv {\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla v=\nu {\nabla }^2\mathbf{v}-\lambda \sqrt{2D_h}\nabla \xi ,$$

$$h^{\prime }=h+ϵ\cdot \mathbf{x},{\mathbf{x}}^{\prime }=\mathbf{x}-\lambda ϵt,t^{\prime }=t,$$

$${\xi }_{\perp }\sim \delta {\tau }^{-{\nu }_{\perp }},{\xi }_{\parallel }\sim \delta {\tau }^{-{\nu }_{\parallel }},\rho \sim \delta {\tau }^{\beta },P_{\mathrm{perc}}\sim \delta {\tau }^{{\beta }^{\prime }}.$$

$$\mathcal{P}[f,h]={⟨\delta (f-f[\eta ])\delta (h-h[\xi ])⟩}_{\xi ,\eta },$$

$$\mathcal{P}[\xi ,\eta ]\propto \exp \left\{-\int d^d\mathbf{x}dt\left[\frac{\xi (\mathbf{x},t{)}^2}{2}+\frac{\eta (\mathbf{x},t{)}^2}{2}\right]\right\}.$$

$$\mathcal{P}[f,h]\propto \int \mathcal{D}[\tilde{f}]\mathcal{D}[\tilde{h}]\exp \left\{-\mathcal{J}[\tilde{f},f,\tilde{h},h]\right\}.$$

$$\begin{split}{\mathcal{J}}_{0}=\int d^{d}{\bf x}dt\ \Big{[}&{\tilde{f}}\left(\partial_{t}f-K\nabla^{2}f-\tau f\right)\\ &+{\tilde{h}}(\partial_{t}h-\nu\nabla^{2}h-\alpha f)-D_{h}{\tilde{h}}^{2}\Big{]},\end{split}$$

$$\begin{split}{\mathcal{J}}_{\rm int}=\int d^{d}{\bf x}dt\ \bigg{[}&\frac{g}{2}{\tilde{f}}f^{2}-D_{f}f{\tilde{f}}^{2}-\beta_{1}{\tilde{f}}\nabla f\cdot\nabla h\\ &-\beta_{2}{\tilde{f}}f\nabla^{2}h-\frac{c}{2}{\tilde{f}}f(\nabla h)^{2}-\frac{\lambda}{2}{\tilde{h}}(\nabla h)^{2}\bigg{]}.\end{split}$$

$${\partial }_{t}h$$

Full text

I Introduction

Directed percolation (DP) serves as the paradigmatic example of a nonequilibrium phase transition between an active state (of coexisting components) and an absorbing state (with one component) Henkel et al. ; Ódor (2004); Täuber (2014). Quantitative understanding of the critical behavior at this transition is rooted in numerous theoretical Janssen and Täuber (2005); Täuber et al. (2005) and numerical Hinrichsen (2000); Henkel et al. analyses. As depicted in Fig. 1, DP can be visualized as the irreversible and directed growth of active sites on a lattice, with new sites added to the growing cluster one layer at a time. However, it is its interpretation as a single-species reaction-diffusion process that offers its most comprehensive formulation and has lead to the suggestion of its widespread prominance in the form of the DP hypothesis: the critical properties near an active-to-absorbing phase transition should generically be governed by the DP scaling exponents, provided the stochastic process is Markovian, the order parameter is decoupled from any other slow variable, there is no quenched disorder and there are no special symmetries Henkel et al. ; Täuber et al. (2005); Täuber (2014).

Recently, a model of competition between bacteria in an expanding population has been proposed as a new manifestation of the DP universality class Korolev et al. (2010); Korolev and Nelson (2011). Such bacterial range expansions, however, differ from traditional DP in one important aspect: the growing population front can be rough, an effect we highlight in Fig. 1 with sample realizations of a toy model. It is now widely appreciated that such a rough frontier can have significant influence on the evolutionary/ecological dynamics Hallatschek et al. (2007); Hallatschek and Nelson (2009); Korolev et al. (2012); Farrell et al. (2017). Experimental evidence reinforced by simple scaling arguments further suggests that the influence of rough growth leads to inter-species domain boundary fluctuations not consistent with the DP universality class Hallatschek et al. (2007). This analysis has been corroborated by a handful of numerical simulations of various modifications of the two-species Eden model allowing for cooperation, competition and/or mutations Saito and Müller-Krumbhaar (1995); Kuhr et al. (2011); Lavrentovich and Nelson (2014); Kuhr and Stark (2015). All have confirmed that the scaling exponents of the absorbing phase transition differ from those at the DP transition. In addition, when the two species have different expansion rates, the measured scaling exponents for the frontier fluctuations appear distinct from those predicted by the Kardar-Parisi-Zhang (KPZ) equation – a classic minimal model of surface growth Kardar et al. (1986).

In this article, we complement these studies with a field-theoretic analysis based on symmetry considerations. The resulting phenomenological model is composed of a pair of stochastic partial differential equations that describe the coupled time evolutions of the growth front and the DP composition field. Our template for this approach are previous studies in Refs. Drossel and Kardar (2000, 2003), where the effects of surface roughening on phase ordering were modeled by coupling an Ising order parameter to the KPZ equation for surface growth, as well as earlier studies of coupled KPZ equations Ertas and Kardar (1993). The construction of the model is detailed in Sec. II, where we start by recounting field theoretic formulations of the DP transition, and KPZ growth, and then indicate how different types of symmetry allow and constrain couplings between these fields. The phenomenology of the minimal set of couplings is discussed in Sec. III. In particular, we elucidate the macroscopic manifestations of the couplings on the shape of the growth fronts, and on the dependence of various velocities and transition points on the overall slope of the front. Section V discusses hidden symmetries of the model, and its mapping to other problems of interest. In particular, it is shown that in certain limits the deterministic version of the equations is exactly solvable (following linearization by a Cole-Hopf transformation), highlighting the emergence of fixation starting from a rough initial condition. The renormalization group analysis, carried out in Sec. VI, then reveals that the coupling between the DP and KPZ fields is relevant, modifying the scaling exponents of both the DP transition and the KPZ equation. Previous renormalization group calculations in this direction considered the effects of surface fluctuations with correlated spatial noise on the absorbing phase transition, but ignored feedback from the DP order parameter onto the surface Sarkar and Basu (2012), which we explicitly include.

II Coupling composition and roughness at a front

The DP active-to-absorbing phase transition can be analyzed in the continuum limit by considering the variations of a field $f({\bf x},t)$ that classically represents the density of active particles at time $t$ and position ${\bf x}$ in a $d$-dimensional space. It is then the disappearance of active particles due to death or coalescence that signals the transition to the absorbing state, $f=0$. Here, our motivation comes from range expansions of immotile bacteria Hallatschek et al. (2007), where we instead consider the competition between a pair of bacterial species: a wild type and a mutant. In this case, the order parameter $f$ characterizes the fraction of wild-type bacteria, with the transition to the absorbing state mediated by the irreversible takeover of the mutants. From this point of view, $f$ reflects the composition of the growing bacterial colony.

In this section, we introduce various limiting dynamics for this bacterial composition field and its coupling to surface fluctuations, working up stepwise from simpler to more complex, terminating with the full coupled roughness/DP equations. The goal is to pinpoint the new phenomenology introduced with each new level of complexity. To aid in the discussion, it is helpful to have in mind a representative microscopic dynamics, like the ones pictured in Figs. 1$(a)$ and $(b)$. In the evolutionary biology community, the underlying “stepping stone” model was introduced by Kimura and Weiss Kimura and Weiss (1964) and has been used extensively since then to examine competitive growth Korolev et al. (2010, 2010); Lavrentovich et al. (2013). The model is equivalent to the Domany-Kinzel (DK) cellular automaton rules for flat growth Korolev et al. (2010); Lavrentovich et al. (2013), and the rough variant is based on the “brick wall” solid-on-solid model used previously for the study of phase transitions in the growth of binary films Drossel and Kardar (2000, 2003), or to the “single step” growth model introduced by Meakin *et al. *Meakin et al. (1986) and Plischke *et al. *Plischke et al. (1987). We will discuss their specifics in more detail in Sec. III, but, briefly, different species – wild-type (gray bricks) and mutant (blue bricks) – grow out on a tilted square lattice by replicating into adjacent empty lattice sites with rates that can depend on their local neighborhood Korolev et al. (2010). The growth can be either restricted to the addition of one full new layer at a time [Fig. 1$(a)$], or growth can evolve organically allowing for a roughening of the surface [Fig. 1$(b)$].

II.1 Flat growth without mutation

The simplest variant of the model starts with a mixture of two species that grow (possibly at different rates) one layer at a time. In the absence of mutations, the distinction between wild-type and mutant disappears, and either species can potentially take over the population (i.e. there are two absorbing states). Figure 1$(a)$ (in the absence of mutation) suggests a mean-field dynamics (ignoring fluctuations) where the fraction of the first species $f$ grows in proportion to the difference in growth rates $\tau^{\prime}$ as $\partial_{t}f=\tau^{\prime}f(1-f)$. Allowing for demographic noise and migration (diffusion) leads to a Langevin equation, describing the evolution of the field $f({\bf x},t)$ as Korolev et al. (2010); Janssen (2005)

[TABLE]

Here, $K$ is the spatial diffusion (or migration) rate, and $\eta({\bf x},t)$ is delta-correlated (in space and time) noise, with intensity (variance) $D_{f}$.

Equation (1) describes a critical, albeit discontinuous, phase transition as $\tau^{\prime}$ goes through zero Henkel et al. . Indeed, for $\tau^{\prime}<0$ the mean-field steady-state value of the order parameter is $f^{*}=0$; whereas, for $\tau^{\prime}>0$, we have $f^{*}=1$. Thus, the order parameter experiences a discontinuous jump at $\tau^{\prime}=0$, from one absorbing phase ($f^{*}=0$) to another ($f^{*}=1$). The choice of absorbing phase at $\tau^{\prime}=0$ is a symmetry-breaking consequence of noise and initial conditions. This indistinguishability of the two species is reflected in Eq. (1) by the symmetry $f\to 1-f$ and $\tau^{\prime}\to-\tau^{\prime}$: interchanging species labels is equivalent to flipping the relative growth rate. Models that display this qualitative behavior fall into the Compact Directed Percolation (CDP) universality class, which is a special limit of DP Henkel et al. ; the nomenclature originates in the observation that the domains of growth are compact, as evidenced in Fig. 1$(c)$. Furthermore, as a result, the domain boundaries in a one-dimensional realization, undergo simple diffusion.

II.2 Flat growth with mutations

Adding mutations to the model changes the phenomenology markedly. Mutations break the symmetry between the species, allowing the mutants to appear within compact regions of the wild-type, as seen in Fig. 1$(e)$. In this case, the irreversible takeover of mutants is an absorbing state in which the wild-type cannot invade. There is now a transition between an active phase with wild-type bacteria, and an absorbing state without wild-type bacteria, at the point when the relative growth rate of mutants is balanced by the mutation rate. Near this critical point the fraction of wild-type bacteria will be small $f\approx 0$, justifying a small $f$ expansion to capture the long-time, large-lengths critical dynamics with the effective Langevin equation Henkel et al. ; Täuber et al. (2005); Täuber (2014)

[TABLE]

The active-to-absorbing phase transition is now controlled by the parameter $\tau$, which quantifies the difference in growth and mutation rates. Indeed, within the mean-field approximation, $\tau<0$ leads to a steady-state solution in the absorbing state $f^{*}=0$, which is continuously connected to an active phase when $\tau>0$ with steady-state value $f^{*}=\tau/2g$, whose magnitude is capped through $g$. As a minimal model for an absorbing phase transition, Eq. (2) is the field-theoretic Langevin equation capturing the critical behavior of the DP universality class Henkel et al. . It is worth noting that enforcing $g=\tau$, does not reimpose the CDP symmetry, and thus some care is required in interpreting the connection between universality classes at this level of description.

The DP hypothesis suggests that the critical behavior embodied by Eq. (2) is quite general. However, by allowing rough growth, thereby coupling the order parameter to the slow dynamics of the surface, we expect novel critical behavior to emerge.

II.3 Reflection asymmetry in rough front growth

Since the bacteria reproduce stochastically (not simultaneously in discrete generations) the growth front becomes rough. Ignoring overhangs and islands, the height of the front is capured by a single-valued height function $h({\bf x},t)$. For a single bacterial species the variations in the height of the front can then be expressed, in the spirit of a gradient expansion, as

[TABLE]

The first term corresponds to the uniform advance of the flat front, the last term $\xi({\bf x},t)$ is uncorrelated (in space and time) noise that captures the stochasticity of the growth steps, while the intermediate terms describe the dependence of smoothening and roughening on the local curvature and the slope of a rough front.

In the context of interface roughening, distinct scaling behaviors and universality classes are associated with the cases of $\lambda=0$ (Edwards-Wilkinson Edwards et al. (1982)) and $\lambda\neq 0$ (KPZ universality). The distinction is enforced by an underlying symmetry, in that fluctuations (in the moving frame with $v_{0}=0$) are symmetric under reflection $h\to-h$ for $\lambda=0$, with finite $\lambda$ breaking this symmetry. Thus, at this point we distinguish between two types of growth that we designate as reflection symmetric or reflection asymmetric. One way to distinguish between the two classes is to consider growth on an initially sloped surface with $\nabla h=\sigma$; Eq. (3) then predicts an average growth velocity $v({\bf\sigma})=\langle\partial_{t}h\rangle=v_{0}+\lambda\sigma^{2}/2+\cdots$. If such average for the bacterial film is independent of slope $\sigma$, then $\lambda$ (and higher order terms depending on slope) must be zero.

Note that Eq. (3) describes the height of the film measured from a flat surface. Thus if the film grows with the same speed $v_{0}$ in all directions (i.e. independently of slope $\sigma$), the component along the “vertical” direction increases as $v_{0}\sqrt{1+\sigma^{2}}\approx v_{0}+v_{0}\sigma^{2}/2+\cdots$. This simple geometric argument indicates that isotropic growth corresponds to $\lambda=v_{0}>0$. However, it is quite possible that due to the manner by which nutrients are supplied the growth velocity is not isotropic with a different dependence on $\sigma$. Indeed, for some growth models (see below) $\lambda$ can even have a negative sign.

II.4 Profile dependent competition without feedback

We next extend the above description of rough front growth to the case of wild-type and mutant bacterial species. As a first step in understanding the emergent behavior, we make the simplifying assumption that the two species are symmetric with respect to their affect on the front, i.e., variations in height of the film are still described by Eq. (3) with no feedback from $f({\bf x},t)$. However, we allow variations of the height profile $h({\bf x},t)$ to influence the bacterial field $f$, which again (in the limit of small $f$ with a gradient expansion) leads to the coupled equations

[TABLE]

Note that, with primary focus on overall shapes and fluctuations, we have set $v_{0}=0$. The effect of height fluctuations on the field $f$ is now captured by the three terms proportional to $\beta_{1}$, $\beta_{2}$, and $c$ whose interpretation will be explained in Sec. III. We emphasize that these are the only (possibly relevant) terms appearing at the lowest order in a gradient expansion. Their effect is to feed the surface fluctuations into the DP equation as an additional source of correlated (multiplicative) noise. A comparable set of equations was recently analyzed in order to understand how environmental fluctuations (here surface fluctuations) modify the DP transition, except the authors took $\beta_{2}=c=0$ and allowed for correlated noise in the KPZ equation, which lead to a somewhat distinct class of behaviors Sarkar and Basu (2012). With finite $\beta_{1}$ and $c$, evolution of the bacterial population is on average dependent on the slope of the initial front. We designate this case as profile dependent competition.

II.5 Profile dependent competition with feedback

Lastly, if the two species expand into the front at different rates as they grow, we expect the surface dynamics to be altered by the bacterial composition field. Indeed, in models of range expansions of immotile bacteria, differential expansion velocities find their origin in differing growth rates of the bacterial populations Kuhr et al. (2011); Farrell et al. (2017): whenever bacteria reproduce they replicate into nearby open space, thereby locally advancing the frontier. Thus, a selective growth advantage ($\tau\neq 0$) inherently couples the field $f$ to the surface height. Allowing for this feedback leads to our most general model of the rough DP transition,

[TABLE]

This full set of couplings is expected to lead to novel scaling behaviors. To the best of these authors’ knowledge this is the first time such a phenomenological model has been proposed.

In the following, we shall distinguish cases with ($\alpha\neq 0$) and without ($\alpha=0$) feedback; with profile independent ($\beta_{1}=\beta_{2}=c=0$) or dependent ($\beta_{1}$, $\beta_{2}$, $c$  non-zero) competition; and with ($\lambda\neq 0$) or without ($\lambda=0$) front reflection symmetry. Each combination may in principle belong to a distinct universality classes. Note that full reflection symmetry $h\to-h$ in the coupled equations would require $\beta_{1}=\beta_{2}=0$, while $c$ may still be non-zero.

III Macroscopic manifestations of coupling parameters

In this section, we discuss and illustrate the macroscopic manifestations of the phenomenological parameters introduced to describe coupling of bacterial composition to the surface shape. While some results of this section, such as the dependence of surface morphology on $\lambda$, are well known in the surface-growth community, they may be less familiar in the biological range-expansion circle, and are therefore included for completeness. The numerical illustrations in this section were generated using a modification of the Domany-Kinzel (DK) cellular automaton Korolev et al. (2010); Lavrentovich et al. (2013) (modified “brick wall” solid-on-solid model Drossel and Kardar (2000, 2003), or “single step” model Meakin et al. (1986); Plischke et al. (1987).), whose elementary steps are illustrated in Fig. 1$(b)$. Growth proceeds on a lattice of transverse extension $L$ with periodic boundary conditions. New sites are added randomly as long as the neighboring sites in the layer below are occupied to prevent overhangs. An added site is assigned one of two ‘colors’ based on the colors of its local neighborhood through three microscopic parameters: Briefly, they are $\tau_{\mu}$ which biases the identity of newly added sites based on the selective differences of the colors in the layer below, $\mu_{\mu}$ a microscopic mutation rate, $\alpha_{\mu}$ which biases the height growth rate based on the composition field, and $\beta_{\mu}$ which couples the slope to the composition field by allowing the active neighbors on the same layer to affect the identity of the newly added site. The role and implementation of these (microscopic) parameters is described in more detail in the following paragraphs with additional specifics in Appendix A.

III.1 Dependence of front velocity on slope

In our model, occupation of new sites is only possible when both neighbors in the layer below are occupied. This causes growth to be slower on sloped fronts where it is harder to find occupied neighbors, suggesting $\lambda<0$ for this model Drossel and Kardar (2003). Indeed, for the equivalent “single step” model Meakin et al. (1986); Plischke et al. (1987) an exact result by Krug et al. Krug et al. (1992) indicates that $v=(1-\sigma^{2})/2$ for such discrete models. The KPZ equation maintains the average slope of the front: If the initial profile has slope $\nabla h(x,t)=\sigma$, then Eq. (3) predicts a constant growth velocity $v=\langle\partial_{t}h\rangle=v_{0}+\lambda\sigma^{2}/2$. Matching the discrete and continuum results than leads to $v_{0}=1/2$ and $\lambda=-1$. We verify this behavior by initializing our simulation with a fixed slope, as depicted in the inset of Fig. 2.

By then measuring the growth velocity $v$ for various slopes $\sigma$, we observe a quadratic dependence, and confirm a negative value of the nonlinearity $\lambda$, as demonstrated in Fig. 2.

III.2 Influence of front profile on composition

The couplings proportional to $\beta_{1}$, $\beta_{2}$ and $c$ encapsulate how height fluctuations alter the selective pressure between the two species, leading to composition field dynamics that depend locally on height variations.

The effect of the first nonlinear coupling parameter $\beta_{1}$ is most dramatic on domain boundaries between wild-type and mutant species ($\nabla f\neq 0$). Like for the KPZ nonlinearity, if the height field has a fixed slope at any time $\nabla h(x,t)=\sigma$, the $\beta_{1}$ term induces a steady drift of the domain boundary with a magnitude proportional to $\beta_{1}\nabla f\cdot\sigma$. For $\beta_{1}>0$ the boundary will tend to move down a slope, and vice versa for $\beta_{1}<0$.

In order to highlight this effect, we allow for a microscopic parameter $\beta_{\mu}$ in our numerics that biases the identity of newly added occupied sites based on neighbors on the same level Drossel and Kardar (2000, 2003). Specifically, the identity of a newly occupied site is determined by randomly copying the identity from either of the neighbors in the layer below or from one of the neighbors in the same layer, the latter weighted by a factor $\beta_{\mu}$. Since on slopes we are more likely to have same-layer neighbors, this accentuates the coupling between the composition field and the height fluctuations.

In Figs. 3$(a)$ and $(b)$, we highlight the influence of $\beta_{\mu}\neq 0$ with two exaggerated examples obtained using sloped initial conditions and a 50/50 split of active wild-type and inactive mutant sites.

For $\beta_{\mu}>0$, we enhance the drift of the domain boundary down the slope; whereas, by taking a sufficiently negative $\beta_{\mu}$, we can cause the domain boundary to move up the slope. We can quantify the macroscopic effect of the model parameter $\beta_{\mu}$ by measuring the dependence of the domain boundary velocity on slope as in Fig. 3$(c)$.

It is important to recognize that even in the absence of an explicit microscopic slope to composition coupling ($\beta_{\mu}=0$), the modified DK algorithm still causes any domain boundary to drift with a nonzero speed, i.e., the microscopic model has an emergent nonzero macroscopic coupling $\beta_{1}$ as anticipated on symmetry grounds in Eq. (II.5). We can measure its magnitude, by ramping up small negative $\beta_{\mu}$ until the mean drift goes away. We see this in Fig. 3$(d)$ where the mean domain boundary drift was determined from the average of 250 realizations initialized as in Figs. 3$(a)$ and $(b)$. For $\beta_{\mu}=0$, the mean domain boundary drift is very small and nonzero. By making $\beta_{\mu}$ more negative, we see the mean domain boundary drift passes through zero. Unfortunately, due to the expense of the simulations our statistics are very poor. Thus, determining the precise crossing point is problematic, but it is reasonable to conclude that its value is in the interval $[-0.005,-0.002]$.

The remaining nonlinear coupling terms in Eq. (II.5) are proportional to the composition field $f$. We can interpret their impact by recognizing that they combine to act as an effective control parameter

[TABLE]

The above equation implies that on a sloped surface, the location of the DP transition shifts from that on a flat surface by an amount that is quadratic in slope. A fit to $\tau_{c}(\sigma)$ can then be used to estimate the macroscopic parameter $c$. This is reminiscent of depinning of an interface in anisotropic random media in which the quadratic dependence of critical force on slope can be used to identify a KPZ-type nonlinearity Tang et al. (1995). It is difficult to propose a similar procedure for obtaining the parameter $\beta_{2}$ as curvature changes during the growth process. Variations in composition along a corrugated surface may provide a means of macroscopically accessing this parameter.

The limit $\beta_{1}=\beta_{2}=\beta$ for $c=0$ is interesting, as the terms combine to form an advection, $\beta\nabla(f\nabla h)$, of the field $f({\bf x},t)$ with local velocity proportional to the slope $\nabla h({\bf x},t)$. Variations in the integrated composition are then not affected by the shape of the surface, although its dynamics can be modified.

III.3 Dependence of front velocity on composition

The term proportional to $\alpha$ in Eq. (6) represents how differences in expansion rates between active wild-type sites and inactive mutant sites lead to variations in the growth velocity of the film height. To unravel this effect, we implement in our numerical model the microscopic parameter $\alpha_{\mu}$, which biases the addition of new occupied sites based on the composition of the local neighborhood. Specifically, depending on the number of mutant species in the layer below $m_{\rm s}=\{0,1,2\}$, we reject the addition of a new occupied lattice site with probability $m_{\rm s}\alpha_{\mu}$. The outcome of an exaggerated implementation with a rather large $\alpha_{\mu}$, starting with flat initial conditions and a bounded region of active wild-type sites, is presented in Fig. 4$(a)$.

As time goes on, the preference to add new active sites results in a bulge that quickly assumes a triangular shape, while the front above inactive sites remains almost flat. The finite slope of the triangular region originates in the balance between increased growth speed $\alpha$ with the reduction $\lambda\sigma^{2}/2$ (since $\lambda<0$) due to slope, leading to $\sigma^{2}=-2\alpha/\lambda$. Thus the macroscopic parameter $\alpha$ can be estimated from this slope. For example, in the realization presented in Fig. 4, we find $\alpha\approx 0.04$. Note that the above argument relies on a negative value of $\lambda$. To visualize the $\lambda>0$ case, we numerically solved the coupled DP/KPZ stochastic partial differential equations [Eqs. (6) and (II.5)] with $\lambda=1$. The result is presented in Fig. 4$(b)$. Here, the faster island remains roughly flat, with triangular dips in the slower domains left behind 111For $\lambda=0$, the bulges and dips will have parabolic shapes..

Clearly, differential expansion rates lead to bulges (or dips) in the front. The proliferation of such protrusions when both species coexist is expected to increase the roughness of the front, due to greater variation in the height of the jagged domains. Indeed, this effect has been previously identified in both an agent-based model of the mechanical expansion of reproducing bacteria Farrell et al. (2017), as well as in a two-species Eden Model with mutations Kuhr et al. (2011), both of which tie the expansion rate to the replication rate.

IV Genealogical paths in deterministic growth

Additional insight is obtained by examining the deterministic limit. While much of this section recapitulates known results, by reanalyzing them in this new context we reveal interesting connections between fitness and shape. Thus, this section can be safely skipped by readers more interested in the scaling analysis.

Our analysis in this section is based on unraveling a hidden connection between surface growth and populations dynamics in the absence of noise. To this end, it will prove beneficial to look at a more general model of competition and mutations in biological populations without noise, the Eigen model Eigen (1971); Eigen and Schuster (1977); Eigen et al. (1988, 1989), whose properties and solution are well documented Leuth ausser (1986, 1987); Domingo and Schuster (2016). However, to make a connection to the scenario at hand, we here first consider a variant of the model that includes stochasticity and allows for migration of the populations in space, before returning to the special limit of deterministic growth.

Consider quite generality a population composed of $N_{i}({\bf x},t)$ individuals of $i=1,2,\dots,n$ species at location ${\bf x}$ at time $t$. The populations change in time according to

[TABLE]

where $\eta_{i}$ are independent, delta-correlated white noises, $\alpha_{i}$ denotes the fitness (reproductive capacity) of species $i$, and the matrix elements $M_{ij}$ (with $\sum_{i}M_{ij}=0$) indicate the rates of mutation from $j$ to $i$. Equation (10) assumes equal rates of diffusion $\nu$ (migration) and demographic noise $D$ for all species. Note that the overall population size is not fixed and may grow exponentially, or decay to zero. Indeed, the total population $N({\bf x},t)=\sum_{i}N_{i}({\bf x},t)$ at location ${\bf x}$ changes according to

[TABLE]

where $\overline{\alpha}=\sum_{i}\alpha_{i}N_{i}/N$ is the average fitness over the population and $\eta$ is delta-correlated, white noise.

We can reformulate Eqs. (9) and (10) in terms of the dynamics of composition fields, as well as an effective height field. This will allow us to tweeze apart the effects of topography on fitness, in the deterministic limit. To this end, we form population fraction fields $f_{i}({\bf x},t)=N_{i}({\bf x},t)/N$, akin to our composition field, that dynamically couple to what we will designate as a ‘height field’ $h({\bf x},t)=\frac{2\nu}{\lambda}\ln N$ (which is recognizable as the Cole-Hopf transformation). Now, using Ito’s Lemma we find for the variation of the composition fields

[TABLE]

whose stochasticity can be represented with correlated multiplicative noise $\langle\eta^{\prime}_{i}(t)\eta^{\prime}_{j}(s)\rangle=(f_{i}\delta_{ij}-f_{i}f_{j})\delta(t-s)$, and

[TABLE]

The stochastic element in Eq. (9), proportional to $\sqrt{N_{i}}$, is due to random reproduction events (demographic noise). Interestingly, this type of noise goes to zero in Eqs. (11) and (12) for large population size, in the limit $N\to\infty$. In this limit, Eq. (12) becomes the deterministic KPZ equation Kardar (2007)

[TABLE]

and Eq. (11) simplifies to

[TABLE]

a noiseless, multi-component version of Eq. (II.5) with $K=\nu$, $\beta=\lambda$ and $\beta_{2}=c=0$. Indeed for two species ($n=2$) with $f_{1}=f$ and $f_{2}=1-f$, $\overline{\alpha}=\alpha_{1}f+\alpha_{2}(1-f)$ and mutations only from 1 to 2 at rate $M_{11}=-\mu$, we recover Eq. (II.5) with $\tau=\alpha_{1}-\alpha_{2}-\mu$ and $g/2=\alpha_{1}-\alpha_{2}$; also $\alpha=2\nu(\alpha_{1}-\alpha_{2})/\lambda$ in Eq. (6). Thus, within the Eigen model we find a special limit of the competition dynamics we have been studying.

The Cole-Hopf transformation maps the noiseless (deterministic) but nonlinear Eqs. (13) and (14) for $h$ and $f_{i}$ into a set of linear equations that are easily solved Kardar (1999). Given an initial population of composition $f_{i}({\bf x},0)$ distributed on a landscape characterized by a profile $h({\bf x},0)$ as in Fig. 5, we can construct ‘populations’ $N_{i}({\bf x},0)=f_{i}({\bf x},0)\exp\left[\frac{\lambda}{2\nu}h({\bf x},0)\right]$. The deterministic version of Eq. (9) can now be easily integrated to

[TABLE]

where the matrix $G_{ij}=M_{ij}+\alpha_{i}\delta_{ij}$ combines the effects of mutation and fitness.

This mapping through Eq. (IV) illustrates the importance of inhomogeneous initial conditions, e.g., uniformly inoculating a non-flat surface with height profile $h({\bf x},0)$, or a flat surface with a non-uniform population number $N({\bf x},0)$. In the simplest case of non-competitive growth of two species ($\alpha_{1}=\alpha_{2}$) with no mutations, i.e., for $G=0$, the form of the front at long times (alternatively, for $\nu\to 0$) can be extracted from Eq. (IV) using the saddle point approximation, as

[TABLE]

In this limit, the surface profile (for $\lambda>0$) is a collection of parabolic mounds, each centered on a local peak, say at ${\bf x_{0}^{\prime}}$ of the initial profile. Equation (IV) then implies $f_{i}({\bf x},t)=f_{i}({\bf x_{0}^{\prime}},0)$, i.e., that the entire “mound” is populated by decedents of the quasi-species member residing at the initial local peak ${\bf x_{0}^{\prime}}$. The advantage conferred by the initial profile in this model is equivalent to an increased population size (on a flat surface) proportional to $\exp[\frac{\lambda}{2\nu}h({\bf x},0)]$. As a result, quasi-species at the height of a local maximum for $\lambda>0$, or a minimum for $\lambda<0$, slowly take over neighboring regions as the coarsening progressively spreads out.

The emergence of parabolic mounds, as well as the benefits of initial location to survival, are illustrated in Fig. 5. Here, the species starting at the highest peak is in fact the less fit ($\alpha_{\rm blue}<\alpha_{\rm gray}$), and would disappear with a flat starting profile. However, starting from a high peak, it actually expands at the expense of the fitter species at lower height, establishing a niche within the growing parabola.

More generally, the combined diffusion/growth/mutation kernel in Eq. (IV) can also be cast as an integral over paths,

[TABLE]

The first (diffusion) component of the kernel is dominated by paths that proceed along a straight line from ${\bf x}^{\prime}$ to ${\bf x}$ in time $t$, with ‘velocity’ ${\bf v}=d{\bf x}/dt=({\bf x}-{\bf x}^{\prime})/t$. For $G=0$, as in Fig. 5, such straight paths join points on a parabolic segment to the corresponding local maximum, indicating the most likely genealogical history for the species at the front. In the presence of mutations and fitness variations, genealogical tracks are composed of segments $i=1,2\cdots,n$. In the $m^{th}$ segment, from $t_{m-1}$ to $t_{m}$, ancestry is of quasi-species $i(m)$, with $i(0)=j$ and $i(n)=i$. The first factor in the third line of Eq. (17) represents the increase in population size according to the fitness within each segment, the second factor accounts for the mutations probabilities between segments.

The path integral formulation can be used to explain the jagged profiles observed in Fig. 4. Consider two species starting from a flat profile as in this figure, with initial separation point marked as $x=0$, with different fitness values $\alpha_{1}>\alpha_{2}$. The optimal genealogical track terminating at $(x,t)$ in quasi-species 2 may actually correspond to ancestry starting with quasi-species 1 at $x=0$, proceeding vertically along $x=0$ for time $t-\tau$, followed by a portion that proceeds with velocity $v=x/\tau$ to the target point in the remaining time $\tau$ after transforming to quasi-species 2. The vertical segment gains weight $\alpha_{1}(t-\tau)$ in Eq. (17), the sloped segment has cost $\alpha_{2}\tau-{x^{2}\over 4\nu\tau}$. The resulting height profile is obtained by optimizing the sum of the two contributions as

[TABLE]

where we have noted $\alpha=2\nu(\alpha_{1}-\alpha_{2})/\lambda$, and set $v_{1}=2\nu\alpha_{1}/\lambda$. There is an optimum at $\tau^{*}={|x|\over\sqrt{4\nu(\alpha_{1}-\alpha_{2})}}$, which upon substitution in the above equation leads to a slope of $\sigma^{2}=h^{2}/x^{2}=|2\alpha/\lambda|$, consistent with the result in Sec. III.3.

V Special limits and symmetries

The remainder of this paper focuses on universal scale-invariant properties of Eqs. (6) and (II.5). However, for certain choices of parameters these equations posses symmetries, providing additional constraints and insights on their general properties. While this information is familiar knowledge from earlier studies on interface roughening and directed percolation, it is nevertheless valuable to present the results in this context.

V.1 ‘Galilean invariance’ and scaling of roughness

There is no characteristic (microscopic) scale limiting the extent of roughness, which thus assumes scale-invariant characteristics. For a single-valued height function, this leads to self-affine scaling, and dynamic fluctuations characterized by the scaling form Kardar (2007)

[TABLE]

where $\chi$ is the roughness exponent (later also designated as $y_{h}$), and $z_{h}$ is a dynamic scaling exponent, measuring the relative scalings in space and time. For the one dimensional KPZ equation, $\chi=1/2$ and $z_{h}=3/2$ are exact results Kardar et al. (1986), while in two dimensions, the estimates $\chi_{h}\approx 0.4$ and $z_{h}\approx 1.6$ come from numerical simulations Halpin-Healy and Takeuchi (2015).

In the absence of the coupling to the composition field, the KPZ equation’s scaling exponents are not independent, due to a hidden symmetry, which is a remnant of the rotational symmetry of the gradient expansion leading to the KPZ equation. Defining a velocity field related to the slope by ${\bf v}=-\lambda{\bf\nabla}h$ Kardar et al. (1986); Frey and Täuber (1994), maps Eq. (6) (for $\alpha=0$) to the Burgers equation Burgers (1974)

[TABLE]

which describes a randomly stirred fluid. The significance of this mapping is that the fluid flow equation is (statistically) invariant under a Galilean transformation to a moving frame (${\bf v^{\prime}}({\bf x^{\prime}},t^{\prime})={\bf v}({\bf x},t)+{\bf u}$ with ${\bf x^{\prime}}={\bf x}+{\bf u}t$ and $t^{\prime}=t$). In terms of the KPZ equation, this transformation corresponds to a tilt

[TABLE]

with ${\bf u}=-\lambda{\bf\epsilon}$. The significance of the Galilean symmetry is that it constrains the scaling of ${\bf v}$ and hence ${\bf\nabla}h$ in terms of those of space and time (as ${\bf x}/t$), leading in Eq. (19) to the exponent identity $\chi=2-z_{h}$. Formally, since $\lambda$ appears both in the symmetry mapping of Eq. (21) and as the nonlinearity in the KPZ equation, any renormalization scheme that maintains the symmetry cannot renormalize $\lambda$ (at any order).

An infinitesimal transformation as in Eq. (21), leads to $\nabla_{\bf x}=\nabla_{\bf x^{\prime}}$ and $\partial_{t}=\partial_{t}^{\prime}-\lambda{\bf\epsilon\cdot\nabla_{\bf x^{\prime}}}$ with $\nabla_{\bf x}h=\nabla_{\bf x^{\prime}}h^{\prime}-{\bf\epsilon}$. It is easy to check that these transformations also leave Eq. (II.5) invariant if $\beta_{1}=\lambda$ and $c=0$ (for all values of $K$ and $\beta_{2}$). With $\beta_{1}=\lambda$, both the composition and front–slope fluctuations are transported along the front with the same speed and maintained in phase. The dynamics of composition fluctuations should thus mirror those of the front, when there is rotational symmetry.

V.2 Scaling in directed percolation

The transition to an absorbing state at a critical value of $\tau_{c}$ is accompanied by a number of singularities. These include divergences of the correlation lengths in space $\xi_{\perp}$ – perpendicular to the direction of growth – and time $\xi_{\parallel}$ – parallel to the growth – as well as, the vanishing of density of particles in the active phase $\rho$, and the probability of a single active seed percolating through space $P_{\rm perc}$. In the limit of $\delta\tau=\tau-\tau_{c}\to 0$, the scaling of these quantities are characterized by critical exponents as Henkel et al.

[TABLE]

A dynamic critical exponent, relating the scalings in the two directions, is defined by $z_{f}=\nu_{\parallel}/\nu_{\perp}$. The values of these exponents have been obtained previously from a renormalization group analysis near the critical dimension $d_{c}=4$, perturbatively in $\epsilon=4-d$, and are listed in Table 1 Täuber et al. (2005); Täuber (2014). In addition for comparison, we have listed values at $d=1$, which were calculated using a series expansion Jensen (1996).

The DP equation also satisfies a form of time-reversal symmetry known as rapidity-reversal  Henkel et al. . It is manifest only at the level of the dynamic action, which we will discuss in the following sections. Briefly, if we introduce the response field ${\tilde{f}}(x,t)$ conjugate to the order parameter $f(x,t)$, then rapidity-reversal is a symmetry under the exchange ${\tilde{f}}(x,t)\longleftrightarrow-f(x,-t)$ of the dynamic action. This symmetry has two important consequences. The first is that it imposes the exponent identity $\beta=\beta^{\prime}$: the disappearance of a cluster of active wild-type sites is the time reverse of the percolation of a cluster of active wild-type sites from one active seed Henkel et al. . Second, on a more practical level, this symmetry forces $g$ and $D_{f}$ to renormalize together, implying that there is actually only one relevant parameter $u\propto gD_{f}$. Since rapidity-reversal is a relative of time-reversal symmetry – a symmetry not shared by the KPZ equation –, we do not expect this symmetry to remain for our equations, as manifested in $g$ and $D_{f}$ renormalizing differently with unequal $\beta$ and $\beta^{\prime}$ exponents.

VI Critical exponents and renormalization

We now turn to the analysis of the universal scaling behaviors at the active-to-absorbing phase transition in the presence of surface fluctuations.

VI.1 Setup

Our renormalization group (RG) analysis is based on the Martin-Siggia-Rose-Janssen-de Dominici (MSRJD) dynamic functional Täuber (2014), which plays the role of the Hamiltonian in a nonequilibrium renormalization scheme. To this end, we determine the probability distribution of the space-time trajectories of $f$ and $h$ by noise-averaging $\delta$-functions that enforce the equations of motion [Eqs. (6) and (II.5)]

[TABLE]

where the average is over Gaussian noise

[TABLE]

The $\delta$-functions are then replaced with their generalized Fourier representations by introducing conjugate (response) fields ${\tilde{h}}$ and ${\tilde{f}}$, allowing us to cast the probability distribution into the form

[TABLE]

The MSRJD function naturally splits as ${\mathcal{J}}={\mathcal{J}}_{0}+{\mathcal{J}}_{\rm int}$, into a solvable Gaussian contribution

[TABLE]

and a piece containing all higher-order interactions

[TABLE]

The first step in our scaling analysis is to obtain the naive scaling dimensions of all the parameters. To this end, we rescale space, time, height and order parameter according to $x\to bx$, $t\to b^{z}t$, $f\to b^{y_{f}}f$, $h\to b^{y_{h}}h$, ${\tilde{f}}\to b^{{\tilde{y}}_{f}}{\tilde{f}}$, and ${\tilde{h}}\to b^{{\tilde{y}}_{h}}{\tilde{h}}$, and then require the invariance of the Gaussian contribution ${\mathcal{J}}_{0}$. This procedure leads to the mean-field values $z=2$, $y_{f}={\tilde{y}}_{f}=-d/2$, $y_{h}=2-d/2$ and ${\tilde{y}}_{h}=-2-d/2$. With these scalings, all nonlinearities acquire a scaling dimension $2-d/2$, with a critical dimension $d_{c}=4$. However, if $\alpha\neq 0$ the noise in the KPZ equation $D_{h}$ scales as $d+z+2{\tilde{y}}_{h}=-2$ and is irrelevant in any dimension, and will be ignored when this is the case. Furthermore, the nonlinear coupling proportional to $c$, while possibly relevant, is higher order in the sense of not renormalizing to one-loop order, and we shall also ignore it.

To carry out the RG calculation, we employ Wilson’s momentum shell approach Täuber (2012); Kardar (2007). This is accomplished by first integrating out short wavelength (high momentum) modes $\Lambda/b<|{\bf q}|<\Lambda$, where $\Lambda$ is the wave vector cutoff and $b>1$ is a scaling factor. We then rescale space, time, height and order parameter as before. This requires treating the nonlinearities perturbatively by expanding $e^{-{\mathcal{J}}_{\rm int}}$ in powers of the couplings. Finally, taking $b=e^{l}$ with $l\ll 1$, leads to the differential RG equations for the model parameters. To be consistent with the perturbation expansion, we have to focus on dimensions near the critical dimension, so we set $d=4-\epsilon$ with $\epsilon\ll 1$. This technique is standard Täuber (2012); Kardar (2007), so we do not go into more detail here, simply collecting the results in the Appendix B.

Before turning to the analysis of the fully coupled dynamics, we will first examine two limiting cases with unidirectional coupling between the fields that allow us to highlight the subtle phenomenology involved in this nonequilibrium critical transition.

VI.2 Profile independent competition, rough growth front, with feedback

We first consider profile independent competition, where surface dynamics does not influence the DP order parameter:

[TABLE]

In this case, the absorbing phase transition is unaltered by the height field, though the presence of order-parameter fluctuations is expected to alter the scaling properties of the height field, with the composition field acting as a source of correlated noise Medina et al. (1989). The general RG equations presented in the Appendix B then reduce to the following

[TABLE]

Above the DP critical dimension, $d>4$, simple dimensional analysis of the DP equation yields $y_{f}=-d/2$ and $z=2$. The choice of $\chi=(2-d/2)<0$ corresponds to the Gaussian fixed point for the field $h$, at which small $\lambda$ is irrelevant and flows to zero. There is also, however, an unstable fixed point with $\chi=0$ for $(\lambda\alpha)^{*}\propto(d-4)>0$. Larger values of $(\lambda\alpha)$ lead to a rough phase that is not perturbatively accessible. While a similar transition exists for white noise in the KPZ equation (for $d>2$), the difference is that here the transition is controlled by the product $(\lambda\alpha)>0$.

We may anticipate that after the unstable fixed point collides with the Gaussian fixed point at $d=4$, a stable fixed point emerges for $d<4$. However, the scaling exponents at the DP fixed point (listed in Table 1 Täuber et al. (2005); Täuber (2014)) are $z_{f}=2-\epsilon/12$ and $y_{f}=-2+7\epsilon/12$. This value of $z$ in Eq. (30) pushes $\nu$ towards zero leading to a divergence of the correction term in Eq. (32).

The complexity of this situation is best understood by first considering the limit $\lambda=0$. In this case, a perturbation of the field $f$ at length scale $\ell$ will persist for a time $\tau_{f}(\ell)\sim A_{f}\ell^{z_{f}}$. The amplitude $A_{f}$ is modified from the Gaussian value of $1/K$ by renormalization of non-linear couplings to $A_{f}\sim\frac{1}{K}\left(\frac{gD_{f}}{K^{2}}\right)^{2-z_{f}}$. With $\lambda=0$, the response of the height field to the perturbation ($\alpha f$) is linear, and characterized by a time-scale $\tau_{h}\sim\ell^{2}/\nu$. The coupled equations are thus characterized by two time scales that are matched at a non-universal length scale $\ell_{c}\sim(\nu A_{f})^{\frac{1}{2-z_{f}}}$. For $\ell\ll\ell_{c}$ height perturbations are smoothed out by the Laplacian term in the front equation, while distortions at scales $\ell>\ell_{c}$ are created and removed rapidly by the field $f$ of critical DP. In each case, we can associate a roughness for the profile (by setting $d\alpha/dl=0$) as $\chi=z+y_{f}$, leading to $\chi=7\epsilon/12$ for $\ell<\ell_{c}$ ($z=2$), and $\chi=\epsilon/2$ for $\ell>\ell_{c}$ ($z_{f}=2-\epsilon/12$). This possibility where dynamic scaling exponents between fields can be different is known as weak dynamic scaling – to be distinguished from strong dynamic scaling when the exponents are the same. Such an effect is not unusual when coupling multiple fields together Folk and Moser (2006) and has been observed in Models C and D Akkineni and Täuber (2004), kinematic waves Das et al. (2001), active-to-absorbing phase transitions in fluctuating environments Sarkar and Basu (2012, 2014), as well as for conserved KPZ equation coupled to critical fields Banerjee and Basu (2018).

For $\lambda\neq 0$, the choice of $z=2$ and $\chi=0$ indicates that a small $\alpha$ is relevant, and the second order term in Eq. (32) has to be considered. Now for $(\lambda\alpha)<0$ a stable fixed point appears with the above exponents. For $(\lambda\alpha)>0$ the coupling runs to large values, making it impossible to assign effective exponents to the regime $\ell<\ell_{c}$ by the current analysis.

VI.3 Profile dependent competition, rough growth front, without feedback

Whereas the previous section considered the field $f$ evolving without influence from the front profile, here we consider that the dynamics of profile $h$ gets no feedback from the composition field while affecting its variations according to

[TABLE]

We have included the factor proportional to $c$. Even though we do not consider it in the full RG treatment; the following conclusions nevertheless apply.

As far as the perturbative RG analysis is concerned, the behavior of the above coupled equations are either trivial or inaccessible. This is because in dimensions $d>2$, the nonlinearity in Eq. (33) is either irrelevant leading to roughness exponent $\chi=(2-d)/2<0$, or relevant but flowing to uncontrolled large values. In the latter case perturbative RG in not useful, leaving numerical simulations as a promising approach Kuhr et al. (2011). In the former case, negative values of $\chi$ imply that the three non-linear couplings $(\beta_{1},\beta_{2},c)$, involving factors of $h$ or $h^{2}$, in Eq. (34) are irrelevant. With all non-linearities eventually disappearing at large scale, the behavior then mimics uncoupled roughness and DP.

VI.4 Full renormalization group analysis

Having discussed how unidirectional coupling affects the scaling properties of the DP transition, or the surface roughness, we now consider an RG analysis of the fully coupled rough DP equations in $d=4-\epsilon$.

The RG equations turn out to only depend on the following reduced set of dimensionless effective parameters, whose analysis is substantially simpler:

[TABLE]

with $A_{4}=S_{4}/(2\pi)^{4}$ where $S_{4}$ is the surface area of a 4-dimensional sphere and is included for convenience. In the first row, $\theta$ measures the relative relaxation times of the DP equation to the KPZ equation, and $\bar{\tau}$ is the dimensionless control parameter for the DP transition. The remaining parameters are the dimensionless combinations of nonlinear couplings. Importantly, by dividing out the $K$ and $\nu$ dependencies, they capture solely their effects on spatial variations, leaving the relative temporal dynamics encoded in $\theta$.

The fixed point structure further is controlled by the value of the time-scale ratio $\theta=K/\nu$. Since $\theta$ is not a nonlinear coupling, it need not be perturbatively small. Thus, there are three possible phenomenologically distinct fixed points corresponding to $\theta=0$, $\theta\to\infty$, or with a nonzero finite value at the fixed point. We address each possibility in turn.

VI.4.1 $\theta\to\infty$ sector

This limit represents the situation where $K\gg\nu$, which signifies that the DP dynamics evolve on a much faster time scale than the KPZ equation. Put another way, from the perspective of surface fluctuations, the DP equation relaxes instantaneously to its probabilistic stationary state conditioned on the surface configuration.

In this limit, the RG equations for the nonlinear couplings reduce to

[TABLE]

There are two fixed point solutions to these equations, where all derivatives are zero. The first is the Gaussian fixed point: $u^{*}=\zeta^{*}=\delta_{1}^{*}=\delta_{2}^{*}=0$. This corresponds to the uncoupled ‘mean-field’ behavior, which is stable above the critical dimension $d_{c}=4$, but unstable below. The second solution corresponds to the standard DP fixed point: $u^{*}=2\epsilon/3$ and $\zeta^{*}=\delta_{1}^{*}=\delta_{2}^{*}=0$ (cf. Table 1). Here, the surface fluctuations are driven by the DP fluctuations, but not vice versa. As in Sec. VI.2, we can explore different scalings and responses of height fluctuations at short and long scales; however, the couplings ($\beta_{1}$ and $\beta_{2}$) are relevant, rendering the standard DP fixed point unstable.

VI.4.2 $\theta$ zero sector

When $\theta=0$, we have the reverse situation where $K\ll\nu$, and the surface dynamics evolve much faster than the DP composition field. In this limit, the RG equations for the dimensionless nonlinear couplings become

[TABLE]

Using a symbolic Mathematics package, we have found 12 fixed points for these equations, including the Gaussian fixed point and the uncoupled DP fixed point (cf. Sec. VI.2).

Of all the fixed points, a linear stability analysis reveals only one to be both stable and in the physical parameter regime ($u>0$) 222There is one fixed point for which all nonlinear couplings are stable and $\theta$ is marginally stable ($d\theta/dl|_{\theta=0}=0$), but flows $\theta\to\infty$ for any finite value of $\theta$: $u=0$, $\zeta=2\epsilon$, $\delta_{1}=-2\epsilon$, and $\delta_{2}=0$. To verify the stability of this fixed point, higher order terms need to be considered., with

[TABLE]

The accompanying scaling exponents are obtained by demanding scale invariance of the remaining parameters, whose RG equations at this fixed point are

[TABLE]

combined with the algebraic constraints on the field exponents

[TABLE]

As discussed in Sec. VI.2, we have a weak dynamic scaling fixed point with two possible choices for the dynamic scaling exponent, $z=2$ or $2+\epsilon/6$, both of which are consistent with the requirement that $\theta=K/\nu\to 0$ is a stable fixed point. If we choose $z=2$, then $\nu$ is scale invariant ($d\nu/dl=0$), and $K$ decays as $K\sim e^{-\epsilon l/6}$. However, this behavior is at odds with the underlying premise of RG that we isolate the scale-invariant behavior that emerges on the longest length and time scales. With this choice of $z$, the RG procedure ends up holding the dynamics on the time scale of $\sim 1/\nu$ fixed. But in so doing, the RG procedure commensurately keeps slowing down the DP dynamics controlled by the time-scale $\sim 1/K$. This is a consequence of the the fact there is actually a slower time scale controlled by $K$. Thus, this choice of $z$ does not focus us on the slowest time scale, where the appropriate choice is $z=2+\epsilon/6$. We should note that as observed in Sec. VI.2, this choice does not rule out another dynamic scaling exponent that controls the response to perturbations on shorter wavelengths.

With this caveat in mind, the rest of the scaling exponents now follow by making $D_{f}$ and $\alpha$ [Eqs. (49) and (50)] scale invariant, coupled with the algebraic constraints in Eqs. (51) and (52), as well as extracting the growth rate of $\bar{\tau}$ from Eq. (45); leading to

[TABLE]

However, the exponents most often measured are those defined in Eqs. (19) and (22). They are obtained through the scaling relationships: $\nu_{\perp}=1/y_{\tau}$, $\beta=-y_{f}/y_{\tau}$, $\beta^{\prime}=-{\tilde{y}}_{f}/y_{\tau}$, and $\chi=y_{h}$ Henkel et al. .

We have collected our results in Table 1. Included there for comparison are values obtained from a field-theoretic analysis of the uncoupled DP phase transtion in $d=4-\epsilon$ Henkel et al. , as well as exponents measured in simulations of a two-species Eden Model in $d=1$ Kuhr et al. (2011), which is expected to be in the DP/KPZ universality class.

VI.4.3 $\theta$ nonzero and finite sector

The full set of RG equations are sufficiently complex to make obtaining fixed points with $\theta\neq 0$ quite challenging. To the best of our knowledge, there are no more stable fixed points.

VI.4.4 Discussion

Within the perturbative approach here, we have found a single stable fixed point and extracted its scaling behavoir. To access the phenomenology of the dynamics near this fixed point, we can write the associated Langevin equations by keeping only terms that have nonzero fixed point values:

[TABLE]

with the important proviso that $\alpha$ and $\beta\equiv\beta_{1}=\beta_{2}$ have opposite signs ($\delta_{j}<0$).

Interestingly, at this fixed point $\theta=0$, suggesting that the relaxation dynamics of the KPZ equations are very fast $\nu\gg K$ (with $\alpha/\nu$ flowing to a finite value), which we interpret as the height dynamics ‘instantaneously’ relaxing to match the DP order parameter configuration. Another point of note is that the coupling of the height back into the order parameter enters as an advection term with effective flow field $\nabla h$: this advective term pushes the active wild-type sites to the maxima or minima of surface undulations depending on the sign of $\beta$. In either case, however, the effects of $\alpha$ and $\beta$ counteract each other. For example, suppose $\alpha>0$ and $\beta<0$. The positive $\alpha>0$ would cause regions of dense active wild-type sites to bulge out, creating regions of negative curvature. However, $\beta<0$ then drives domain boundaries up the slopes condensing the region of active wild-type sites, thereby slowing the surface growth. In effect, washing out surface undulations. This picture is consistent with a negative roughening exponent $\chi=-\epsilon/6$ and a flat front.

Given a configuration $f({\bf x},t)$, the profile of the front is formally obtained from Eq. (56) as $h({\bf x},t)=(\nu/\alpha)\nabla^{-2}f$. Substituting this in Eq. (57) indicates that coupling to the height profile results in a non-local advection of the composition front through a term $(|\alpha\beta|/\nu)\nabla\left(f\nabla^{-1}f\right)$. The exponents in Table 1 thus correspond to those of such a field theory.

VII Summary

Using symmetry arguments, we have proposed a phenomenological model for the co-evolution of surface roughness and the evolutionary/ecological dynamics of bacterial range expansions through a coupling of the Directed Percolation and KPZ equations. In addition to the usual parameters of DP and KPZ, our model includes only four new parameters ($\alpha$, $\beta_{1}$, $\beta_{2}$, and $c$) which have well defined macroscopic manifestations. In particular, we have explored how these parameters modify the profile of a growing film and can in turn be computed through macroscopic measurements. We have further explored how the coupling of the DP composition to a rough front profile modifies the critical behavior at the DP transition (in turn modifying the scaling of KPZ fluctuations) through an RG analysis. In fact all couplings become relevant below the upper critical dimension of four, suggesting a perturbative analysis in $d=4-\epsilon$, which we carry out. Unfortunately (but perhaps not surprisingly) the RG analysis is not conclusive: the most promising candidate exhibits ‘weak dynamic scaling,’ with exponents that cannot be profitably compared to existing numerical estimates in $d=1$. We may hope that extensive numerical simulations (including in higher dimensions), along with further analysis of the coupled equations, can shed more light on this problem.

Through various mappings the KPZ equation can be related to the Burgers equation for a randomly stirred fluid and to meanderings of a directed path in a random medium. Through similar steps, we show that certain limits of our proposed equations can be mapped (by a generalized Cole-Hopf transformation) to genealogical paths following ancestry of growing/mutating species, and that this mapping also enables exact solution of the co-evolving composition and front profile in the noiseless limit. We expect that variants of these equations are relevant to a host of other problems, including range expansions of chiral George and Korolev or mutualistic Lavrentovich and Nelson (2014) bacteria, or even growth in a weakly compressible Plummer et al. or turbulant fluid flow Antonov et al. (2009).

Acknowledgements.

The authors grateful acknowledge conversations with Alexandre Solon, Sherry Chu, and David R. Nelson. JMH is supported by the Gordon and Betty Moore Foundation as a Physics of Living Systems Fellow through Grant No. GBMF4513. MK acknowledges support from the NSF through grant number DMR-1708280.

Appendix A Microscopic dynamics growth rules

A.1 Flat growth

Illustrations with flat growth are implemented using the Domany-Kinzel (DK) cellular automaton rules or equivalently the “stepping stone” model Korolev et al. (2010); Lavrentovich et al. (2013). Simulations are performed on a lattice $L$ sites wide and $T$ sites deep. Each lattice site on layer $n\in[1,T]$ and position $i\in[1,L]$ is labeled by a composition value $f_{n,i}$ taking values $-1$ or $1$ for the two bacterial species, pictured respectively as blue or gray, and [math] on empty sites.

Each new layer is determined sequentially based on the composition of the previous layer. In updating the $n$-th layer, the identity in each position $i$ is determined by two ‘neighbors’ in the layer below, $n-1$. To implement the interactions specified by the tilted square lattice depicted in Fig. 1, the layer-below neighbors for even layers are $i$ and $i+1$, whereas the layer-below neighbors for odd layers are $i-1$ and $i$. The composition value at each lattice site $f_{n,i}$ along layer $n$ is then determined probabilistically from layer $n-1$ according to the rules

[TABLE]

for $n$ even, with symmetric rules if $n$ is odd. Here, $\mu_{\mu}$ is a microscopic mutation probability from species $1\to-1$, and $\tau_{\mu}$ is a microscopic selection bias (or growth rate differential) favoring the reproduction of species 1.

A.2 Rough growth

Rough growth simulations are implemented with a rough variant of the DK rules, sometimes refereed to as the modified “brick wall” solid-on-solid model Drossel and Kardar (2000, 2003), or “single step” model Meakin et al. (1986); Plischke et al. (1987). As in the the stepping stone model, simulations are performed on a lattice $L$ sites wide and $T$ sites deep. Each lattice site on layer $n\in[1,T]$ and position $i\in[1,L]$ is labeled by a composition value $f_{n,i}$ taking values $-1$ or $1$ for the two bacterial species, pictured respectively as blue or gray, and [math] on empty sites. However, here the composition value of lattice sites are not updated layer by layer, but sites are updated randomly.

Updates proceed by first randomly choosing a position $i$ and then determining the last nonzero layer where the composition value has not been updated, say $n$. The addition is then accepted depending on whether the two ‘neighbors’ in the layer below are nonzero. Again, to implement the interactions specified by the tilted square lattice depicted in Fig. 1, the layer-below neighbors for even layers are $i$ and $i+1$, whereas the layer-below neighbors for odd layers are $i-1$ and $i$. If either of the neighbors in the layer below is zero then the addition is rejected. If they are both full, rejections are composition dependent with probabilities, say for an even layer,

[TABLE]

If the addition of a new site is accepted, the composition is determined in a manner akin to the flat growth DK rules, except allowing for the influence of same level neighbors weighted by a factor $\beta_{\mu}$. The algorithm proceeds by first calculating the weighted fraction of neighbors with composition value $1$ of a newly added site $(n,i)$:

[TABLE]

where $\Theta(x)$ is the heaviside step function taking value $1$ for $x>0$ and [math] otherwise. Then based on this value the composition at site $(n,i)$ is determined from the probability distribution

[TABLE]

which are well defined probabilities for small $\beta_{\mu}$, $\tau_{\mu}$ and $\mu_{\mu}$, which is our region of interest. Apart from the illustrations in Fig. 1, all numerical calculations are performed with $\tau_{\mu}=\mu_{\mu}=0$.

Updates are performed a fixed number of times. All simulations finish before the addition of any new sites in the last layer $T$.

Appendix B Renormalization group calculation

The RG program operates by tracking the effect on the parameters of repeatedly integrating out short wavelengths (high momenta). However, this is only analytically tractable for nonlinear theories perturbatively around a solvable Gaussian theory near a critical dimension. Here, the solvable theory is provided by the Gaussian part of the MSRJD action [Eq. (26)] with nonlinear terms in the action [Eq. (27)] treated as small corrections.

To facilitate our presentation, let us introduce the vector of fields $S=({\tilde{h}},{h},{\tilde{f}},{f})$. Furthermore, as we have already argued that $D_{h}$ is irrelevant under RG while $c$ does not renormalize at the lowest order, allowing us to ignore them. That leaves for the Gaussian action [Eq. (26)] in Fourier space [with for example $k=(\omega,q)$ and $h(x,t)=\int\frac{d^{d}xd\omega}{(2\pi)^{d}}\ e^{iqx-i\omega t}h(q,\omega)$]

[TABLE]

where we have introduced the shorthand $d^{\prime}=d+1$ and the correlation matrix

[TABLE]

whose elements determine the unperturbed (connected) correlation functions or cumulants, sometimes called the bare propagators, $\langle S_{m}(k)S_{n}(l)\rangle^{c}=(2\pi)^{d}{\mathcal{A}}_{mn}(k)\delta(k+l)$:

[TABLE]

Renormalization proceeds by considering how the behavior of these correlation functions change when one perturbatively includes the short wavelength effects of the nonlinear interactions in the action [Eq. (27)], which in Fourier space read

[TABLE]

A formally instructive way to organize this perturbation expansion is to consider the generating function for the cumulants of the distribution Täuber (2014)

[TABLE]

whose derivates with respect to the dummy variables $j$ and ${\tilde{j}}$ generate the cumulants. For example, the two-point correlation function is

[TABLE]

The perturbative analysis of the cumulants can be simplified by constructing an alternative set of functions known as vertex functions Täuber (2014). They are obtained, like the cumulants, from a generating function as follows. We first define the conjugate variables to $j$ and ${\tilde{j}}$

[TABLE]

and from them obtain the vertex generating function as the Legendre transform of the cumulant generating function

[TABLE]

Its derivatives at the origin define the vertex functions

[TABLE]

Their relationship to the cumulants can be obtained by successive derivatives. For example, the first derivatives are

[TABLE]

as is expected for a Legendre transform. Another important example is obtained from two derivatives

[TABLE]

demonstrating that the two point vertex is the matrix inverse of the cumulants (or matrix of propagators). Successive derivatives lead to higher-order vertices. Importantly, as we turn to the perturbative calculation, the corrections to the vertex functions are determined from the class of one-particle irreducible diagrams.

With this formalism in hand, we now calculate the corrections to the vertex functions due to the effects of the short wavelength behavior propagated through the nonlinearities. This proceeds by perturbatively expanding the vertex function in the nonlinearities and then integrating over a small annulus of short wavelengths, $\Lambda/b<q<\Lambda$, with $b\gtrsim 1$. We will see that this perturbation theory is consistent – that is all the nonlinearities are indeed small – only near the critical dimension $d=4-\epsilon$. In addition, the critical value of the control parameter near $d=4-\epsilon$ is also expected to be small, allowing us to simplify the presentation by noting that within the order of our perturbation expansion $\tau$ will not affect the renormalization of any other parameter, but itself. With this in mind, we consider the lowest order perturbative corrections to the vertex functions, which can be graphically represented using Feynman diagrams containing one loop.

The first vertex function we consider is the inverse of the height response function

[TABLE]

which has no corrections to one-loop order. Thus, $\nu$ does not renormalize. (Note that in the standard KPZ equation, there is renormalization of $\nu$ due to finite $D_{h}$, which is absent here.)

Next, we consider the vertex function corresponding to the inverse of the DP response function, whose corrections are listed in Fig. 6:

[TABLE]

where

[TABLE]

where integrals were performed by first expanding to first order in $\omega$ and $q$, since we are only interested in the long time, long wavelength behavior. The procedure for performing such integrals has been described previously, for example in the Appendix of Ref. Medina et al. (1989).

The next two-point vertex to renormalize comes from the coupling of the response field ${\tilde{h}}$ to the DP composition field $f$. Its perturbative expansion depicted in Fig. 7, corresponds to

[TABLE]

where the only correction surviving with the external momenta set to zero $q=\omega=0$ is

[TABLE]

There are four three-point vertices corresponding to each of the nonlinear parameters:

• •

We first address the KPZ nonlinearity $\lambda$ whose renormalization corresponds to the vertex function

[TABLE]

with no corrections at one-loop order, implying that $\lambda$ does not renormalize at this order.

• •

On the other hand, the DP nonlinearities both renormalize. The corrections to the noise vertex

[TABLE]

depicted in Fig. 8

correspond to the integrals evaluated at zero external momenta

[TABLE]

• •

The other DP nonlinearity’s renormalization is captured with the vertex

[TABLE]

where the Feynman diagrams for the corrections are in Fig. 9.

The corresponding integrals at zero external momenta are

[TABLE]

• •

Finally, the nonlinear coupling of the height field to the DP order parameter occurs through the two beta couplings, whose renormalization can be handled together through the vertex

[TABLE]

where the corrections, depicted in Fig. 10, are

[TABLE]

Next, we rescale space, time and the fields as $x\to bx$, $t\to b^{z}t$, $f\to b^{y_{f}}f$, $h\to b^{y_{h}}h$, ${\tilde{f}}\to b^{{\tilde{y}}_{f}}{\tilde{f}}$, and ${\tilde{h}}\to b^{{\tilde{y}}_{h}}{\tilde{h}}$. The corrections are then absorbed into the model parameters, so that the formal structure of the theory remains invariant. The new effective model parameters now have changed by a magnitude that depends on $b$. By letting $b=e^{l}$, with $l\ll 1$, and setting $\Lambda=1$ for convenience, the change (or flow) in these parameters is described by the collection of RG equations

[TABLE]

coupled with the algebraic constrains on the scaling exponents

[TABLE]

which come from demanding the terms proportional to $\omega$ in the propagators are scale invariant.

A close examination of the RG equations reveals that they only depend on particular collections of variables, the dimensionless effective parameters in Eq. (VI.4). The associated RG equations for these effective parameters are

[TABLE]