Particle algorithms for population dynamics in flows

Particle algorithms for population dynamics in flows This article has been downloaded from IOPscience. Please scroll down to see the full text article. 2011 J. Phys.: Conf. Ser. 333 012013 (http://iopscience.iop.org/1742-6596/333/1/012013) Download details: IP Address: 147.83.119.146 The article was downloaded on 30/08/2012 at 11:15 Please note that terms and conditions apply. View the table of contents for this issue, or go to the journal homepage for more Home Search Collections Journals About Contact us My IOPscience Particle algorithms for population dynamics in flows Prasad Perlekar1, Roberto Benzi2, Simone Pigolotti3, and Federico Toschi1 1Department of Physics and Department of Mathematics and Computer Science, Eindhoven University of Technology, Eindhoven 5600MB, The Netherlands. CNR-IAC, Via dei Taurini 19, 00185 Rome, Italy. 2Dipartimento di Fisica and INFN, Universit`a “Tor Vergata”, Via della Ricerca Scientifica 1, I-00133 Roma, Italy. 3Departament de Fisica i Enginyeria Nuclear, Universitat Polit`ecnica de Catalunya Edif. GAIA, Rambla Sant Nebridi s/n, 08222 Terrassa, Barcelona, Spain. E-mail: [email protected] Abstract. We present and discuss particle based algorithms to numerically study the dynamics of population subjected to an advecting flow condition. We discuss few possible variants of the algorithms and compare them in a model compressible flow. A comparison against appropriate versions of the continuum stochastic Fisher equation (sFKPP) is also presented and discussed. The algorithms can be used to study populations genetics in fluid environments. 1. Introduction Bacterial colonies growing on the surface of a Petri dish form intriguing patterns that can be mathematically modeled and studied in terms of a stepping stone model of population genetics [1]. Oceans cover a vast amount of the Earth’s surface supporting an incredibly diversity of species, playing a key role since early stage of life on our planet. Many microorganisms must have found ways to survive in the flows generated in ocean. Therefore, though not much explored, it is important to understand the role of flows on models for the dynamic of populations. To understand this question two approaches can be taken: (a) Use continuum models and couple them with fluid flows [3, 4]. This is an approach used commonly in oceanography models to study plankton blooms [2]; (b) Use particle based models coupled to flows. The advantage of this approach is that it allows us to pose in a simpler way questions related to the populations to number fluctuations and the extinction of a mutant species, which in turn is central to the understanding of population genetics [1, 10]. A similar approach has been employed earlier to study the role of discrete effects in the propagation of a Fisher wave [5]. In this paper we discuss algorithms to study particle-based population dynamics in flows, following the above approach (b). An entity based algorithm with different variants is discussed in section 2. We then validate our algorithm for the Moran process in section 3, where we first present the results for the case of uniformly mixed populations with zero diffusivity and compare our results against theoretical predictions. Next we focus on the spatial diffusivity. The continuum limit of Moran model with spatial diffusivity is the stochastic Fisher-KolmogorovPetrovsky-Piscounov equation (sFKPP) [7, 8, 11]. We study the role of population fluctuations on the speed of a Fisher wave and compare our results against the asymptotic predictions of Ref. [11]. Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 Published under licence by IOP Publishing Ltd 1 Finally, in section 4 we study the reaction A→A+A,A+A→Ain a flow, which also has the Fisher equation as a continuum limit, but with a different noise term [11, 15]. We show that in presence of a simple, deterministic, compressible flow the average population size (carrying capacity) decreases, in agreement with earlier investigations [4, 15]. 2. Numerical algorithm One of the example of population dynamics is the growth of bacteria on a Petri dish. The motion of the frontier of a bacterial colony can be understood as a solution of the stochastic Fisher equation. Another way to study frontier growth is by encoding the birth and death rules at an individual entity level. The advantage of this approach is that biologically relevant questions about the birth and extinction of species can be easily addressed. We now describe the two algorithms used by us to study population genetics with flows. 2.1. Algorithm 1 The numerical algorithm that is used to solve the population dynamics is similar to the algorithm used to study the front propagation of Fisher equation in Ref. [5]. We consider each individual in the population as a particle that is advected by the flow, diffused by the Brownian motion, and that can reproduce itself, die, or compete with other particles. Birth, death and competition can be modeled in terms of binary reactions. In order to efficiently implement these binary reactions (which necessarily happen when the two particles sufficiently close to each other) we introduce a spatial mesh. More precisely, we divide our one dimensional domain of size Linto Msubintervals of width δ=L/M and where δis the interaction distance. As a typical starting condition one may consider particles uniformly distributed at positions x∈L. The population is evolved according to the following two main steps: •Advection and spatial diffusion: the particles are advected and diffused according to xi(t+ ∆t) = xi(t) + u∆t+√2D∆tΓi(t), i= 1,2,...n. Here Dis the spatial diffusion coefficient, uis the advecting velocity, ∆tis the time step, nis the total number of particles, and Γ(t) is a Gaussian white noise with zero mean and unit variance. •Population dynamics: the particle population is assumed homogeneous inside the subdomain δand interact with each other using a prescribed set of rules determining the population dynamics. Given the framework that models advection, diffusion and locality of interaction, one must specify the proper reactions which encode for the biological processes in the populations. Here below we describe how to implement three of the building blocks at the base of the different reaction models used in this manuscript. Supposing that our ecosystem is composed by two distinct populations, Aand B, we denote the number of particles within the interaction radius, δ, of type Aas NAand those of type Bas NB. The total number of particles in the interval δ is indicated with N=NA+NBwhereas the total number of particles in Lis indicated with n. (i) A+Bk −→ B+B: In this type of reaction, a particle of one type, A, interacts with a particle of the other type, B, and at a rate kconverts itself into B. To implement this rule we first count the number of Bparticles present in δ. If kNB∆t < rtthen the particle Aconverts itself into Botherwise nothing happens (rtdenotes a random number uniformly distributed in [0,1]). (ii) Ak −→ A+A: In this reaction a particle of type Agives birth to another particle of the same type at a rate k. If k∆t < rtthen the particle multiplies into two; otherwise nothing happens. Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 2 (iii) A+Ak −→ A: In this reaction a particle of type Ainteracts with another particle of its own kind and one of the two dies with a rate k. Again we first count the number of Aparticles present in δ, if k(NA−1)∆t < rtthen the particle annihilates; otherwise nothing happens. 2.2. Algorithm 2 This type of algorithm constitutes a variant over algorithm 1. Here the advection and diffusion parts are implemented in the same way. Two main differences are implemented in the population dynamics part of the algorithm: Grid size At each timestep, the particle distribution is binned on a finer resolution. The bin size is chosen as δ/m, where δis the interaction distance as before and mis an odd natural number (in the following we will choose m= 11). Binary reactions occur at a rate which depends on the total number of neighbors in the row of mcells centered on the bin where each particle is. Rates The timestep ∆tis taken small enough that no more than one reaction of each type can occur within each timestep. Then, a reproduction A→2Aoccurs with a probability kn(t)∆twhere n(t) is the total number of particles. A death by competition, A+A→A occurs with a rate k∆tPjNj(t), where Nj(t) is the total number of particles in the row centered in the bin where particle jis and the sum runs on all the particles present in the system. In the case of reproduction, the reproducing particle is chosen at random with uniform probability, while in the case of death by competition the choice is weighted with the number of neighbors, i.e. particle ihas a probability to be chosen equal to Ni/PjNj. The first modification allow to have competition in an area δmore precisely centered in the particle position. We found that the second modification slows down simulation considerably for our parameter values. Therefore, it should be implemented only in cases where the timescales of population dynamics are much longer than those of advection diffusion, so that one is not forced to choose a very small time step in order to avoid having more than one reaction per time step. 3. Validation The validation is divided in two parts aimed at testing, respectively, the reaction algorithm alone and the diffusion-reaction implementation. We first focus on the reaction part and validate our population dynamics algorithm for the case of Moran process [6, 9] of population genetics. Here we compare our results with the theoretical prediction on the fixation probability from Kimura [9] (see Subsection 3.1). In subsection 3.2 we add the effect of the spatial diffusion: with this addition the Moran process becomes equivalent to a stochastic Fisher equation [7, 8, 11]. We compare the results of our discrete simulations against the analytical predictions based on the stochastic Fisher equation, both in strong and weak noise limit [11, 14]. 3.1. Moran model We consider a homogeneous population of two species Aand Bwith NAindividuals of type-A and NBindividuals of type Bin a domain of size δ. When an individual of type A(B) interacts with an entity of the other type B(A) then B(A) converts into A(B) with a rate kA(kB). These reaction processes can be written in form the following form: (A+BkA −−→ A+A A+BkB −−→ B+B During the above processes the total number of entities N≡NA+NBremain conserved for all times. The continuum limit of the above process is the stochastic equation dc = Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 3 µc(1 −c) + pσ2c(1 −c)dW (FKPP equation with noise) where c=NA/N is the concentration of the type-Aspecies, µ≡(kA−kB)N/δ and σ2≡(kA+kB) δ,dW is a Wiener process. One of the central question in population genetics is the fixation probability of a mutation. For what concerns the Moran model this question can be formulated as: if Bis the parent species and Ais the mutant, what is the probability that the mutant converts all the parent species B into itself A(i.e. Agets fixated)? The two control parameters are the initial mutant population concentration p=cA(t= 0) = NA/N and the selective advantage µ. For convenience, in this section we set δ= 1. For the case of the Moran process Kimura [9, 10] predicted that the fixation probability is: Pfix =1−e−αNp 1−e−αN (1) where α≡kA−kB kA+kBand pis the initial concentration of the mutant population. Note that for kA=kBor µ= 0 (case corresponding to no selective advantage), the fixation probability Pfix =p. The plot in Figure 1 shows the fixation probability for two different populations of size Nwith p= 1/N and varying µ. The negative values of the selective advantage implies that the parent species has a selective advantage over the mutant. We find that, in agreement with the theoretical prediction, the fixation probability increases with the selective advantage of the mutant. For the special case of zero selective advantage (kA=kB) we find, in agreement with theory, Pfix =p. In our simulations the fixation probability, Pfix , is defined as one minus the fraction of realizations for which the mutant dies. Our measurements were obtained by averaging over an ensemble of 2000 realizations for each data point. 3.2. Stochastic Fisher equation To study the effect of spatial diffusivity, Fisher, Kolmogorov et al. [7, 8] added a diffusion term to the Fisher equation: dc = [µc(1 −c) + D∇2c]dt +pσ2c(1 −c)dW. (2) In absence of noise, for the case of a localized an initial condition, the Fisher equation exhibits a traveling wave solution with a front speed given by vF= 2√Dµ [7, 8]. More recently, it was shown that the speed of the Fisher wave is reduced in presence of noise. In the regime of weak noise [11, 12, 13, 14]: v∼pDµ 2−π2 (log N)2(3) where Nis the total population within the interaction radius δand σ2∝1/N. In the strong noise regime the Fisher speed is drastically reduced [11, 14]: v∼2Dµ σ2.(4) In order to validate our algorithm (1) we conducted a series of discrete particle simulations of the Moran process with spatial diffusion, with diffusivity D. The continuum limit is the sFKPP equation for different values of µand σ2. In Figure 2 we plot the normalized speed, v/vF, of the Fisher wave versus the dimensionless noise strength σ/(Dµ)1/4as obtained from our discrete particle simulations. In the simulations the size of the domain was L= 100 with an interaction radius of δ= 1. We observe, in agreement with [11], a cross-over in the Fisher wave speed around the value unity for the normalized noise strength. The asymptotic predictions for the weak and strong noise limits, see Equations (3) and (4), are also captured with our discrete particle simulations. Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 4 0 0.1 0.2 0.3 0.4 0.5 0.6 -1 -0.5 0 0.5 1 1.5 2 FIXATION PROBABILTY µ N=10 N=10, analytical N=40 N=40, analytical Figure 1. The behaviour of the fixation probability as a function of the selective advantage µ, for population sizes N= 10 and N= 40. A comparison with the analytical prediction is shown (see text for more details). We have chosen the initial allele frequency as p=NA/N = 1/N in all simulations. 4. Results on the carrying capacity in a model toy flow A biologically relevant quantity is the carrying capacity (average population size in an ecosystem). In a recent paper [3] it was shown that a one-dimensional turbulent compressible velocity field leads to a reduction in the carrying capacity. This reduction of carrying capacity was later verified in a more realistic two-dimensional surface flow model of turbulence [4]. To investigate the behavior of our algorithms in presence of flow field we consider a very simple sinusoidal flow u=Usin(x). This flow is compressible and leads to a reduction in the carrying capacity on increasing the flow strength [15]. Below we consider a very simple model of population dynamics (the birth-coagulation process) whose mean-field limit also gives the Fisher equation. (A+Aχ −→ A Aγ −→ A+A where χdenotes the death (coagulation) rate and γis the birth rate. The continuum limit of the birth-coagulation process gives: ∂tc+∂x(uc) = D∇2c+µc(1 −c) + pσ2c(1 + c)Γ(t),(5) where cis the local population concentration and Γ(t) is a Gaussian white noise with zero mean and unit variance. We use the algorithm described in Section 2 to study the population dynamics Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 5 0.1 1 0.1 1 10 Normalized front speed Noise strength Weak noise Strong noise: Discrete particle simulation Figure 2. The normalized front speed of the Fisher wave (made dimensionless by vF) versus the dimensionless noise strength. Data from our particles ssimulations are the red filled circles, while the asymptotic conjectures for the weak noise (blue line, v∼√Dµ h2−π2 (log N)2i) and the strong noise (black line, v∼2Dµ σ2) limits are also shown. Each data point is obtained by doing an average over 50 independent noise realizations. in this simple flow with the above population rules. An important point to notice is that this model, at variance with the Moran model, allows for fluctuations in the total population size. In absence of any flow the population reaches a steady state with a population size N0, while the population size can changes drastically due to the presence of the sinusoidal flow field. We define the carrying capacity as: Zp(t)≡Nav N0 (6) where Nav is the average population size in the presence of the flow. Figure 3 shows that the carrying capacity Zpdecreases strongly at increasing the forcing strength. This is consistent with earlier 1dand 2dsimulations of population dynamics in turbulence flow field [3, 4] as well as with the stochastic simulations from [15]. 5. Conclusion We have introduced and discussed simple algorithms to study the dynamics of populations subject to an advecting flow. We benchmarked the numerical implementation amongst variants of the algorithm itself as well as against theoretical results. We show that the compressible flow reduces the total population size and thereby the carrying capacity. The proposed algorithms are a natural choice to investigate population genetics in flows. Future studies would look at the effect of population genetics in higher dimensional flows which are close to realistic settings Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 <Zp> U Algorithm1 Algorithm2 Stochastic Figure 3. Carrying capacity versus the forcing strength of the sinusoidal flow field, U, for our one-dimensional compressible toy flow. We observe that the carrying capacity dramatically decrease with increasing forcing strength and saturates for U≥0.3. Algorithm 1 and algorithm 2 are the two algorithms described above, while stochastic refers to the results obtained from the numerical integration of Equation 5. for e.g., surface flows. This would aid in our understanding of fixation of species in natural environments where fluid flow plays an important role. 6. Acknowledgement We thank M. H. Jensen and D.R. Nelson for discussions. We acknowledge the COST Action MP0806 for support. FT and PP acknowledge the Kavli Institute of Theoretical Physics for hospitality. This research was supported in part by the National Science Foundation under Grant No. NSF PHY05-51164. 7. Reference [1] K.S. Korolev, M. Avlund, O. Hallatschek, and D.R. Nelson, Genetic demixing and evolution in linear stepping stone models, Rev. Mod. Phys. 82,1691 (2010). [2] Mathias Sandulescu, Crist´obal L´opez, Emilio Hern´andez-Garc´ ia, and Ulrike Feudel, Biological activity in the wake of an island close to a coastal upwelling, Ecological Complexity 5, 228 (2008). [3] R. Benzi and D.R. Nelson, Fisher equation with turbulence in one dimension, Physica D (Amsterdam) 238, 2003 (2009). [4] P. Perlekar, R. Benzi, D.R. Nelson, and F. Toschi, Population dynamics at high Reynolds number, Phys. Rev. Lett. 101, 144501 (2010). [5] S. Berti, C. Lopez, D. Vergni, and A. Vulpiani, Discreteness effects in a reacting system of particles with finite interaction radius, Phys. Rev. E 76, 031139 (2007). [6] P. Moran, The Statistical Processes of Evolutionary Theory. Clarendon Press (1962). Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 7 [7] R.A. Fisher, The wave of advance of advantageous genes, Ann. Eugenics 7, 353 (1967). [8] A. Kolmogorov, N. Petrovsky, N. Piscounov, Etude de l’equation de la diffusion avec croissance de la quantite de la matiere et son application a un plobleme biologique Moscow Univ. Math Bull. 1, 1 (1937). [9] M. Kimura, On the probability of fixation of mutant genes in a population, Genetics 47, 713 (1962). [10] O. Hallatschek and D.R. Nelson, Population genetics and range expansion 62,42 (2009). [11] C. R. Doering, C. Mueller, and P. Smereka, Interacting particles, the stochastic Fisher-Kolmogorov-PetrovskyPiscounov equation and duality Physica A (Amsterdam) 325, 243 (2003). [12] E. Brunet, B. Derrida, Shift in the velocity of a front due to a cutoff, Phys. Rev. E 56 2597 (1997). [13] D.A. Kessler, Z. Ner, L.M. Sanderm Front propagation: precursors, cutoffs, and structural stability, Phys. Rev. E 58 107 (1998). [14] O. Hallatschek and K.S. Korolev, Fisher Waves in the Strong Noise Limit, Phys. Rev. Lett. 103, 108103 (2009). [15] S. Pigolotti, R. Benzi, M.H. Jensen, and D.R. Nelson, Population genetics in compressible flows, arXiv:1106.3506v1. Particles in Turbulence 2011 IOP Publishing Journal of Physics: Conference Series 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 8