# Emerging predictable features of replicated biological invasion fronts

^{a}Laboratory of Ecohydrology, School of Architecture, Civil and Environmental Engineering, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland;^{b}Department of Aquatic Ecology, Eawag: Swiss Federal Institute of Aquatic Science and Technology, CH-8600 Dübendorf, Switzerland; and^{c}Dipartimento di Ingegneria Civile, Edile ed Ambientale, Università di Padova, I-35131 Padua, Italy

See allHide authors and affiliations

Contributed by Andrea Rinaldo, November 12, 2013 (sent for review August 10, 2013)

## Significance

Biological dispersal is a key driver of several fundamental processes in nature, crucially controlling the distribution of species and affecting their coexistence. Despite its relevance for important ecological processes, however, the subject suffers an acknowledged lack of experimentation, and current assessments point at inherent limitation to predictability even in the simplest ecological settings. We show, by combining replicated experimentation on the spread of the ciliate *Tetrahymena* sp. with a theoretical approach based on stochastic differential equations, that information on local unconstrained movement and reproduction of organisms (including demographic stochasticity) allows reliable prediction of both the propagation speed and range of variability of invasion fronts over multiple generations.

## Abstract

Biological dispersal shapes species’ distribution and affects their coexistence. The spread of organisms governs the dynamics of invasive species, the spread of pathogens, and the shifts in species ranges due to climate or environmental change. Despite its relevance for fundamental ecological processes, however, replicated experimentation on biological dispersal is lacking, and current assessments point at inherent limitations to predictability, even in the simplest ecological settings. In contrast, we show, by replicated experimentation on the spread of the ciliate *Tetrahymena* sp. in linear landscapes, that information on local unconstrained movement and reproduction allows us to predict reliably the existence and speed of traveling waves of invasion at the macroscopic scale. Furthermore, a theoretical approach introducing demographic stochasticity in the Fisher–Kolmogorov framework of reaction–diffusion processes captures the observed fluctuations in range expansions. Therefore, predictability of the key features of biological dispersal overcomes the inherent biological stochasticity. Our results establish a causal link from the short-term individual level to the long-term, broad-scale population patterns and may be generalized, possibly providing a general predictive framework for biological invasions in natural environments.

What is the source of variance in the spread rates of biological invasions? The search for processes that affect biological dispersal and sources of variability observed in ecological range expansions is fundamental to the study of invasive species dynamics (1⇓⇓⇓⇓⇓⇓⇓⇓–10), shifts in species ranges due to climate or environmental change (11⇓–13), and, in general, the spatial distribution of species (3, 14⇓–16). Dispersal is the key agent that brings favorable genotypes or highly competitive species into new ranges much faster than any other ecological or evolutionary process (1, 17). Understanding the potential and realized dispersal is thus key to ecology in general (18). When organisms’ spread occurs on the timescale of multiple generations, it is the byproduct of processes that take place at finer spatial and temporal scales that are the local movement and reproduction of individuals (5, 10). The main difficulty in causally understanding dispersal is thus to upscale processes that happen at the short-term individual level to long-term and broad-scale population patterns (5, 18⇓–20). Furthermore, the large fluctuations observed in range expansions have been claimed to reflect an intrinsic lack of predictability of the phenomenon (21). Whether the variability observed in nature or in experimental ensembles might be accounted for by systematic differences between landscapes or by demographic stochasticity affecting basic vital rates of the organisms involved is an open research question (10, 18, 21, 22).

Modeling of biological dispersal established the theoretical framework of reaction–diffusion processes (1⇓–3, 23⇓–25), which now finds common application in dispersal ecology (5, 14, 22, 26⇓⇓⇓–30) and in other fields (17, 23, 25, 31⇓⇓⇓⇓–36). Reaction–diffusion models have also been applied to model human colonization processes (31), such as the Neolithic transition in Europe (25, 37, 38). The classical prediction of reaction–diffusion models (1, 2, 24, 25) is the propagation of an invading wavefront traveling undeformed at a constant speed (Fig. 1*E*). Such models have been widely adopted by ecologists to describe the spread of organisms in a variety of comparative studies (5, 10, 26) and to control the dynamics of invasive species (3, 4, 6). The extensive use of these models and the good fit to observational data favored their common endorsement as a paradigm for biological dispersal (6). However, current assessments (21) point at inherent limitations to the predictability of the phenomenon, due to its intrinsic stochasticity. Therefore, single realizations of a dispersal event (as those addressed in comparative studies) might deviate significantly from the mean of the process, making replicated experimentation necessary to allow hypothesis testing, identification of causal relationships, and to potentially falsify the models’ assumptions (39).

Here, we provide replicated and controlled experimental support to the theory of reaction–diffusion processes for modeling biological dispersal (23⇓–25) in a generalized context that reproduces the observed fluctuations. Firstly, we experimentally substantiate the Fisher–Kolmogorov prediction (1, 2) on the existence and the mean speed of traveling wavefronts by measuring the individual components of the process. Secondly, we manipulate the inclusion of demographic stochasticity in the model to reproduce the observed variability in range expansions. We move from the Fisher–Kolmogorov equation (*Materials and Methods*) to describe the spread of organisms in a linear landscape (1, 2, 24, 25). The equation couples a logistic term describing the reproduction of individuals with growth rate *r* and carrying capacity *K* and a diffusion term accounting for local movement, epitomized by the diffusion coefficient *D* . These species’ traits define the characteristic scales of the dispersal process. In this framework, a population initially located at one end of a linear landscape is predicted to form a wavefront of colonization invading empty space at a constant speed (1, 2, 24, 25), which we measured in our dispersal experiment (Fig. 1*D* and *SI Text*).

## Results

In the experiments, we used the freshwater ciliate *Tetrahymena* sp. (*Materials and Methods*) because of its short generation time (16) and its history as a model system in ecology (16, 40, 41). The experimental setup consisted of linear landscapes, filled with a nutrient medium, kept in constant environmental conditions and of suitable size to meet the assumptions about the relevant dispersal timescales (*Materials and Methods*). Replicated dispersal events were conducted by introducing an ensemble of individuals at one end of the landscape and measuring density profiles throughout the system at different times, through image analysis (*Materials and Methods*). Density profiles are shown in Fig. 2, in six replicated dispersal events (Fig. 2 *A*–*F*). Organisms introduced at one end of the landscape rapidly formed an advancing front that propagated at a remarkably constant speed (Fig. 3 and Table S1). The front position at each time was calculated as the first occurrence, starting from the end of the landscape, of a fixed value of the density (Fig. 2). As for traveling waves predicted by the Fisher–Kolmogorov equation, the mean front speed in our experiment is notably constant for different choices of the reference density value (Fig. 3*C*).

The species’ traits *r*, *K*, and *D* were measured in independent experiments (Table 1). In the local-growth experiment, a low-density population of *Tetrahymena* sp. was introduced evenly across the landscape, and its density was measured locally at different times. Recorded density measurements were fitted to the logistic growth model, which gave the estimates for *r* and *K* (Table 1 and Table S2). In the local unimpeded movement experiment, we computed the mean-square displacement (*SI Text* and Fig. S1) of individuals’ trajectories (42⇓–44) to estimate the diffusion coefficient *D* in density-independent conditions (Table 1 and *Materials and Methods*). The growth and movement measurements were performed in the same linear landscape settings as in the dispersal experiment and therefore are assumed to accurately describe the dynamics at the front of the traveling wave in the dispersal events.

The comparison of the predicted front speed to the wavefront speed measured in the dispersal experiment, , yields a compelling agreement. The observed speed in the dispersal experiment was (mean ± SE) (Table S1), which we compare with the predicted one cm/d (mean ± SE). The two velocities are compatible within one SE. A *t* test between the replicated observed speeds and bootstrap estimates of gives a *P* value of (, ). Thus, the null hypothesis that the mean difference is 0 is not rejected at the 5% level, and there is no indication that the two means are different. As the measurements of *r* and *D* were performed in independent experiments, at scales that were orders of magnitude smaller than in the dispersal events, the agreement between the two estimates of the front velocity is deemed remarkable.

Although the Fisher–Kolmogorov equation correctly predicts the mean speed of the experimentally observed invading wavefront, its deterministic formulation prevents it from reproducing the variability that is inherent to biological dispersal (21). In particular, it cannot reproduce the fluctuations in range expansion between different replicates of our dispersal experiment (Fig. 3*A*). We propose a generalization of the Fisher–Kolmogorov equation (*Materials and Methods*) accounting for demographic stochasticity that is able to capture the observed variability. The strength of demographic stochasticity is embedded in an additional species’ trait *σ* . In this stochastic framework, the demographic parameters *r*, *K*, and σ were estimated from the local growth experiment with a maximum-likelihood approach (Table 1 and *Materials and Methods*) whereas the estimate of the diffusion coefficient *D* was left unchanged (*SI Text*). We then used these local independent estimates to numerically integrate the generalized model equation (45, 46), with initial conditions as in the dispersal experiment, and found that the measured front positions are in accordance with simulations (Fig. 3*A* and Fig. S2). In particular, most experimental data are within the 95% confidence interval for the simulated front position, and the observed range variability is well-captured by our stochastic model (Fig. 3*B*). Accordingly, the estimate for the front speed and its variability in the experiment are in good agreement with simulations (*SI Text*). Demographic stochasticity can therefore explain the observed variability in range expansions.

## Discussion

To summarize, we suggest that measuring and suitably interpreting local processes allows us to accurately predict the main features of global invasions. The deterministic Fisher–Kolmogorov equation is shown to correctly predict the mean speed of invasion but cannot capture the observed variability. Instead, characterizing the inherent stochasticity of the biological processes involved allows us to predict both the mean and the variability of range expansions, which is of interest for practical purposes, such as the delineation of worst-case scenarios for the spread of invasive species. Our phenomenological approach allows us to make predictions on the spread of organisms without the need to introduce all details on the movement behavior, biology, or any other information. Such details are synthesized in three parameters describing the density-independent yet stochastic behavior of individuals riding the invasion wave. The parsimony of the model allows generalization to organisms with different biology (e.g., growth rates and diffusion coefficients are available for several species in the literature) (6) and supports the view that our protocol may possibly provide a general predictive framework for biological invasions in natural environments.

In conclusion, we have shown that, at least in the simple ecological settings investigated here, predictability remains, notwithstanding biological fluctuations, owing to the stochastic treatment devised. We confirm that deterministic models can be applied to describe ecological processes and show that additional information on the stochasticity acting at the mesoscopic scale allows us to estimate fluctuations at the macroscopic scale. We believe that our results might have implications for the dynamics of phenomena other than species’ invasions, such as morphogenesis (23, 47), tumor growth (23, 25, 36), and the spreading of epidemics (23, 30, 34, 35), which have been traditionally modeled with reaction–diffusion equations.

## Materials and Methods

### Study Species.

The species used in this study is *Tetrahymena* sp. (Fig. 1*B*), a freshwater ciliate, purchased from Carolina Biological Supply. Individuals of *Tetrahymena* sp. have typical linear size (equivalent diameter) of 14 μm (41). Freshwater bacteria of the species *Serratia fonticola*, *Breviacillus brevis*, and *Bacillus subtilis* were used as a food resource for ciliates, which were kept in a medium made of sterilized spring water and protozoan pellets (Carolina Biological Supply) at a density of 0.45 gL^{−1}. The experimental units were kept under constant fluorescent light for the whole duration of the study, at a constant temperature of 22 °C. Overall, experimental protocols are well-established (16, 41, 48⇓–50), and the contribution of laboratory experiments on protists to the understanding of population and metapopulation dynamics proved noteworthy (48).

### Experimental Setup.

Experiments were performed in linear landscapes (Fig. 1*A*) filled with a nutrient medium and bacteria of the three species above mentioned. The linear landscapes were 2 m long, 5 mm wide, and 3 mm deep, respectively, and 10^{5}, 350, and 200 times the size of the study organism (41). Landscapes consisted of channels drilled on a Plexiglas sheet. A second sheet was used as lid, and a gasket was introduced to avoid water spillage (Fig. 1*A*). At one end of the landscapes, an opening was placed for the introduction of ciliates. The Plexiglas sheets were sterilized with a 70% (vol/vol) alcohol solution, and gaskets were autoclaved at 120 °C before filling the landscape with medium. As Plexiglas is transparent, the experimental units could be placed under the objective of a stereomicroscope, to record pictures (for counting of individuals) or videos (to track ciliates). Individuals were observed to distribute mainly at the bottom of the landscape, whose length was three orders of magnitude larger than its width (*w*) and depth and two orders of magnitude larger than the typical length scale of the process . The population was thus assumed to be confidently well-mixed within the cross section after a time , which in our case is of the order of a minute after introduction of the ciliates in the landscape.

### Experimental Protocol.

We performed three independent and complementing experiments, specifically: (*i*) a dispersal experiment was carried out to study the possible existence and the propagation of traveling invasion wavefronts in replicated dispersal events; (*ii*) a growth experiment was run to obtain estimates of the demographic species’ traits, which are *r* and *K* in the deterministic framework of Eq. **1** and *r*, *K*, and σ in the stochastic framework of Eq. **2**; (*iii*) a local movement experiment was performed to study the local unimpeded movement of *Tetrahymena* sp. over a short timescale (in a time window ), to estimate the diffusion coefficient *D* for our study species, independently from the dispersal and growth experiments.

#### Dispersal experiment.

We performed six replicated dispersal events in the linear landscapes. After filling the landscapes with medium and bacteria, a small ensemble of *Tetrahymena* sp. was introduced at the origin. Subsequently, the density of *Tetrahymena* sp. was measured at 1-cm intervals, five times in the first 48 h and twice in the last 48 h. The whole experiment lasted for about 20 generations of the study species.

#### Local growth experiment.

We performed five replicated growth measurements in the linear landscapes, to measure the demographic species’ traits, in the same environmental conditions as in the dispersal experiment, but independently from it. A low-density culture of *Tetrahymena* sp. was introduced in the whole landscape, and its density was measured by taking pictures and counting individuals, covering a region of 7 cm along the landscape. Density measurements were performed at several time points for each of the five replicates, in a time window of 3 d.

#### Local movement experiment.

We performed four additional, replicated dispersal events in the linear landscapes, initialized in the same way as in the dispersal experiment, to measure the diffusion coefficient of *Tetrahymena* sp. The diffusion coefficient *D* is the proportionality constant that links the mean square displacement of organisms’ trajectories to time (42, 44) (*SI Text*). Macroscopically, it relates the local flux to the density of individuals, under the assumption of steady state (44). To estimate the diffusion coefficient, we recorded several videos of individuals moving at the front of the traveling wave (at low density), reconstructed their trajectories (42, 43), and computed their mean square displacement .

##### Video recording.

We recorded videos of *Tetrahymena* sp. at the front of the traveling wave in four replicated dispersal events, at various times over 4 d. The area covered in each video was of 24 mm in the direction of the landscape and 5 mm orthogonal to it. Each video lasted for 12 min.

##### Trajectories reconstruction.

For each recorded video, we extracted individuals’ spatial coordinates in each frame and used the MOSAIC plugin for the software ImageJ to reconstruct trajectories (43). The goodness of the tracking was checked on several trajectories by direct comparison with the videos. Examples of reconstructed trajectories can be seen in Fig. 1*C* or in Movie S1.

##### Diffusion coefficient estimate.

For each video, the square displacement of each trajectory in the direction parallel to the landscape was computed at all time points and then averaged across trajectories. Precisely, for each trajectory *i* we computed the quantity , where is the 1-dimensional coordinate of organism *i* at time *t* in the direction parallel to the landscape and is its initial position. The mean square displacement in a video was then computed as the mean of across all trajectories, that is, (where *N* is the total number of trajectories). A typical measurement of is shown in Fig. S1. As shown in the figure, there exists an initial correlated phase, which we discuss in *SI Text*. To estimate the diffusion coefficient from the mean square displacement, we fitted the measured to the function (*SI Text*) with the two parameters *D* (diffusion coefficient) and τ (correlation time). The total number of recorded videos was 28, that is, 7 for each replica.

### Mathematical Models.

#### Deterministic framework.

The Fisher–Kolmogorov equation (1, 2). reads:

where is the density of organisms, *r* the species’ growth rate, *D* the diffusion coefficient, and *K* the carrying capacity. Eq. **1** is known to foster the development of undeformed traveling waves of the density profile. Mathematically, the existence of traveling wave solutions implies that , where *v* is the speed of the advancing wave. Fisher (1) proved that traveling wave solutions can only exist with speed , and Kolmogorov (2) demonstrated that, with suitable initial conditions, the speed of the wavefront is the lower bound.

The microscopic movement underlying the Fisher–Kolmogorov Eq. **1** is brownian motion (25, 51). Investigation of the movement behavior of *Tetrahymena* sp., instead, shows that individuals’ trajectories are consistent with a persistent random walk with an autocorrelation time . The corresponding macroscopic equation for the persistent random walk should thus be the reaction–telegraph equation (25) (*SI Text*). Nonetheless, as the autocorrelation time for our study species is much smaller than the growth rate *r* , Eq. **1** provides an excellent approximation to the reaction–telegraph equation. See *SI Text* for a detailed discussion.

#### Stochastic framework.

The stochastic model equation reads:

where is a Gaussian, zero-mean white noise (i.e., with correlations , where δ is the Dirac’s delta distribution) and is constant. We adopt the Itô’s stochastic calculus (51), as appropriate in this case. Note, in fact, that the choice of the Stratonovich framework would make no sense here, as the noise term would have a constant nonzero mean (22, 51), which would allow an extinct population to possibly escape the zero-density absorbing state. The square-root multiplicative noise term in Eq. **2** is commonly interpreted as describing demographic stochasticity in a population (46) and needs extra care in simulations (45, 52). In particular, standard stochastic integration schemes fail to preserve the positivity of ρ. We adopted a recently developed split-step method (45) to numerically integrate Eq. **2**. This method allows us to perform the integration with relatively large spatial and temporal steps maintaining numerical accuracy.

Data from the growth experiment were fitted to the equation:

where is the local density, is a Gaussian, zero-mean white noise (i.e., with correlations ), is constant, and *l* is the size of the region over which densities were measured (*SI Text*). Eq. **3** describes the time-evolution of the density in a well-mixed patch of length *l* (*SI Text*). The likelihood function for Eq. **3** can be written as:

where *n* is the total number of observations in the growth time series, is the vector of demographic parameters, and is the transitional probability density of having a density of individuals ρ at time *t*, given that the density at time was (for a given θ). The transitional probability density satisfies the Fokker–Planck equation associated to Eq. **3** (*SI Text*), which was solved numerically for all observed transitions and choices of parameters (*SI Text*), adopting the implicit Crank–Nicholson scheme. The best fit parameters were those that maximized the likelihood function Eq. **4** (Table 1 and *SI Text*).

## Acknowledgments

We thank Amos Maritan and Enrico Bertuzzo for insightful discussions, Janick Cardinale for support with the MOSAIC software, and Regula Illi for the ciliates’ picture. We acknowledge the support provided by the discretionary funds of Eawag: Swiss Federal Institute of Aquatic Science and Technology; by the European Research Council advanced grant program through the project “River networks as ecological corridors for biodiversity, populations and waterborne disease” (RINEC-227612); and by Swiss National Science Foundation Projects 200021_124930/1 and 31003A_135622.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: andrea.rinaldo{at}epfl.ch.

Author contributions: A.G., A.R., and F.A. designed research; A.G. and F.C. performed the experiments; A.G. analyzed data; A.R., F.C., and F.A. assisted with data analysis; A.G., A.R., and F.A. wrote the paper; and A.G. developed the stochastic theoretical framework.

The authors declare no conflict of interest.

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

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- Kolmogorov AN,
- Petrovskii IG,
- Piskunov NS

- ↵
- Skellam JG

- ↵
- Elton CS

- ↵
- ↵
- ↵
- Shigesada N,
- Kawasaki K

- ↵
- Clobert J,
- Danchin E,
- Dhondt AA,
- Nichols JD

- ↵
- Okubo A,
- Levin SA

- ↵
- ↵
- ↵
- Bell G,
- Gonzalez A

- ↵
- ↵
- Bertuzzo E,
- Maritan A,
- Gatto M,
- Rodriguez-Iturbe I,
- Rinaldo A

- ↵
- ↵
- Carrara F,
- Altermatt F,
- Rodriguez-Iturbe I,
- Rinaldo A

- ↵
- ↵
- ↵
- ↵
- ↵
- Melbourne BA,
- Hastings A

- ↵
- ↵
- Murray JD

- ↵
- ↵
- Méndez V,
- Fedotov S,
- Horsthemke W

- ↵
- ↵
- ↵
- Solé RV,
- Bascompte J

- ↵
- Roques L,
- Garnier J,
- Hamel F,
- Klein EK

- ↵
- Pybus OG,
- et al.

- ↵
- ↵
- Hallatschek O,
- Hersen P,
- Ramanathan S,
- Nelson DR

- ↵
- ↵
- Bertuzzo E,
- Casagrandi R,
- Gatto M,
- Rodriguez-Iturbe I,
- Rinaldo A

- ↵
- Righetto L,
- et al.

- ↵
- ↵
- Ammerman AJ,
- Cavalli-Sforza LL

- ↵
- Fort J

- ↵
- ↵
- ↵
- Giometto A,
- Altermatt F,
- Carrara F,
- Maritan A,
- Rinaldo A

- ↵
- Berg HC

- ↵
- ↵
- Polin M,
- Tuval I,
- Drescher K,
- Gollub JP,
- Goldstein RE

*Chlamydomonas*swims with two “gears” in a eukaryotic version of run-and-tumble locomotion. Science 325(5939):487–490. - ↵
- ↵
- Bonachela J,
- Muñoz MA,
- Levin SA

- ↵
- Turing AM

- ↵
- ↵
- ↵
- ↵
- Gardiner C

- ↵