## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# The noisy edge of traveling waves

Edited by* Pierre C. Hohenberg, New York University, New York, NY, and approved November 30, 2010 (received for review September 12, 2010)

### Related Articles

- Leading the dog of selection by its mutational nose- Feb 02, 2011

## Abstract

Traveling waves are ubiquitous in nature and control the speed of many important dynamical processes, including chemical reactions, epidemic outbreaks, and biological evolution. Despite their fundamental role in complex systems, traveling waves remain elusive because they are often dominated by rare fluctuations in the wave tip, which have defied any rigorous analysis so far. Here, we show that by adjusting nonlinear model details, noisy traveling waves can be solved exactly. The moment equations of these tuned models are closed and have a simple analytical structure resembling the deterministic approximation supplemented by a nonlocal cutoff term. The peculiar form of the cutoff shapes the noisy edge of traveling waves and is critical for the correct prediction of the wave speed and its fluctuations. Our approach is illustrated and benchmarked using the example of fitness waves arising in simple models of microbial evolution, which are highly sensitive to number fluctuations. We demonstrate explicitly how these models can be tuned to account for finite population sizes and determine how quickly populations adapt as a function of population size and mutation rates. More generally, our method is shown to apply to a broad class of models, in which number fluctuations are generated by branching processes. Because of this versatility, the method of model tuning may serve as a promising route toward unraveling universal properties of complex discrete particle systems.

- nonequilibrium statistical mechanics
- population genetics
- solitary waves
- speed of adaptation
- stochastic processes

The wave-like spread of discrete entities pervades our everyday life. For example, the spread of ions, pathogens, and beneficial mutations control the human heart beat, the yearly threat of influenza, and evolutionary progress (1). A thorough understanding of how these waves form and spread has numerous applications ranging from the control of chemical reactions (2, 3) to the prediction of epidemic outbreaks (4, 5). Great research effort has therefore been made on the question of how traveling waves emerge in complex systems from the multiplicity of relatively simple interactions, in particular the random dispersal of particles and reactions between particles. To simplify the analysis, most theoretical studies have been neglecting number fluctuations, which are inevitable in systems of discrete particles. However, when stochastic simulations became feasible, it was found that those previously ignored fluctuations can have a strong impact on the dynamics of waves (6, 7, 8). Simple models of biological evolution (described below) serve as the prime example for this drastic sensitivity on noise because they are dominated by the few most fit individuals in a population, and they break down if number fluctuations are neglected (6). With the added difficulty of particle discreteness, the analysis of traveling waves became one of the important challenges of statistical physics and mathematical biology, which still defies systematic analytical techniques.

Here, we show that an exact analysis of noisy traveling waves is feasible when the nonlinear details of the underlying model are chosen appropriately. The main idea is to *tune* the nonlinearities to obtain the least difficult math, while retaining the universal features of the model. Specifically, we show that it is possible to close the hierarchy of moment equations using a suitably designed nonlinear constraint on the dynamics. The tuned model can then be used to extract the universal features of noisy traveling waves, for which only heuristic approaches have been available so far. The method of model tuning can be naturally generalized to a wide range of unsolved stochastic nonlinear problems, and therefore provides a promising tool to unraveling universal features of nonequilibrium systems.

## Noisy Wave Models

Simple models of evolution (6, 9–13) provide spectacular examples of noisy traveling waves because of their drastic sensitivity to rare fluctuations in the wave tip. As illustrated in Fig. 1, these waves describe the continual increase of growth rate, also called fitness, due to spontaneous mutations in a finite population. The wave speed, which is a measure for how quickly populations adapt, increases without bound if number fluctuations are neglected (6). To reproduce a *steady* state with constant wave speed observed in stochastic simulations, any analysis has to incorporate number fluctuations, at least heuristically by introducing an ad hoc “cutoff” in the noisy tip of the wave (6). The relationship between traveling speed and population size has been investigated extensively in the recent literature (13–16), with some controversy regarding the universal behavior (17, 18), in order to interpret growth rate measurements in microbial evolution experiments (19, 20).

We now take a closer look at the models underlying these fitness waves to elucidate their characteristic structure. This will lead us to a general model of noisy traveling waves, which forms the basis of our analytic approach. Most evolution models implement the reproduction of individuals by a standard branching process. The growth rate is subject to small variations due to random mutations, which are modeled through a random walk, for instance, by standard diffusion (6). In addition, all evolution models contain a nonlinearity that limits the total population size. This accounts for the fact that natural populations must stay finite due to resource limitations.

The generic mathematical structure of such evolution models can be framed as follows: The state of the population at time *t* is described by a function *c*_{t}(*x*) that represents the number density of individuals having a growth rate *x*. The population is assumed to evolve in discrete time steps of size *ϵ*, which is eventually sent to zero in order to obtain a continuous time Markov process. Any such time step consists of two substeps. The first one, illustrated in Fig. 2*A*, embodies the mutational process, which leads to slight variations in growth rates, and the process of reproduction, which causes individuals with growth rate above the mean to increase in relative number. In general, the combined effect of both processes can be described by the equation [1]which takes the concentration field from *c*_{t} to an intermediate value . The term represents the *expected* change in density due to reproduction and mutations. This term is linear in the concentration field because the number of offspring and mutants per time step is proportional to the current population density. The linear operator is similar to a Hamiltonian in physics. A specific example will be discussed later on. The stochastic term in Eq. **1** accounts for all random factors that influence the reproduction process. The function *η*(*x*) represents standard white noise, i.e., a field of delta correlated random numbers, , where denotes the ensemble average of a random variable *f*. The amplitude of the noise term in Eq. **1** is typical for number fluctuations: Because of the law of large number, the expected variance from one time step to the next is proportional to the number *ϵc*_{t} of expected births or deaths during one time step.

Because the reproduction step changes particle numbers, another substep of population control is required to enforce a constant population size. In most models and experiments (19, 20), this step is realized by a random culling of the population: Individuals are eliminated at random from the population until the population size constraint is restored; see Fig. 2*B*. Mathematically, the selection step can be cast into the form [2]where *λ* represents the fraction of the population that has to be removed to comply with the population size constraint. The second substep completes the computational time step and takes the concentration field from the intermediate state to the properly constrained state *c*_{t+ϵ}.

The combination of noise and nonlinearity (via the population constraint) has led many theoretical studies to engage in approximations that are often hard to justify or to control^{†}. Here, we take a different approach by optimizing the *form* of the nonlinearities in order to obtain an exactly solvable model. To this end, we first generalize the above model in the following way: Instead of enforcing a fixed population size, we allow for a whole class of constraints. Specifically, we assume that the selection step enforces a constant value of 1 for the inner product of the concentration field *c*_{t}(*x*) and a new function *u*(*x*), [3]Here, we have introduced the “bra-ket” notation commonly used in quantum mechanics. The function *u*(*x*) defines the selection rule and will be called *selection function* from now on. It will serve as a “tuning wheel” in the following. If one chooses *u* = *N*^{-1}, one obviously recovers the fixed population size constraint. For any other choice, the population size will not be fixed. At best, one obtains a steady state with a population size fluctuating around its mean value, .

## Results

### Model Tuning.

Eqs. **1**–**3** define a class of models that generate branching random walks subject to a global constraint. As we argue in *Discussion*, these features not only characterize fitness waves but furnish the essence of most noisy traveling waves. The fundamental difficulty in analyzing such models becomes apparent when we try to determine the typical dynamics by averaging over the stochastic noise term. This can be done by eliminating *λ* using the constraint, Eq. **3**, sending *ϵ* to zero and carrying out the ensemble average. This straightforward calculation, detailed in *SI Appendix*, yields [4]for the mean concentration field . Eq. **4** reflects the usual “horror” inherent to nonlinear stochastic problems: Moment equations do not close in general. The mean depends on the second moment, the second on the third, and so on. A whole hierarchy of moment equations would have to be solved self-consistently to make progress.

However, if we consider *u*(*x*) as a tunable function, Eq. **4** allows for different kind of simplification. Suppose, the solution *u*_{∗}(*x*) of [5]exists, and we choose *u*_{∗}(*x*) as the selection function. For this particular model, the nonlinearity in Eq. **4** disappears identically. Thus, the dynamics of the first moment becomes linear, [6]and hence solvable. In particular, if a steady state exists, the steady-state concentration field is in the null space of . These observations suggest an algorithm to construct for a given linear operator , a solvable constrained branching walk model: First, identify the selection function *u*_{∗}(*x*) for which the averaged dynamics becomes simple (i.e., linear). To this end, one has to solve Eq. **5**, which is in general a deterministic partial differential equation. Second, solve the corresponding moment Eq. **6**, which is guaranteed to be linear for the chosen selection function *u*_{∗}(*x*).

At this point it is not clear whether Eqs. **5** and **6** have stationary solutions, and, if so, whether these solutions can be used to extract universal features of traveling waves. Therefore, we apply our recipe to a simple model of fitness waves. This model is defined by the operator [7]which corresponds to a branching random walk on a reaction rate gradient (22). It has been proposed as the simplest model of asexual evolution in order to explain the fitness growth observed in directed evolution experiments with large populations of RNA viruses (6). The term (*x* - *x*_{0})*c*_{t} in Eq. **7** represents the reproduction with *x*_{0} being the mean growth rate of the population. The mutational process is modeled as a diffusion process with diffusivity *D*. This diffusion approximation is justified when mutations occur at a higher rate than the typical growth rate difference they confer, which is an appropriate assumption for viruses with large mutation rates. For microbes with lower mutation rates, it is necessary to discretize the fitness landscape and model mutations as fitness jumps of finite size. The discretized scenario is slightly more complex as it has one additional parameter, namely, the characteristic fitness effect of mutations, but can be discussed in complete analogy to the present case of a continuous fitness landscape (*SI Appendix*). Finally, the term *v*∂_{x}*c*_{t} appears in Eq. **7** because the dynamics of the fitness wave is described in a reference frame moving with the average speed *v* of the fitness wave.

Using in Eq. **5**, we find that the equation [8]specifies the selection function *u*_{∗}(*x*). Once *u*_{∗} is determined, the mean concentration field follows as the stationary solution of the linear Eq. **6**. This second task requires no special effort if the random walks are simple diffusion processes as in Eq. **7**. Then, and *u*_{∗} are simply related by [9]which generates Eq. **8** when inserted into Eq. **6**. Thus, follows from *u*_{∗} directly by multiplying with a decaying exponential. The prefactor in **9** is set by the constraint Eq. **3**.

The set of Eqs. **8** and **9** is controlled by just one parameter *vD*^{-2/3}, as can be seen by introducing scaled variables. In Fig. 3, we depict representative numerical solutions for both regimes of large and small values of this control parameter and compare with simulations of the tuned model. We find that both regimes exhibit strongly different wave profiles. Whereas the mean concentration exhibits for the most part a Gaussian shape in regime of large wave speeds (*vD*^{-2/3} > 1), it is strongly skewed for small speeds. Notice the nearly perfect agreement between numerics and stochastic simulations, which confirms our theoretical framework and demonstrates the feasibility of our model tuning approach.

### Interpretation of the Tuned Models.

Among the three graphs depicted in Fig. 3, only the function has an obvious interpretation as the fitness wave profile. Using the results obtained in ref. 23, one can also give an intuitive interpretation of the functions *u*_{∗}(*x*) and . Both functions relate to the phenomenon of fixation. Imagine sampling an individual at position *x* and labeling it with an inheritable label (neutral mutation). As the dynamics proceeds, the abundance of this label will change due to number fluctuations and the fitness of its carriers. Eventually, this label will either go extinct or become *fixed* in the population. The latter case occurs if the descendants of the labeled individual take over the population.

Fixation events are much more likely if the initially labeled individual belongs to the fitter part of the population. We thus expect the probability of fixation to be a steeply increasing function of *x*, similar to the function *u*_{∗}(*x*). Indeed, the interpretation of *u*_{∗}(*x*) is precisely that of a fixation probability of a particle at position *x*, which is derived in *SI Appendix*. It is interesting to note about Eq. **5** determining *u*_{∗}(*x*), that a very similar equation is well-known to describe the survival probability of an *unconstrained* branching random walk (12). The only difference is the prefactor of 2 instead of 1 in the nonlinearity of Eq. **5**.

The product also has an intuitive interpretation in terms of a probability density^{‡}. It represents the positional distribution function of the individual whose descendants eventually will take over the population (23). Even though there must surely be such a “lucky” individual at any time, it cannot be described deterministically, simply because the fixation event depends on random future events. Thus, the position of the “common ancestor of future generations” can be described only probabilistically, similar to the position of quantum mechanical particle.

### Wave Speed.

Next, we compare the evolution model with its tuned counterpart, considering at first the relation between speed and population size. From Fig. 4, it can be seen that deviations are significant only for intermediate values of the control parameter *vD*^{-2/3} and disappear in the limit of large and small values. The scaling *v* ∼ ln^{1/3}*N* for large speeds has been observed previously, but could be modeled only by introducing a cutoff into the mean-field equations (6, 22). Our tuned model, which exhibits the same scaling, can be used to rationalize and refine the heuristic cutoff idea. To this end, consider the closed equation for the mean population density, [10]which is obtained by eliminating *u*_{∗} in Eq. **8** using **9** and the constraint (Eq. **3**). The closed moment Eq. **10** has the form of the deterministic approximation except for the last term, which acts (for large population densities) similar to a cutoff: It is negligible for small *x* and completely dominates for sufficiently large *x* because of the exponential “amplification” factor in the numerator. This leads to an asymptotic wave profile ∼*e*^{-vx} in the tip of the wave, just as found from a heuristic cutoff in the reaction rate (7). The location of the cutoff can be determined by balancing the last term in Eq. **10** with the reaction term . In Fig. 3, this crossover point can be identified as the x value where *u* crosses over from exponential to linear. Notice that the functional form of the cutoff is quite different from a step function, which has been used as a cutoff in ref. 6 and most subsequent studies. The peculiar nonlinear form of the cutoff turns out to be essential for determining the leading order corrections to the wave speed (and its fluctuations) in the limit of large speeds, as we show in the *SI Appendix*. In the opposite limit of small speeds *v* → 0, we find that the tuned model exhibits *v* ∼ *N*/2, which can be derived from a perturbation analysis of the moment equation (24). Surprisingly, the exact same asymptotic is measured for the original evolution model, suggesting universal behavior not only for large but also for small speeds.

### Fluctuations.

Another telling comparison concerns the fluctuations of waves, rather than the mean density field. At first sight, fluctuations in both models seem to be very different because they arise from different statistical ensembles: Whereas wave speeds fluctuate in the evolution model by Tsimring et al. (6), they are perfectly constant in the tuned model. On the other hand, population sizes are constant in the evolution model but fluctuate in the tuned model. One might expect, however, that both ensembles become equivalent in the universal large population size limit, similar to different ensembles in the thermodynamic limit of statistical mechanics. On the basis of this equivalence hypothesis, one can try to infer the fluctuations in wave speed of the evolution model (and thus wave diffusivities) based on the fluctuations of the mortality rate *λ* in the solvable tuned model (*SI Appendix*). This reproduces correctly the known diffusivity scaling *D*_{wave} ∼ ln^{-3}*N* of Fisher–Kolmogorov waves (25) and predicts a unique scaling *D*_{wave} ∼ ln^{-1}*N* for evolutionary waves. We stress, however, that these analytic results rest on the equivalence of the two descriptions in the large population size limit, which remains to be proven.

## Discussion

As a case study for noisy traveling waves in general, our analysis concentrated on simple, yet unsolved, models of evolution. These models effectively describe branching random walks subject to a global constraint (12): The branching random walk has the effect of modeling the growth of the population and the mutations through diffusion along the fitness axis. The global constraint fixes the population size accounting for the fact that natural population must stay finite due to resource limitations.

Even these simplest evolution models could not be solved previously because of the nonlinearity associated with the rigid constraint. However, although a limit on the population size is natural, there is no obvious biological reason for demanding a nonfluctuating population size. After abandoning the idealized fixed population size constraint, we could show that a solvable evolution model with a closed moment equation can be obtained when the form of the constrained is tailored to satisfy Eq. **5**. The solvable model was found to reproduce the properties of the fixed population size model very well, in particular, for large and small population sizes (or wave speeds).

Our approach of model tuning not only applies to evolutionary waves but to any solitary wave that is generated by a branching random walk under a global constraint, or some other nonlinearity, to keep particle numbers finite. Such models, which we may term constrained branching random walks, share universal features that depend only on the form of the linear part of the model, as has been found in extensive computer simulations (8). The universality class is defined by the linear part of the model, which describes a branching process and a random walk. For instance, the important class of Fisher–Kolmogorov waves (26, 27) has a reaction term that saturates at a constant value in the tip of the wave, in contrast to the linear reaction rate gradient that characterizes the above evolutionary waves. This difference results in a different scaling of the wave speed and its dependence on number fluctuations, which is correctly reproduced by our approach beyond the leading order cutoff approximation (7) (*SI Appendix*). Wave speeds are also strongly dependent on the stochastic particle motion, which is diffusive only in the simplest cases. The spread of epidemic waves in the modern human population, for instance, is dominated by the scale-free air transportation network (5). This leads to superdiffusion (28), which can be incorporated in our framework by the use of a fractional diffusion operator in . Apart from these obvious extensions of our approach, our method can also be applied to problems with time-dependent scenarios (*SI Appendix*). This allows to account for dynamic driving forces and to study, for instance, waves in temporally changing environments.

In summary, we have developed and applied an exact method to determine the universal dynamics of noisy traveling waves. By tuning nonlinear details of traveling wave models without changing their universal features, we could close the hierarchy of moment equations for a large class of noisy traveling waves. These solvable models are amenable to a simple probabilistic interpretation in terms of fixation processes. For large population sizes (or weak noise), our description resembles mean-field theory supplemented with a cutoff and thus provides a strong rationale for the heuristic cutoff approach (6, 7). The peculiar, nonlocal form of the nonlinearity sets the precise location of the cutoff, which is crucial for the correct prediction of the wave speed and its fluctuations. More generally, the analysis presented in this article is applicable to those reaction–diffusion problems in which a global nonlinear interaction suffices to keep particle numbers finite. This encompasses simple models of microbial evolution and range expansions (Fisher–Kolmogorov waves), for which we explicitly demonstrated the utility of our approach. It will be interesting to see whether model tuning can be extended to situations where interactions are inherently local, such as in the problem of directed percolation (29).

## Acknowledgments

We are grateful to many useful discussions with Daniel S. Fisher. We thank Richard Neher and Kirill Korolev for stimulating discussions and comments on the manuscript.

## Footnotes

^{1}E-mail: oskar.hallatschek{at}ds.mpg.de.Author contributions: O.H. designed research, performed research, analyzed data, and wrote the paper.

The author declares no conflict of interest.

*This Direct Submission article had a prearranged editor.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1013529108/-/DCSupplemental.

↵

^{†}A remarkable exception is the exactly solvable “exponential” model by Brunet et al. (21)↵

^{‡}Note that, as required for a probability density, the integral of over*x*is equal to 1 by virtue of the constraint (Eq.**3**).

Freely available online through the PNAS open access option.

## References

- ↵
- Murray JD

- ↵
- Winfree A

- ↵
- Kuramoto Y

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Rouzine IM,
- Wakeley J,
- Coffin JM

- ↵
- Park S,
- Krug J

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Fisher R

- ↵
- Kolmogorov A,
- Petrovsky I,
- Piskunov N

- ↵
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Physical Sciences

- Biological Sciences
- Evolution