Population genetics in compressible flows
Abstract
We study competition between two biological species advected by a compressible velocity field. Individuals are treated as discrete Lagrangian particles that reproduce or die in a densitydependent fashion. In the absence of a velocity field and fitness advantage, number fluctuations lead to a coarsening dynamics typical of the stochastic Fisher equation. We then study three examples of compressible advecting fields: a shell model of turbulence, a sinusoidal velocity field and a linear velocity sink. In all cases, advection leads to a striking drop in the fixation time, as well as a large reduction in the global carrying capacity. Despite localization on convergence zones, one species goes extinct much more rapidly than in wellmixed populations. For a weak harmonic potential, one finds a bimodal distribution of fixation times. The longlived states in this case are demixed configurations with a single boundary, whose location depends on the fitness advantage.
pacs:
87.23.Cc, 47.27.EChallenging problems arise when spatial migrations of species are combined with population genetics. The population dynamics of a single species expanding into new territory was first studied in the pioneering works of Fisher, Kolmogorov, Petrovsky and Piscounov (FKPP) (1); (2); (3). Later, Kimura and Weiss studied individualbased counterparts of the FKPP equation (4), revealing the important role of number fluctuations. In particular, stochasticity is inevitable at a frontier, where the population size is small and the discrete nature of the individuals becomes essential. Depending on the parameter values, fluctuations can produce radical changes with respect to the deterministic predictions (3); (5). If is the population fraction of, say, a mutant species and that of the wild type, the stochastic FKPP equation reads in one dimension (6):
(1) 
where is the spatial diffusion constant, is the genetic diffusion constant (inversely proportional to the local population size ), is the genetic advantage of the mutant and is a Gaussian noise, deltacorrelated in time and space that must be interpreted using Ito calculus (6); (7).
However, many species, from the distant past (8) up to the present, have competed in liquid environments, such as lakes, rivers and oceans. Interesting new phenomena arise when population dynamics is coupled to hydrodynamic flows (9). For example, satellite observations of chlorophyll concentrations have identified long lived, segregated patches of marine microorganisms off the eastern coast of the southern tip of South America (10), where species domains are largely determined by the tangential velocity field obtained from satellite altimetry. In cases such as these, the dynamics takes place in the presence of advecting flows, some of them at high Reynolds numbers (11). It is often appropriate to consider a compressible velocity field, both because of inertial effects (12) associated with fairly large microorganisms (diameter 5500m), and because photosynthetic bacteria and plankton often control their buoyancy to stay close to the ocean surface (13). In the latter case, the coarsegrained velocity field advecting the microorganisms will contain a compressible component to account for downwellings (14). Recent works studied the dynamics of a single population in the presence of a compressible turbulent velocity field in one and two dimensions (15); (16). Here, the interplay between turbulent dynamics and population growth leads to âquasilocationâ on convergence zones and a remarkable reduction in the carrying capacity.
The study of competition in a hydrodynamics context, where both a compressible velocity field and stochasticity due to finite population sizes are present, calls for a nontrivial generalization of Eq. (1). One complication is that, because of compressibility, the sum of the concentrations of the two species is no longer invariant during the dynamics. Thus, we must clarify the definition of , the fraction of one particular species. A biologically important issue arises from overshooting: the density can exceed unity near velocity sinks, resulting in an unphysical imaginary noise amplitude in Eq. (1).
In this Letter, we overcome these problems by introducing a new model, designed to explore, for the first time, how compressible velocity flows affect the competition between two different species which can die and reproduce via cell division. In absence of an advecting field, the model reproduces the coarsening dynamics of spatial population genetics, as shown for one spatial dimension in Fig. 1a. Of particular importance in population genetics is the fixation time , defined as the average time needed for one of the two species to reach extinction. The value of is proportional to the total population size in mean field approximation (i.e. assuming uniform spatial mixing) and grows to for a one dimensional system (6). Conversely, Fig. 1b shows a typical simulation with our model using a high Reynold number synthetic velocity field extracted from a shell model of turbulence (17). In striking contrast to the result without advecting flows (Fig. 1a and Ref. (6)), one now finds competitions carried out on transient but persistent sinks in the velocity field. We find that turbulence not only compresses and reduces the overall population density, it also greatly reduces the fixation time . Despite quasilocalization onto velocity sinks, the typical fixation time is much shorter than if particles were well mixed. We have also considered a compressible sinusoidal velocity field and a confining sink arising from a linear velocity profile. Remarkably, we find that these low Reynolds number compressible flows can also dramatically reduce fixation times. In the following, we investigate this phenomenon and propose a phenomenological explanation.
Eq. (1) is derived by assuming a fixed total concentration of individuals, so that if is the concentration of one species, the other will have a concentration exactly equal to . To describe cases in which the total concentration can change due to an advecting flow, we introduce an offlattice model in which two different organisms, and , advect and diffuse in space, while undergoing duplication (i.e. cell division) and densitydependent annihilation (death). Specifically, we implement the following reactions: each particle of species duplicates with rate and annihilates with a rate , where is the number of neighboring particles (of both types) in an interaction range , and is the total number of organisms that can be accomodated in the unit interval with total density . To reduce the number of parameters, we set , but take to allow the possibility of a selective advantage (faster reproduction rate) of species . We will start by analyzing in depth the neutral case and consider the effect of in the end of the Letter. In one dimension, our macroscopic coupled equations for the densities and of individuals of type and in an advecting field read
(2) 
with and . As discussed in (18), for an interval of size , let , where is the local carrying capacity of a region of size . and are independent deltacorrelated noise sources with an Itocalculus interpretation as in Eq. (1), and we have rescaled the particle densities so that the local carrying capacity is when . Although may seem nongeneric in the context of dynamical systems theory, this is a case of particular interest in population genetics, where it corresponds to the neutral theory of Kimura (4).
The above equations follow from a microscopic master equation via the KramersMoyal method (18), at this order equivalent to Van Kampen inverse system size expansion (19). Multispecies versions of the FKPP equations have been already considered in Refs (20), but not in the presence of an advecting field or number fluctuations.
We first simulated the particle model directly for , finding a dynamics similar to that of the stochastic FKPP equation, Eq. (1), for , see Fig.(1a). In this simple limit, our model can be considered as a grandcanonical generalization of Eq (1), where the total density of individuals is now allowed to fluctuate around an average value . Details of the numerical implementation and the role of density fluctuations are presented in (18). Hereafter, we fix the following parameters: , , and where is the one dimensional domain endowed with periodic boundary conditions. With these parameters, would be for the one dimensional FKKP equation, and for the wellmixed case.
Introducing a compressible velocity field , as shown in Fig.(1b), leads to radically different dynamics. Individuals tend to concentrate at sinks in the velocity field. Further, competition is enhanced and the total number of individuals present at time is on average significantly smaller than the total carrying capacity . We implemented three different velocity fields: 1) a velocity field generated by a shell model of compressible turbulence (17); (15), reproducing the power spectrum of a high Reynolds number turbulent velocity field of forcing intensity (18), 2) a static sine wave, , and 3) a linear velocity profile, , confining particles near the origin. Periodic boundary conditions on the unit interval are used in cases (1) and (2), while in case (3) the system size is so large that particles never reach the boundaries.
Fig.(2) shows the average fixation time for in the first two cases, at varying the intensity of advection. In the left panel, we plot the fixation times as a function of the time averaged reduced carrying capacity , where is the carrying capacity reduction, i.e. the ratio between the actual number of particles and the average number of particles observed in absence of the velocity field. Plotting vs. allows comparisons with the mean field prediction, , valid for well mixed systems (black dashed line) (6). For the shell model, we include simulations of the macroscopic equations (Population genetics in compressible flows) with different resolutions ( and lattice sites on the unit interval), obtaining always similar results for vs. .
In all cases, the presence of a spatially varying velocity field leads to a dramatic reduction of , compared to the mean field theory. The fixation time drops abruptly as soon as , even for very small . Simulations suggest a singular limit of zero intensity () of the velocity field (and consequently in Fig. 2), as we discuss later.
These observations can be understood by assuming that global fixation is determined by coalescence of allele boundaries like those shown in Fig. (1a). Although in the stochastic FKPP with the boundaries perform a random walk (6), here they are also advected by the velocity field. In particular, the average position of one boundary , neglecting diffusion and number fluctuations, will simply evolve according . The fixation time is then determined by the time needed for the boundaries to reach the center of a sink and annihilate. Because periodic boundary conditions are imposed, the number of interfaces is necessarily even, so all of them annihilate pairwise. Upon expanding the velocity field around a sink position , , we find a fixation time . Hence, we model the average time to fixation vs. the velocity field intensity as =, where represents the typical time to reach fixation at large strain rate and is an dimensionless constant of order unity. This phenomenological theory describes well the data of Fig. 2, right panel (dotdashed lines), where in the turbulent case we take the forcing as a proxy for ( for the sine wave and for turbulence).
We now study particles subject to a converging linear velocity, . When is large (and is small), the fixation time is comparable to the mean field prediction. However, as the forcing decreases, it becomes much larger than the mean field prediction (Fig. 3a, dashed line). In this regime, the fixation time probability distribution is bimodal (Fig. 3a, inset): roughly half of the realizations have a short fixation time (faster than mean field), while the other half maintain coexistence for an extended period. The spacetime evolution of these configurations in Fig. 3b reveals that the two mutant species are demixed, one on the left and one on the right of the sink. Correspondingly, the average total heterozigosity , defined as the probability of two random individuals to be of different type (6), decays exponentially for very large (consistent with mean field theory (6)), but tends to a constant, nonzero value when smaller values of are chosen (Fig 3c).
This remarkable behavior arises because of the different boundary conditions for the linear sink, such that the number of domain boundaries is no longer always even. Hence, the initial number of boundaries will be odd approximately half the time for random initial conditions, leading to the single surviving boundary as shown on the right of Fig. 3a. This configuration corresponds to a stable stationary solution of the deterministic version of Eqs Population genetics in compressible flows with . We expect that this demixed solution, with and nonzero on opposite sides of the sink, should have a lifetime (inaccessible in numerical simulations) which grows exponentially with (7); (19).
Finally, what happens for a selective advantage ? Fig. 4 shows a comparison between the neutral case of Fig. 3 and a simulation in which red particles reproduce faster (). Even with such a large selective advantage, genetic interfaces are still present. In this case, we estimate the shift in the interface as , by equating the outward Fisher genetic wave velocity (1) with the inward advection velocity at distance from the origin.
To conclude, we have introduced a new model which allows studies of how compressible advection affects Darwinian competition. A compressible flow leads to a radically new scenario, with the fixation time being largely determined by how species boundaries collapse into the sinks of the velocity field, rather than by diffusive annihilation. Our phenomenological theory suggests a significative effect even for very weak compressible fields. Compressible advection becomes irrelevant for the fixation time only when the time to drift into the sinks (of order ) is much larger than the diffusion time (of order , where is a typical linear size of the population). An obvious application of the model is the study of compressible flows in higher space dimensions, including flows relevant to biological oceanography of photosynthetic organisms near the water surface. Recent advances in experimental techniques could also lead to tests of the one dimensional results. Motility of bacteria has been studied in microstructures of diameter comparable or even smaller to that of the organisms (21). In this setting, a compressible flow could be created by pumping liquid nutrient from the two ends of the tube and extracting it from the center via a semipermeable membrane. Another possibility could be to study floating microorganisms (13); (22) and make use of the compressible effects caused by the vertical component of a convecting velocity field (14).
Acknowledgements.
We are grateful to F. Toschi and P. Perlekar for interesting discussions. Work by SP and MHJ was supported by the Danish National Research Foundation through the Center for Models of Life. Support for DRN was provided by the National Science Foundation in part through Grant No. DMR1005289 and by the Harvard Materials Research Science and Engineering Center through NSF Grant No. DMR0820484.I Supplementary Materials
In this note we present the technical derivations of the results in the Letter “Population Genetics in compressible flows” (“main text”). First, we describes the derivation of the macroscopic equations for the particle model introduced in the main text. Then, we describe details of the continuum simulations and of the shell model of turbulence. Finally, we discuss the relation between the model introduced in the main text and the simpler stochastic FisherKolmogorovPetrovskyPiskounov (FKPP) equation (Eq. (1) of the main text).
Ii Macroscopic equations for the two species model
To derive the mean field version of the model discussed in the Letter. we consider first positionindependent quantites and which represent the number of particles of species and respectively. This “zerodimensional” limit corresponds to a single well mixed “deme” or island in a stepping stone model. The birth and death processes are characterized by the following transition rates:
(3) 
where is the probability per unit time of population to increase its size by one unit when the total populations are and , and similarly for and .
There exist several methods to derive macroscopic equations for particle models (7); (23); (19), most popular being the Van Kampen system size expansion and the KramersMoyal expansion. While the former is more rigorous when studying problems beyond the linear noise approximation, the latter is normally used because it leads to more transparent expressions. Moreover, it can be shown (see e.g. (7), page 251) that at the order of the Fokker Planck equation the two methods are completely equivalent. The KramersMoyal method works as follows: consider a master equation, and call the transition rate from particles to particles. When is typically large and relatively small ( in the case of birthdeath process considered here) a Taylor series expansion of the master equation in for the probability of having particles at time is appropriate:
(4) 
with
(5) 
Truncating the Taylor expansion at the second order in leads to a FokkerPlanck equation.
It is straightforward to apply the above procedure to two species and and obtain the zerodimensional FokkerPlanck equation for the joint probability distribution function . Substituting the expression of the transition rates of Eq. (II) into (4 5) leads to:
(6) 
corresponding to the stochastic differential equations:
(7) 
where and are deltacorrelated Gaussian processes in time with unit amplitude, e.g. and are uncorrelated with each other, . For consistency we have to adopt the Ito prescription for interpreting the nonlinear noise term (7).
From Eq. (II) we see that typical values of near the steady state are of order . Upon defining the rescaled particle densities
(8) 
we obtain
(9) 
These coupled Langevin equations approximate well the master equation when is large. If is not large, the noise cannot be approximated by a Gaussian noise.
The model is generalized to include advection and diffusion in space as follows: Particles belonging to the two species advect and diffuse in space in a Lagrangian way, i.e. the position of a given particle evolves according to:
(10) 
where is a unit amplitude Brownian noise source that implements diffusion with diffusion constant . At each time step , a given particle duplicates with probability or annihilates due to competition with rate , where are now the number of particles of type in a region of size centered in the particle itself.
The macroscopic equations in the spatial case can be derived using arguments similar to the well mixed case. We assume that binary reactions can occur only when both particles are within a length , and also take this as the size of our discretization volume. Although the total particle number in this region fluctuates, plays a role similar to an island spacing in a stepping stone model. We first write the discretized stochastic differential equation for the and , the number of particles in the cell of size :
(11) 
where the noise terms now obey .
Here, we assume the ration is adjusted to match the overall carrying capacity of the cell, and that is big enough to insure the local carrying capacity . When is sufficiently large, the diffusion and advection terms enter only in the deterministic part of our equations.
There are two possible choices for passing to the concentration equations:

Change variables to the integrated dimensionless concentration inside each cell: .

Change variables to a normalized real concentration (number of particles per unit volume): in dimensions (we focus here on ).
In the first case, one is left with a discretized set of well defined Langevin equations, a natural generalization of the discrete stepping stone model with dimensionless species concentrations. The first procedure, after setting as before , results in:
(12) 
which is a set of well defined stochastic differential equations, valid for a particular lattice spacing .
The second procedure, however, leads to a more convenient passage to the continuum limit,
(13) 
Starting with Eq.(II), a formal way to achieve the continuous limit is to set first by rescaling the densities. This choice corresponds to an average number of particles equal to in a cell in the absence of fluid flow when the total density is constant and equal to . Then, we take the continuum limit at fixed , leading to:
(14)  
(15) 
where now , and . We also defined and . We have chosen units of length such that . To return to dimensional units, simply let in and . To understand the connection with stepping stone models, it is helpful to note that , so that the amplitude of the noise terms in Eqs. (14) and (15) can be written in a similar form as the “generic diffusion constant” discussed in Ref. (6). Eqs. (14) and (15) correspond to the macroscopic equations presented in the Letter in absence of a selective advantage, . The derivation of the nonneutral is a straightforward generalization in which the parameter controls the relative difference of the duplication rate of species with respect to that of species . Mathematically, this correspondly to repeat the expansion but substituting to the first equation in (II) the more general expression:
(16) 
leading to the the expression of Eq. (2) in the main text.
The continuous limit of such reactiondiffusionadvection processes is only a concise summary of the stochastic dynamics; underlying this is always the spacediscretized version (see Ref. (7), page 314). One reason is that advection and diffusion terms can contribute significantly to the variance of the noise when the cells are small ((7), page 313); these terms can be neglected only if the equations will be discretized on a sufficiently coarse lattice. Moreover, the validity of the KramersMoyal (or Van Kampen) expansion relies on . Even in absence of a carrying capacity reduction due to compressible flows, the number of particles in a sufficiently small cell may not be particularly large. The macroscopic continuum equations will then be consistent with results from particle simulations only when there is a sufficiently large number of particles when we coarsegrain to scales small compared to the typical scales of the density gradients. This is the case for the parameter ranges considered in the main text.
Iii Continuous simulations and shell models
In addition to particle simulations, Eq.s (1415) (i.e. Eqs. (2) in the main text) have been simulated numerically by using finite difference methods in space and the EulerCauchy method for time steps with a variety of velocity fields . The effect of number fluctuations (i.e. genetic drift) has been taken into account by using the method proposed by Doering et. al. (25): at each time step we compute the positiondependent quantities . If (i.e. if there is a nonzero number of particles of the appropriate type in the region of size ), then we add the noise term; otherwise we set the noise term equal to . It has been shown by Doering et. al. that this numerical scheme converges, in a suitable norm, to the continuum limit. In our numerical simulations, we found excellent agreement between the results obtained by integrating eq.s (1415) with this technique, compared to those obtained by particle methods, as will be discussed in details in a forthcoming paper.
Next we illustrate how we built the ”turbulent” velocity field used in our high Reynolds number simulations. Although we focus primarily on one spatial dimension, our analysis of the statistical properties of employs velocity fluctuations typical of three dimensional turbulent fluid mechanics. The statistical properties of our advecting velocity field will be characterized by intermittency both in space and in time. To obtain a velocity field with generic space and time dynamics at high Reynolds number, we build up by appealing to a simplified shell model of fluid turbulence (24). Wavenumber space is divided into shells centered around , . For each shell with characteristic wavenumber , we describe turbulence by using the complex Fourierlike variable , satisfying the following equations of motion:
(17)  
The model contains one free parameter, , and it conserves two quadratic invariants (when the force and the dissipation terms are absent) for all values of . The first is the total energy and the second is , where . We fix , as this choice reproduces intermittency features of the real three dimensional Navier Stokes equation with surprising good accuracy (24). Using , we can build the real one dimensional velocity field as
(18) 
where the constant controls the strength of velocity fluctuations relative to other parameters in the model. In all numerical simulations we use a forcing function , i.e. energy is supplied only to the largest scale corresponding to , where is the injection rate of turbulent energy. With this choice, the power injected into the shell model is simply given by , i.e. it is constant in time. To solve Eqs. (1415) and (17) we use a finite difference scheme with periodic boundary conditions in space.
Iv Stochastic FKPP limit
In this section, we clarify the connection between the advectionfree model described in the Letter and in Eqs. (14,15) and the stochastic FKPP (FisherKolmogorovPetrovskyPiscounov) equation. To handle advecting flows, our model generalizates the FKPP equation by relaxing the constraint on the total density, for all values of and , see Fig. (5).
To understand this correspondence, it is helpful to change variables from the densities , to the relative densities where . Eqs. (1415) now become:
(19) 
where now the noise amplitudes are given by and by definition .
Upon substituting and in absence of advection (), each of the fractions evolves according to the stochastic FKPP equations for two neutral alleles.
A further issue is to understand if there is a limit in which the term can be neglected, so that one retrieves Eq. (IV) without having to impose any constraint. An estimate follows from studying the fluctuations of , which obeys a closed equation:
(20) 
where . By defining and assuming , we obtain a linear stochastic differential equation,
(21) 
By means of a Fourier transform we can diagonalize the diffusion operator and find that the Fourier modes are independent Gaussian variables with zero mean and variance . It follows that
(22) 
The fluctuations in will be negligible when is large. In addition, it is interesting to study
(23) 
where we introduced a large cutoff (related to the discretization size, ) in the above integral to avoid an ultraviolet divergence. With this premise, the above integral is small if is small, suggesting that the logarithmic terms in Eqs. (IV) can be neglected in this limit.
References
 R. A. Fisher, Ann. Eugenics 7, 353369 (1937).
 A. Kolmogorov, N. Petrovsky, and N. Piscounov, Moscow Univ. Math. Bull. 1, 125 (1937).
 W. van Saarloos, Phys. Rep. 386, 29222 (2003).
 M. Kimura and G. H. Weiss, Genetics 49, 561576 (1964); J. F. Crow and M. Kimura An Introduction to Population Genetics, Blackburn Press, Caldwell, NJ (2009).
 O. Hallatschek and K. Korolev, Phys. Rev. Lett. 103, 108103 (2009), and references therein.
 For a recent review, see K. Korolev et al. Rev. Mod. Phys. 820, 16911718 (2010).
 C. Gardiner Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences, Springer Series in Synergetics (2004).
 B. A. Whitton and M. Potts The Ecology of Cyanobacteria: Their Diversity in Time and Space eds. Kluwer, Dordrecht, Netherlands (2000).
 T .Tel et al., Physics Reports 413, 91â196 (2005).
 F. D’Ovidio et al., Proc. Natl. Acad. Sci. 107, 1836618370 (2010).
 W. J. McKiver and Z. Neufeld, Phys. Rev. E 79, 061902 18 (2009).
 J. Bec, Phys. Fluids 15, L81L84 (2003).
 A. P. Martin, Progr. Ocean. 57, 125174 (2003).
 J. R. Cressman et al., Europhys. Lett. 66, 219225 (2004); G. Boffetta et al., Phys. Rev. Lett. 93, 134501 (2004).
 R. Benzi, D. R. Nelson Physica D 238 20032015 (2009).
 P. Perlekar et al., Phys. Rev. Lett. 105, 144501 (2010).
 M. H. Jensen et al., Phys. Rev. A 43, 798805 (1991) ; T. Bohr et al., Dynamical Systems Approach to Turbulence Cambridge University Press, Cambridge (1998).
 Supplementary Material are available online.
 N. G. van Kampen Stochastic Processes in Physics and Chemistry, Elsevier, Amsterdam (2007).
 M. O. Vlad et al., Phys. Rev. E 65, 061110 (2002); M. O. Vlad et al., Proc. Natl. Acad. Sci. 13, 1024910253 (2004).
 J. Männik et al. Proc. Natl. Acad. Sci. 35, 1486114866 (2009).
 P. B. Raney and M. Travisano Nature 294, 6972 (1998).
 Risken H (1996) The FokkerPlanck Equation: Methods of Solution and Applications, SpringerVerlag, Berlin.
 Biferale L (2003) Annu. Rev. Fluid Mech. 35, 441.
 Doering CR, Sargsyan KV, Smereka P (2005) Physics Letters A 344: 149D155; see also Doering CR, Mueller C and Smereka P (2003) Physica A 325: 243259.