Particle algorithms for population dynamics in flows

Pa icle algo i hms o popula ion dynamics in lows This a icle has been downloaded om IOPscience. Please sc oll down o see he ull ex a icle. 2011 J. Phys.: Con . Se . 333 012013 (h p://iopscience.iop.o g/1742-6596/333/1/012013) Download de ails: IP Add ess: 147.83.119.146 The a icle was downloaded on 30/08/2012 a 11:15 Please no e ha e ms and condi ions apply. View he able o con en s o his issue, o go o he jou nal homepage o mo e Home Sea ch Collec ions Jou nals Abou Con ac us My IOPscience Pa icle algo i hms o popula ion dynamics in lows P asad Pe leka 1, Robe o Benzi2, Simone Pigolo i3, and Fede ico Toschi1 1Depa men o Physics and Depa men o Ma hema ics and Compu e Science, Eindho en Uni e si y o Technology, Eindho en 5600MB, The Ne he lands. CNR-IAC, Via dei Tau ini 19, 00185 Rome, I aly. 2Dipa imen o di Fisica and INFN, Uni e si `a “To Ve ga a”, Via della Rice ca Scien i ica 1, I-00133 Roma, I aly. 3Depa amen de Fisica i Enginye ia Nuclea , Uni e si a Poli `ecnica de Ca alunya Edi . GAIA, Rambla San Neb idi s/n, 08222 Te assa, Ba celona, Spain. E-mail: . oschi@ ue.nl Abs ac . We p esen and discuss pa icle based algo i hms o nume ically s udy he dynamics o popula ion subjec ed o an ad ec ing low condi ion. We discuss ew possible a ian s o he algo i hms and compa e hem in a model comp essible low. A compa ison agains app op ia e e sions o he con inuum s ochas ic Fishe equa ion (sFKPP) is also p esen ed and discussed. The algo i hms can be used o s udy popula ions gene ics in luid en i onmen s. 1. In oduc ion Bac e ial colonies g owing on he su ace o a Pe i dish o m in iguing pa e ns ha can be ma hema ically modeled and s udied in e ms o a s epping s one model o popula ion gene ics [1]. Oceans co e a as amoun o he Ea h’s su ace suppo ing an inc edibly di e si y o species, playing a key ole since ea ly s age o li e on ou plane . Many mic oo ganisms mus ha e ound ways o su i e in he lows gene a ed in ocean. The e o e, hough no much explo ed, i is impo an o unde s and he ole o lows on models o he dynamic o popula ions. To unde s and his ques ion wo app oaches can be aken: (a) Use con inuum models and couple hem wi h luid lows [3, 4]. This is an app oach used commonly in oceanog aphy models o s udy plank on blooms [2]; (b) Use pa icle based models coupled o lows. The ad an age o his app oach is ha i allows us o pose in a simple way ques ions ela ed o he popula ions o numbe luc ua ions and he ex inc ion o a mu an species, which in u n is cen al o he unde s anding o popula ion gene ics [1, 10]. A simila app oach has been employed ea lie o s udy he ole o disc e e e ec s in he p opaga ion o a Fishe wa e [5]. In his pape we discuss algo i hms o s udy pa icle-based popula ion dynamics in lows, ollowing he abo e app oach (b). An en i y based algo i hm wi h di e en a ian s is discussed in sec ion 2. We hen alida e ou algo i hm o he Mo an p ocess in sec ion 3, whe e we i s p esen he esul s o he case o uni o mly mixed popula ions wi h ze o di usi i y and compa e ou esul s agains heo e ical p edic ions. Nex we ocus on he spa ial di usi i y. The con inuum limi o Mo an model wi h spa ial di usi i y is he s ochas ic Fishe -Kolmogo o - Pe o sky-Piscouno equa ion (sFKPP) [7, 8, 11]. We s udy he ole o popula ion luc ua ions on he speed o a Fishe wa e and compa e ou esul s agains he asymp o ic p edic ions o Re . [11]. Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 Published unde licence by IOP Publishing L d 1 Finally, in sec ion 4 we s udy he eac ion A→A+A,A+A→Ain a low, which also has he Fishe equa ion as a con inuum limi , bu wi h a di e en noise e m [11, 15]. We show ha in p esence o a simple, de e minis ic, comp essible low he a e age popula ion size (ca ying capaci y) dec eases, in ag eemen wi h ea lie in es iga ions [4, 15]. 2. Nume ical algo i hm One o he example o popula ion dynamics is he g ow h o bac e ia on a Pe i dish. The mo ion o he on ie o a bac e ial colony can be unde s ood as a solu ion o he s ochas ic Fishe equa ion. Ano he way o s udy on ie g ow h is by encoding he bi h and dea h ules a an indi idual en i y le el. The ad an age o his app oach is ha biologically ele an ques ions abou he bi h and ex inc ion o species can be easily add essed. We now desc ibe he wo algo i hms used by us o s udy popula ion gene ics wi h lows. 2.1. Algo i hm 1 The nume ical algo i hm ha is used o sol e he popula ion dynamics is simila o he algo i hm used o s udy he on p opaga ion o Fishe equa ion in Re . [5]. We conside each indi idual in he popula ion as a pa icle ha is ad ec ed by he low, di used by he B ownian mo ion, and ha can ep oduce i sel , die, o compe e wi h o he pa icles. Bi h, dea h and compe i ion can be modeled in e ms o bina y eac ions. In o de o e icien ly implemen hese bina y eac ions (which necessa ily happen when he wo pa icles su icien ly close o each o he ) we in oduce a spa ial mesh. Mo e p ecisely, we di ide ou one dimensional domain o size Lin o Msubin e als o wid h δ=L/M and whe e δis he in e ac ion dis ance. As a ypical s a ing condi ion one may conside pa icles uni o mly dis ibu ed a posi ions x∈L. The popula ion is e ol ed acco ding o he ollowing wo main s eps: •Ad ec ion and spa ial di usion: he pa icles a e ad ec ed and di used acco ding o xi( + ∆ ) = xi( ) + u∆ +√2D∆ Γi( ), i= 1,2,...n. He e Dis he spa ial di usion coe icien , uis he ad ec ing eloci y, ∆ is he ime s ep, nis he o al numbe o pa icles, and Γ( ) is a Gaussian whi e noise wi h ze o mean and uni a iance. •Popula ion dynamics: he pa icle popula ion is assumed homogeneous inside he subdomain δand in e ac wi h each o he using a p esc ibed se o ules de e mining he popula ion dynamics. Gi en he amewo k ha models ad ec ion, di usion and locali y o in e ac ion, one mus speci y he p ope eac ions which encode o he biological p ocesses in he popula ions. He e below we desc ibe how o implemen h ee o he building blocks a he base o he di e en eac ion models used in his manusc ip . Supposing ha ou ecosys em is composed by wo dis inc popula ions, Aand B, we deno e he numbe o pa icles wi hin he in e ac ion adius, δ, o ype Aas NAand hose o ype Bas NB. The o al numbe o pa icles in he in e al δ is indica ed wi h N=NA+NBwhe eas he o al numbe o pa icles in Lis indica ed wi h n. (i) A+Bk −→ B+B: In his ype o eac ion, a pa icle o one ype, A, in e ac s wi h a pa icle o he o he ype, B, and a a a e kcon e s i sel in o B. To implemen his ule we i s coun he numbe o Bpa icles p esen in δ. I kNB∆ < hen he pa icle Acon e s i sel in o Bo he wise no hing happens ( deno es a andom numbe uni o mly dis ibu ed in [0,1]). (ii) Ak −→ A+A: In his eac ion a pa icle o ype Agi es bi h o ano he pa icle o he same ype a a a e k. I k∆ < hen he pa icle mul iplies in o wo; o he wise no hing happens. Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 2 (iii) A+Ak −→ A: In his eac ion a pa icle o ype Ain e ac s wi h ano he pa icle o i s own kind and one o he wo dies wi h a a e k. Again we i s coun he numbe o Apa icles p esen in δ, i k(NA−1)∆ < hen he pa icle annihila es; o he wise no hing happens. 2.2. Algo i hm 2 This ype o algo i hm cons i u es a a ian o e algo i hm 1. He e he ad ec ion and di usion pa s a e implemen ed in he same way. Two main di e ences a e implemen ed in he popula ion dynamics pa o he algo i hm: G id size A each imes ep, he pa icle dis ibu ion is binned on a ine esolu ion. The bin size is chosen as δ/m, whe e δis he in e ac ion dis ance as be o e and mis an odd na u al numbe (in he ollowing we will choose m= 11). Bina y eac ions occu a a a e which depends on he o al numbe o neighbo s in he ow o mcells cen e ed on he bin whe e each pa icle is. Ra es The imes ep ∆ is aken small enough ha no mo e han one eac ion o each ype can occu wi hin each imes ep. Then, a ep oduc ion A→2Aoccu s wi h a p obabili y kn( )∆ whe e n( ) is he o al numbe o pa icles. A dea h by compe i ion, A+A→A occu s wi h a a e k∆ PjNj( ), whe e Nj( ) is he o al numbe o pa icles in he ow cen e ed in he bin whe e pa icle jis and he sum uns on all he pa icles p esen in he sys em. In he case o ep oduc ion, he ep oducing pa icle is chosen a andom wi h uni o m p obabili y, while in he case o dea h by compe i ion he choice is weigh ed wi h he numbe o neighbo s, i.e. pa icle ihas a p obabili y o be chosen equal o Ni/PjNj. The i s modi ica ion allow o ha e compe i ion in an a ea δmo e p ecisely cen e ed in he pa icle posi ion. We ound ha he second modi ica ion slows down simula ion conside ably o ou pa ame e alues. The e o e, i should be implemen ed only in cases whe e he imescales o popula ion dynamics a e much longe han hose o ad ec ion di usion, so ha one is no o ced o choose a e y small ime s ep in o de o a oid ha ing mo e han one eac ion pe ime s ep. 3. Valida ion The alida ion is di ided in wo pa s aimed a es ing, espec i ely, he eac ion algo i hm alone and he di usion- eac ion implemen a ion. We i s ocus on he eac ion pa and alida e ou popula ion dynamics algo i hm o he case o Mo an p ocess [6, 9] o popula ion gene ics. He e we compa e ou esul s wi h he heo e ical p edic ion on he ixa ion p obabili y om Kimu a [9] (see Subsec ion 3.1). In subsec ion 3.2 we add he e ec o he spa ial di usion: wi h his addi ion he Mo an p ocess becomes equi alen o a s ochas ic Fishe equa ion [7, 8, 11]. We compa e he esul s o ou disc e e simula ions agains he analy ical p edic ions based on he s ochas ic Fishe equa ion, bo h in s ong and weak noise limi [11, 14]. 3.1. Mo an model We conside a homogeneous popula ion o wo species Aand Bwi h NAindi iduals o ype-A and NBindi iduals o ype Bin a domain o size δ. When an indi idual o ype A(B) in e ac s wi h an en i y o he o he ype B(A) hen B(A) con e s in o A(B) wi h a a e kA(kB). These eac ion p ocesses can be w i en in o m he ollowing o m: (A+BkA −−→ A+A A+BkB −−→ B+B Du ing he abo e p ocesses he o al numbe o en i ies N≡NA+NB emain conse ed o all imes. The con inuum limi o he abo e p ocess is he s ochas ic equa ion dc = Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 3 µc(1 −c) + pσ2c(1 −c)dW (FKPP equa ion wi h noise) whe e c=NA/N is he concen a ion o he ype-Aspecies, µ≡(kA−kB)N/δ and σ2≡(kA+kB) δ,dW is a Wiene p ocess. One o he cen al ques ion in popula ion gene ics is he ixa ion p obabili y o a mu a ion. Fo wha conce ns he Mo an model his ques ion can be o mula ed as: i Bis he pa en species and Ais he mu an , wha is he p obabili y ha he mu an con e s all he pa en species B in o i sel A(i.e. Age s ixa ed)? The wo con ol pa ame e s a e he ini ial mu an popula ion concen a ion p=cA( = 0) = NA/N and he selec i e ad an age µ. Fo con enience, in his sec ion we se δ= 1. Fo he case o he Mo an p ocess Kimu a [9, 10] p edic ed ha he ixa ion p obabili y is: P ix =1−e−αNp 1−e−αN (1) whe e α≡kA−kB kA+kBand pis he ini ial concen a ion o he mu an popula ion. No e ha o kA=kBo µ= 0 (case co esponding o no selec i e ad an age), he ixa ion p obabili y P ix =p. The plo in Figu e 1 shows he ixa ion p obabili y o wo di e en popula ions o size Nwi h p= 1/N and a ying µ. The nega i e alues o he selec i e ad an age implies ha he pa en species has a selec i e ad an age o e he mu an . We ind ha , in ag eemen wi h he heo e ical p edic ion, he ixa ion p obabili y inc eases wi h he selec i e ad an age o he mu an . Fo he special case o ze o selec i e ad an age (kA=kB) we ind, in ag eemen wi h heo y, P ix =p. In ou simula ions he ixa ion p obabili y, P ix , is de ined as one minus he ac ion o ealiza ions o which he mu an dies. Ou measu emen s we e ob ained by a e aging o e an ensemble o 2000 ealiza ions o each da a poin . 3.2. S ochas ic Fishe equa ion To s udy he e ec o spa ial di usi i y, Fishe , Kolmogo o e al. [7, 8] added a di usion e m o he Fishe equa ion: dc = [µc(1 −c) + D∇2c]d +pσ2c(1 −c)dW. (2) In absence o noise, o he case o a localized an ini ial condi ion, he Fishe equa ion exhibi s a a eling wa e solu ion wi h a on speed gi en by F= 2√Dµ [7, 8]. Mo e ecen ly, i was shown ha he speed o he Fishe wa e is educed in p esence o noise. In he egime o weak noise [11, 12, 13, 14]: ∼pDµ 2−π2 (log N)2(3) whe e Nis he o al popula ion wi hin he in e ac ion adius δand σ2∝1/N. In he s ong noise egime he Fishe speed is d as ically educed [11, 14]: ∼2Dµ σ2.(4) In o de o alida e ou algo i hm (1) we conduc ed a se ies o disc e e pa icle simula ions o he Mo an p ocess wi h spa ial di usion, wi h di usi i y D. The con inuum limi is he sFKPP equa ion o di e en alues o µand σ2. In Figu e 2 we plo he no malized speed, / F, o he Fishe wa e e sus he dimensionless noise s eng h σ/(Dµ)1/4as ob ained om ou disc e e pa icle simula ions. In he simula ions he size o he domain was L= 100 wi h an in e ac ion adius o δ= 1. We obse e, in ag eemen wi h [11], a c oss-o e in he Fishe wa e speed a ound he alue uni y o he no malized noise s eng h. The asymp o ic p edic ions o he weak and s ong noise limi s, see Equa ions (3) and (4), a e also cap u ed wi h ou disc e e pa icle simula ions. Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 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, analy ical N=40 N=40, analy ical Figu e 1. The beha iou o he ixa ion p obabili y as a unc ion o he selec i e ad an age µ, o popula ion sizes N= 10 and N= 40. A compa ison wi h he analy ical p edic ion is shown (see ex o mo e de ails). We ha e chosen he ini ial allele equency as p=NA/N = 1/N in all simula ions. 4. Resul s on he ca ying capaci y in a model oy low A biologically ele an quan i y is he ca ying capaci y (a e age popula ion size in an ecosys em). In a ecen pape [3] i was shown ha a one-dimensional u bulen comp essible eloci y ield leads o a educ ion in he ca ying capaci y. This educ ion o ca ying capaci y was la e e i ied in a mo e ealis ic wo-dimensional su ace low model o u bulence [4]. To in es iga e he beha io o ou algo i hms in p esence o low ield we conside a e y simple sinusoidal low u=Usin(x). This low is comp essible and leads o a educ ion in he ca ying capaci y on inc easing he low s eng h [15]. Below we conside a e y simple model o popula ion dynamics ( he bi h-coagula ion p ocess) whose mean- ield limi also gi es he Fishe equa ion. (A+Aχ −→ A Aγ −→ A+A whe e χdeno es he dea h (coagula ion) a e and γis he bi h a e. The con inuum limi o he bi h-coagula ion p ocess gi es: ∂ c+∂x(uc) = D∇2c+µc(1 −c) + pσ2c(1 + c)Γ( ),(5) whe e cis he local popula ion concen a ion and Γ( ) is a Gaussian whi e noise wi h ze o mean and uni a iance. We use he algo i hm desc ibed in Sec ion 2 o s udy he popula ion dynamics Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 5 0.1 1 0.1 1 10 No malized on speed Noise s eng h Weak noise S ong noise: Disc e e pa icle simula ion Figu e 2. The no malized on speed o he Fishe wa e (made dimensionless by F) e sus he dimensionless noise s eng h. Da a om ou pa icles ssimula ions a e he ed illed ci cles, while he asymp o ic conjec u es o he weak noise (blue line, ∼√Dµ h2−π2 (log N)2i) and he s ong noise (black line, ∼2Dµ σ2) limi s a e also shown. Each da a poin is ob ained by doing an a e age o e 50 independen noise ealiza ions. in his simple low wi h he abo e popula ion ules. An impo an poin o no ice is ha his model, a a iance wi h he Mo an model, allows o luc ua ions in he o al popula ion size. In absence o any low he popula ion eaches a s eady s a e wi h a popula ion size N0, while he popula ion size can changes d as ically due o he p esence o he sinusoidal low ield. We de ine he ca ying capaci y as: Zp( )≡Na N0 (6) whe e Na is he a e age popula ion size in he p esence o he low. Figu e 3 shows ha he ca ying capaci y Zpdec eases s ongly a inc easing he o cing s eng h. This is consis en wi h ea lie 1dand 2dsimula ions o popula ion dynamics in u bulence low ield [3, 4] as well as wi h he s ochas ic simula ions om [15]. 5. Conclusion We ha e in oduced and discussed simple algo i hms o s udy he dynamics o popula ions subjec o an ad ec ing low. We benchma ked he nume ical implemen a ion amongs a ian s o he algo i hm i sel as well as agains heo e ical esul s. We show ha he comp essible low educes he o al popula ion size and he eby he ca ying capaci y. The p oposed algo i hms a e a na u al choice o in es iga e popula ion gene ics in lows. Fu u e s udies would look a he e ec o popula ion gene ics in highe dimensional lows which a e close o ealis ic se ings Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 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 Algo i hm1 Algo i hm2 S ochas ic Figu e 3. Ca ying capaci y e sus he o cing s eng h o he sinusoidal low ield, U, o ou one-dimensional comp essible oy low. We obse e ha he ca ying capaci y d ama ically dec ease wi h inc easing o cing s eng h and sa u a es o U≥0.3. Algo i hm 1 and algo i hm 2 a e he wo algo i hms desc ibed abo e, while s ochas ic e e s o he esul s ob ained om he nume ical in eg a ion o Equa ion 5. o e.g., su ace lows. This would aid in ou unde s anding o ixa ion o species in na u al en i onmen s whe e luid low plays an impo an ole. 6. Acknowledgemen We hank M. H. Jensen and D.R. Nelson o discussions. We acknowledge he COST Ac ion MP0806 o suppo . FT and PP acknowledge he Ka li Ins i u e o Theo e ical Physics o hospi ali y. This esea ch was suppo ed in pa by he Na ional Science Founda ion unde G an No. NSF PHY05-51164. 7. Re e ence [1] K.S. Ko ole , M. A lund, O. Halla schek, and D.R. Nelson, Gene ic demixing and e olu ion in linea s epping s one models, Re . Mod. Phys. 82,1691 (2010). [2] Ma hias Sandulescu, C is ´obal L´opez, Emilio He n´andez-Ga c´ ia, and Ul ike Feudel, Biological ac i i y in he wake o an island close o a coas al upwelling, Ecological Complexi y 5, 228 (2008). [3] R. Benzi and D.R. Nelson, Fishe equa ion wi h u bulence in one dimension, Physica D (Ams e dam) 238, 2003 (2009). [4] P. Pe leka , R. Benzi, D.R. Nelson, and F. Toschi, Popula ion dynamics a high Reynolds numbe , Phys. Re . Le . 101, 144501 (2010). [5] S. Be i, C. Lopez, D. Ve gni, and A. Vulpiani, Disc e eness e ec s in a eac ing sys em o pa icles wi h ini e in e ac ion adius, Phys. Re . E 76, 031139 (2007). [6] P. Mo an, The S a is ical P ocesses o E olu iona y Theo y. Cla endon P ess (1962). Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 7 [7] R.A. Fishe , The wa e o ad ance o ad an ageous genes, Ann. Eugenics 7, 353 (1967). [8] A. Kolmogo o , N. Pe o sky, N. Piscouno , E ude de l’equa ion de la di usion a ec c oissance de la quan i e de la ma ie e e son applica ion a un plobleme biologique Moscow Uni . Ma h Bull. 1, 1 (1937). [9] M. Kimu a, On he p obabili y o ixa ion o mu an genes in a popula ion, Gene ics 47, 713 (1962). [10] O. Halla schek and D.R. Nelson, Popula ion gene ics and ange expansion 62,42 (2009). [11] C. R. Doe ing, C. Muelle , and P. Sme eka, In e ac ing pa icles, he s ochas ic Fishe -Kolmogo o -Pe o sky- Piscouno equa ion and duali y Physica A (Ams e dam) 325, 243 (2003). [12] E. B une , B. De ida, Shi in he eloci y o a on due o a cu o , Phys. Re . E 56 2597 (1997). [13] D.A. Kessle , Z. Ne , L.M. Sande m F on p opaga ion: p ecu so s, cu o s, and s uc u al s abili y, Phys. Re . E 58 107 (1998). [14] O. Halla schek and K.S. Ko ole , Fishe Wa es in he S ong Noise Limi , Phys. Re . Le . 103, 108103 (2009). [15] S. Pigolo i, R. Benzi, M.H. Jensen, and D.R. Nelson, Popula ion gene ics in comp essible lows, a Xi :1106.3506 1. Pa icles in Tu bulence 2011 IOP Publishing Jou nal o Physics: Con e ence Se ies 333 (2011) 012013 doi:10.1088/1742-6596/333/1/012013 8