Discrete Eulerian model for population genetics and dynamics under flow

TL;DR

This paper introduces a novel algorithm for simulating population genetics in fluid flows, validated in 1D and 2D, revealing increased fixation probabilities in certain flow conditions, with implications for marine species dynamics.

Contribution

The paper presents a new algorithm that efficiently models population genetics under fluid flow, validated in multiple dimensions, and applied to marine species in weakly compressible flows.

Findings

Organisms at sources in flow have higher fixation probabilities.

Algorithm accurately models advection and fluctuations in population dynamics.

Flow conditions significantly influence genetic fixation outcomes.

Abstract

Marine species reproduce and compete while being advected by turbulent flows. It is largely unknown, both theoretically and experimentally, how population dynamics and genetics are changed by the presence of fluid flows. Discrete agent-based simulations in continuous space allow for accurate treatment of advection and number fluctuations, but can be computationally expensive for even modest organism densities. In this report, we propose an algorithm to overcome some of these challenges. We first provide a thorough validation of the algorithm in one and two dimensions without flow. Next, we focus on the case of weakly compressible flows in two dimensions. This models organisms such as phytoplankton living at a specific depth in the three-dimensional, incompressible ocean experiencing upwelling and/or downwelling events. We show that organisms born at sources in a two-dimensional time-independent flow experience an increase in fixation probability.

Equations

$$P_{\text{fix}}=\frac{1-e^{-sNf}}{1-e^{-sN}}.$$

$$p\equiv \frac{D\Delta t}{\Delta x^2};$$

$$x_{a}(i)=\left(i-\frac{1}{2}\right)\Delta x+\Delta x\bigg{(}\eta-\frac{1}{2}\bigg{)}(1+2p),$$

$$x_a(i,t+\Delta t)=x_a(i)+u(x_a(i))\Delta t$$

$$j\equiv \lfloor \frac{x_a(i,t+\Delta t)}{\Delta x}\rfloor +1$$

$${\tilde{N}}_j={\tilde{N}}_j+1$$

$${\tilde{N}}_j^{(A)}$$

$$r_b(A)$$

$$r_d(A)$$

$${\tilde{N}}_j^{(B)}$$

$$r_b(B)$$

$$r_d(B)$$

$${\partial }_{t}N(x,t)+{\partial }_x(u(x,t)N(x,t))=D\Delta N.$$

$$\frac{d{N}_A(t)}{dt}$$

$$\frac{d{N}_B(t)}{dt}$$

$${\partial }_t{c}_A+{\partial }_x(u{c}_A)=D\Delta c_A+\mu c_A(1-c_A-c_B)$$

$$+s\mu c_A{c}_B+\sqrt{\mu \frac{c_A}{N_0}(1+c_A+c_B)}{\eta }_A(x,t)$$

$${\partial }_t{c}_B+{\partial }_x(u{c}_B)=D\Delta c_B+\mu c_B(1-c_A-c_B)$$

$$-s\mu c_A{c}_B+\sqrt{\mu \frac{c_B}{N_0}(1+c_A+c_B)}{\eta }_B(x,t)$$

$${\partial }_{t}f+u(x,t){\partial }_{x}f=D\Delta f+s\mu f(1-f)+\sqrt{\frac{2\mu f(1-f)}{N_0}}\eta (x,t).$$

$${\partial }_{t}c=D{\partial }_{xx}c+\mu c(1-c),$$

$$c(x,t)=\omega (x-vt),$$

$$D{\omega }^{\prime \prime }+v{\omega }^{\prime }+\mu \omega (1-\omega )=0,$$

$$\text{with}\text{conditions}\omega (-\infty )=1,\omega (\infty )=0.$$

$$v\sim\sqrt{D\mu}\bigg{[}2-\frac{\pi^{2}}{(\ln N)^{2}}\bigg{]}.$$

$$v\sim D\mu N.$$

$$H(x,x^{\prime };t)\equiv ⟨f(x,t)[1-f(x^{\prime },t)]⟩.$$

$$u(x)=u_{0}\sin\bigg{(}x-\frac{\pi}{2}\bigg{)}.$$

$$N_{\text{eff}}=B_1{\rho }_0\sqrt{\frac{D}{u_0}},$$

$$u(x)$$

$$u(y)$$

$$N_{\text{eff}}=2B_1^2{\rho }_0\frac{D}{u_0}$$

$$N_g\sim 4{\delta }^2{\rho }_0=4{\left[\frac{2\sqrt{D\mu s}}{u_0}\right]}^2{\rho }_0.$$

Full text

I I. Introduction

Marine plankton account for roughly half of the total biological production on Earth; they are responsible for most of the transfer of carbon dioxide from the atmosphere to the ocean  guasto2012fluid ; field1998primary . Planktonic organisms are an essential part of the global carbon cycle, and even small changes in their productivity or in the relative abundance of the thousands of species could have a substantial influence on climate change  field1998primary . It is important to understand the variation in physical factors that a population can withstand and how it can continue to thrive in high Reynolds number fluid environments in order to support our oceanic ecosystem  perlekar2010population .

Microorganism populations are carried along the uppermost layer (euphotic zone $\sim$100 m) of the ocean  d2010fluid .

The euphotic zone is characterized by a low quantity of nutrients due to consumption by phytoplankton. Periodic events, such as upwelling and downwelling currents, supply nutrients to the upper water column. The aforementioned mechanisms can trigger the processes of water exchange in the mixed layer of the ocean. The upwelling current leads to a rising up of deep water, where a rich concentration of nutrients resides. Passively traveling organisms, transported by the ocean circulation, experience compressible turbulence  benzi2012population from the convergence or the divergence of water masses.

The study of genetic variation within a population deals with the biological differences affecting reproduction, feeding strategies, disease resistance, and many other factors. Well-adapted individuals with inherited favorable characteristics may survive and grow faster than others, passing on the genes that make them successful; such organisms have a selective advantage.

If we consider two species, one with a selective advantage and one without, called $A$ and $B$, respectively, it is not possible to determinate a priori which one of the two will be the dominant one in the presence of turbulence. However it is possible to calculate a probability. In the absence of advection, Kimura  kimura1962probability derived a theoretical prediction for the fixation probability of one of the species for the well-mixed case,

[TABLE]

This formula describes the fixation probability for a species with selective advantage $s$ in a population of size $N$ that makes up an initial fraction $f$ of all organisms, neglecting any space dependency. This result can be applied to a spatially extended population with simple migration patterns, such as diffusion maruyama1974simple . Several stochastic models for genetic evolution have been developed. Among these, we mention the Moran model moran1958random , a simple approach that takes into consideration the selection of the organisms and the genetic drift, and the Stepping Stone model korolev2010genetic , an extension of the previous one including the migration and reproduction of the individuals. These models share many similarities with the ones used to investigate nonequilibrium phase transitions (see hinrichsen2000non ; dornic2005integration ). The aforementioned models are tailored for lattice rules and do not allow straightforward generalization to take into account an external velocity advection. In pigolotti2013growth , an alternative method has been introduced: each individual is advected by the external velocity and diffusion is implemented by a stochastic noise, while death and reproduction processes are implemented using an interaction distance $\delta$. This requires an extra computational cost to evaluate the individual numbers in each virtual deme of size $\delta$. This method has been recently used in  plummer2018fixation , where the competition between two different species, distributed in continuous space, and under the effect of a compressible flow is examined through an agent-based model. It has been shown that a turbulent flow can dramatically change outcomes and, in particular, it can reduce the effect of selective advantage on fixation probabilities  plummer2018fixation ; pigolotti2012population ; benzi2009fisher .

In this paper, we propose a computational approach which merges the accuracy of working in continuous space with the efficiency of working on a lattice. We assume a uniform lattice of spacing $\Delta x$ with each site occupied by $N_{j}$ individuals. At each time step, we redistribute the $N_{j}$ individuals on a domain $(1+a)\Delta x^{d}$, where $d$ is the dimension of our system and $a$ is suitably chosen to introduce a diffusion process (see next section). Next, we advect the $N_{j}$ individuals in continuous space using the external velocity (if present). After this step, some of the original $N_{j}$ individuals have been moved to different regions of space, i.e., to a different box of size $\Delta x$, changing the number of individuals of the new box. Once we complete the diffusion and advection for all sites, independently one from another, we apply the birth-death processes stochastically according to the prescribed rates. Note that we do not need to remember the exact position in each site from one step to another: it is enough to know how many individuals of each species are present at the prescribed site. In this way, we can efficiently work with an extremely large number of individuals per site without managing the position of each individual. This is actually the reason why we can achieve a significant increase in the computational performance.

The paper is organized as follows. In Sec. II, we discuss the details of our method. We present a systematic comparison of our approach against known analytical and numerical results in one dimension (Sec. III) and two dimensions (Sec. IV). In Sec. V, this approach is used to extend the previous findings of  plummer2018fixation . In particular, we investigate the fixation probability of an advantageous species in a two-dimensional weak compressible flow.

II II. Method

The computational approach is described in this section, where, for simplicity, we start by considering a one-dimensional (1D) system with periodic boundary conditions. Let $L$ be the size of the 1D lattice, which we discretize with $n$ intervals of size $\Delta x=L/n$. Each interval $i=1,....n$ spans the region $x\in[(i-1)\Delta x,i\Delta x]$. We denote by $N_{i}^{(\beta)}$ the number of individuals in the interval $i$, where $\beta=A,B$ refers, in this case, to the two possible species (for different realizations, the number of the species may also be greater). At equilibrium and with no external flow, we can define the density $N_{0}$ of individuals per mesh point corresponding to the overall carrying capacity $N_{0}L/\Delta x$. With this definition, we can also think of $N_{0}$ as the average carrying capacity for a mesh site.

At time $t$, our knowledge is given by the set of numbers $N_{i}^{(\beta)}$ for $i=1,...,n.$ Our task is to compute the evolution of the system at time $t+\Delta t$, where $\Delta t$ represents our time step.

We implement the evolution using four different steps. In step 1, we implement a Markov chain with next-neighbor hopping and periodic boundary conditions, which is known to be consistent with the diffusion equation with diffusivity $D$ once the hopping probability is given by the relation

[TABLE]

with $p\ll 1$.

Step 1. $\it{Diffusion}$. For each interval $i$, we compute the particle positions $x_{a}(i)\,(a=1,...,N_{i})$ according to the rule:

[TABLE]

where $\eta$ is a random number that is uniformly distributed $[0,1]$. In this step, only a small fraction of the $N$ individuals is spread outside of the initial site $i$. Note that we do not assume any knowledge of the previous position of the individuals.

Step 2. $\it{Advection}$. Once step 1 is performed, we can compute the advection and obtain

[TABLE]

where $u(x,t)$ is a prescribed advecting field.

Step 3. $\it{Relabeling}$. For each off-mesh particle $a$, we can now determine the deme index,

[TABLE]

and therefore apply the rule

[TABLE]

to increment the deme occupancy number. Note that before implementing Eq. (3) we put ${\tilde{N}}_{j}=0$ for all demes $j=1,....,n$. Since from Step 1 to Step 3 we repeat the same operation for both species in the sections, we ignore the label $\beta$ for the different species.

Step 4. $\it{Birth\,and\,death\,processes}$. After running step 1 to step 3 for all the intervals, we apply the last step where we execute the rules for stochastic population dynamics for each segment $j$. At this stage, for every $j$ interval, we compute the birth-death process $\tilde{N_{j}}$ times according to the following rules:

[TABLE]

where $s$ is the selective advantage, $s>0$, or disadvantage, $s<0$, of individuals $A$ with respect to $B$. Here, $r_{b}$ and $r_{d}$ denote the $birth$ and the $death$ probability, respectively. Note that for each mesh site, the probability to obtain $k$ new offspring or deaths is binomial and it approximates a Poisson distribution only when the number of individuals considered in the specific process is large enough. This is never the case near the edge of a propagating front and/or near extinction even for large value of $N_{0}$.

At the end of step 4, we can put ${N}_{j}^{(\alpha)}={\tilde{N}}_{j}^{(\alpha)}$ and we can start with a new time step.

Let us now briefly comment about our method. The effect of advection does not change the number of particles, i.e., it is conservative. Thus neglecting, for the time being, the death-birth process, we obtain the equation for each species,

[TABLE]

The birth-death process is the same one implemented in Ref.  plummer2018fixation .

On the other side, ignoring diffusion and advection and neglecting terms of the order of $s/N_{0}$ inside the noise term  pigolotti2013growth , step 4 gives

[TABLE]

where $\eta_{A}$ and $\eta_{B}$ are independent $\delta$ correlated in time Wiener processes. Upon defining $c_{A}=N_{A}/N_{0}$ and $c_{B}=N_{B}/N_{0}$ and introducing the advection the final equations of motion read:

[TABLE]

Finally, assuming that $c_{T}\equiv c_{A}+c_{B}\sim 1$ everywhere and upon denoting $f=c_{A}/c_{T}$, we obtain

[TABLE]

We remark that the statistical properties of the system are invariant upon the scaling: $\mu\rightarrow 1$, $D\rightarrow\frac{D}{\mu}$, $t\rightarrow t\mu$, which is equivalent to working in units of generation time.

III III. Numerical test in one dimension

In this section, we introduce some numerical tests confined to one-dimensional systems. First, we need to solve the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation that describes the space-time evolution of a population in a reaction-diffusion system; in one space dimension, it reads

[TABLE]

where $c(x,t)$ is a continuous variable that identifies the concentration of individuals, $D$ is the diffusion coefficient, and $\mu$ is the growth rate. The uniform solutions of Eq. (17) are $c=1$ and $c=0$ for a stable and unstable state, respectively. In 1995, Mueller and Sowers  mueller1995random showed that for $\mu>0$, the traveling wave solutions to Eq. (17) are always characterized by a $\it{compact\,support\,property}$. We can set up initial conditions that depend on $c$ as follows: $c(x,0)\rightarrow 1$ as $x\rightarrow-\infty$ and $c(x,0)\rightarrow 0$ as $x\rightarrow+\infty$. For this kind of boundary conditions, we can find a continuous family $\omega$ of traveling wave solutions of the form

[TABLE]

where $v$ is the velocity of the traveling wave and $\omega$ is a function that must satisfy the following ordinary differential equation

[TABLE]

Around the unstable state, $c(x,0)$, the velocity of the front approaches the deterministic continuum minimum value $v_{\text{min}}=2\sqrt{D\mu}$. The fronts at this minimum speed are called “pulled fronts”, which are pulled along by the growth and spreading of small perturbations in the leading edge where $c\ll 1$. We therefore expect this velocity to change due to the discreteness of our model: we are in the presence of a discrete process in both time and space and the observed value for the Fisher wave velocity propagation is lower than the deterministic one. Brunet and Derrida  brunet gave an estimation of how far the Fisher wave value has to be from the continuum wave speed as

[TABLE]

From Eq. (20), one can clearly observe that the convergence to the continuum limit is extremely slow as $N\rightarrow\infty$. Fluctuations have been considered by the Doering, et al. conjecture  doering2003interacting adding a noise term to the FKPP equation; for the strong noise regime (or weak growth limit), they found that the speed value goes according to

[TABLE]

In Fig.  1 the normalized Fisher wave speed versus the number of individuals per site $N_{0}$ is shown. There are two theoretical estimates, corresponding to the weak and strong limits. The simulations, identified by dots, are consistent with the theoretical lines: with $10$ particles per site we are in the strong noise regime, where the Fisher velocity is equal to $\approx$ 0.3 times the theoretical expectation. In this work, simulations are performed in the weak noise regime, where the velocity of the genetic wave $v_{g}$ is $\simeq v_{\text{min}}$.

The diversity of a population composed of two genotypes in one dimension is measured by the heterozygosity  korolev2010genetic ,

[TABLE]

This quantity is given by the product of the two fractions $f(x,t)$ and $(1-f(x^{\prime},t))$ and it defines the probability that two selected individuals, chosen at random, are from different species (carry different alleles)  pigolotti2013growth . For homogeneous conditions, $H(x,x^{\prime};t)$ depends on the $r=|x-x^{\prime}|$. The heterozygosity becomes zero when there is fixation of one of the two genotypes. Moreover, it is known that in a one-dimensional system, $H(t)\equiv H(x,x;t)$ decays in time as $t^{-1/2}$. In Fig.  2, we have tested this theoretical prediction using our methods with $N_{0}=50$ on a domain with periodic boundary conditions discretized with 512 mesh points: the result very clearly confirms the theoretical behavior.

Next, to further validate the algorithm, we calculate, in the absence of advection, the fixation probability given by Eq. (1).

In Fig. 3, different panels corresponding to a different number of individuals per box are shown. In our simulations, we focus on small selective advantages, in order to study more realistic cases. Our results are in good agreement with the theoretical predictions (continuous black lines, in the figures).

IV IV. Numerical test in two dimensions

In this section, we implement the method previously introduced (Sec. II) and validated for a one-dimensional system on a two-dimensional configuration. Following the same schematic procedure of the 1D case we start by estimating the heterozygosity parameter. It is known that in two spatial dimensions, the local heterozygosity decay in time is slower compared to 1D: it goes to zero as $H(t)\sim 1/\ln(t)$  korolev2010genetic ; pigolotti2013growth . To check whether our method is able to exhibit such (slow) decay, we specifically perform a set of numerical simulations with $N_{0}=20$ on a domain with periodic boundary conditions and $256^{2}$ mesh point. In Fig.  4, such slow logarithmic decay is appreciable. In this figure, we plot $1/H(t)$ versus time. Note that starting with well mixed conditions, $H(0)=1/4$. Therefore, $1/H(t)$ is $4$ at $t=0$ and grows in time as $\ln(t)$, as shown in the figure. The loss of the genetic variability given by our simulations (purple triangles) is in agreement with the theory (black solid line).

The second step, as in the numerical validations in 1D, is to verify Kimura’s formula, given by Eq. (1), for the two-dimensional system in the absence of advection. The formula for the probability of fixation is still valid for higher dimensions and our results together with the theoretical prediction (solid black line) show an unequivocal agreement in Fig.  5.

V V. Weak compressible flow in D=2

Before adding an advecting velocity field to our two-dimensional system, we briefly discuss the main results achieved by Plummer et al.  plummer2018fixation , where a particular configuration of the velocity field was used, given by

[TABLE]

For small enough $u_{0}$, the flow field in (23) is weakly compressible, i.e., the condition $c_{A}+c_{B}=1$ is valid within a small percentage (up to $4$ percent for $u_{0}=0.05$ on a domain size $2\pi$). We will test whether, as in 1D, the Kimura formula is still valid provided we define $N$ as an effective population size, $N_{\text{eff}}$. For $s\rightarrow 0$, it has been shown in  plummer2018fixation that $N_{\text{eff}}$ depends only on the diffusion constant $D$, $u_{0}$, and on the maximum number $N_{0}$ of individuals per site. The crucial point is to recognize that near to the source, one can define a characteristic scale, $l_{s}=\sqrt{D/u_{0}}$. Any organism that moves significantly farther than $l_{s}$ from the source is unlikely to be able to return and has, therefore, a negligible chance of fixation as it is drawn into the sink. It follows (see  plummer2018fixation for details) that $N_{\text{eff}}$ can be estimated as

[TABLE]

where $B_{1}$ is a constant of the order of unity and $\rho_{0}$ is the density at each point, namely $N_{0}/\Delta x$.

Following plummer2018fixation , one simple way to understand the physical meaning of Eq. (24) is to consider the deterministic case, i.e., Eq. (16) in the limit $N_{0}\rightarrow\infty$, and assume an initial population $f=1$ in a small box $\Delta$ at the location $x_{0}$, and zero otherwise. Then, the population, whose spatially averaged initial ratio is $\bar{f}_{0}\equiv\langle f_{0}\rangle=\Delta/L$, where $L=2\pi$, evolves to an asymptotic value, $f_{\infty}(x_{0})=lim_{t\rightarrow\infty}f_{t}(x_{0})$, which depends on $x_{0}$. The ratio $f_{\infty}(x_{0})/\bar{f}_{0}$ is a function of $x_{0}$ and shows a Gaussian-like behavior in terms of $x_{0}-x_{s}$, where $x_{s}$ is the position of the source with a variance proportional to $l_{s}$ and $f_{\infty}(x_{s})/\bar{f}_{0}\gg 1$. This is equivalent to saying that for $s=0$, there is a significant advantage for the offspring occurring near the source and a strong disadvantage for those occurring downstream. This implies that the effective population size (for small $s$) is the one corresponding to the population size close to the source, i.e., at distance $l_{s}$ from the source. Following lieberman , one simple way to understand this result is to consider a simple toy model on a linear graph where the source is a relative “cold” site (node of the graph) with respect to the downstream “hot” sites, where “cold” and “hot” refer to the probability for an offspring to be advected by the external flow using the same language of lieberman .

The same considerations can be made for the two-dimensional version of the same problem. For this purpose, we consider the following flow

[TABLE]

with periodic boundary conditions and a domain of size $(2\pi)^{2}$. In Fig.  6, the final fraction of the initially localized species $f_{\infty}(x_{0},y_{0})/\bar{f}_{0}$ is shown where, now, $\bar{f}_{0}=\Delta^{2}/(4\pi^{2})$. Two peaks are clearly visible in correspondence with the sources, representing the upwelling regions. The asymptotic value of $f_{\infty}/\bar{f}_{0}$ is increasing in proximity of the sources while being it is reduced moving away from them. In Fig.  7, we show, with a black line, a one-dimensional section (along the $y$ axis) of the two-dimensional behavior of $f_{\infty}/\bar{f}_{0}$. Since for $s=0$ $P_{\text{fix}}=\bar{f}_{0}$, one can consider the black line as the increase in $P_{\text{fix}}$ due to the effect of the velocity field near to the source. To validate this interpretation, as well as the quality of our method, we performed a series of numerical simulations with $N_{0}=2$ at $s=0$ using the same initial conditions of the deterministic simulations. After estimating the fixation probabilities, we compute the increase of $P_{\text{fix}}$ as a function of the initial position, $x_{0}$. The results are shown as symbols in Fig.  7 where an excellent agreement is visible with the deterministic value of $f_{\infty}/\bar{f}_{0}$. This result demonstrates that the mechanism described in  plummer2018fixation , for small enough $s$ should be true for the two-dimensional flow considered here.

Based on the previous results, we can generalize Eq. (24) for the two dimensional case as follows:

[TABLE]

The factor $2$ in Eq. (26) comes from the fact that for our flow field, given by Eq. (25), we have two sources and two sinks. Using a grid resolution of $64^{2}$, with $N_{0}=2$, we have computed $P_{\text{fix}}$ as a function of $s$ as reported in Fig.  8. Two different behaviors can be observed depending on the value of $s$. The small $s$ region is very well fitted by the Kimura formula (1) with an effective population size given by Eq. (26) and with the same value of $B_{1}=3.5$ used in  plummer2018fixation .

From Fig.  8, it is clear that the behavior of $P_{\text{fix}}$, for large enough $s$, is controlled by a different value of the effective populations size, hereafter referred to as $N_{g}$. In one dimension, following  plummer2018fixation , the effective population size is estimated by considering the scale $\delta$ near to a source in $x_{s}$ where $u_{0}\delta\sim B_{2}2\sqrt{D\mu s}$, with $B_{2}$ another constant of the order of $1$: an initial population in $x\in[x_{s}-\delta,x_{s}+\delta]$ can develop a Fisher genetic wave at speed $2\sqrt{D\mu s}$, which is supported by the velocity field. Only Fisher genetic waves that start in this interval are able to cross the system; this provides an estimate $N_{g}=2\delta(s)\rho_{0}$. In two dimensions, the same argument gives:

[TABLE]

Using Eq. (27), we obtain the curve (black) of Fig.  8, which provides an excellent fit of the numerical simulations.

Finally, since both $N_{\text{eff}}$ of Eq. (24) and $N_{g}$ given by Eq. (27) are proportional to $N=64^{2}N_{0}$ for our simulations, we can easily predict that upon increasing $N_{0}$, the fixation probability will follow the same master curve if plotted as a function of $sN$. To demonstrate this and to validate the quality of our method for large $N_{0}$, we show, in Fig.  9, $P_{\text{fix}}$ as obtained for the same flow as Eq. (25) for $N_{0}=2$ and $N_{0}=10$. The red dot-dashed and black continuous curves obtained using the prescriptions discussed above for small $s$ and large $s$, respectively, provide an excellent fit for the numerical results. Overall, the results discussed in this section extend the ones previously obtained in  plummer2018fixation and demonstrate the validity of our method for population dynamics advected by an external compressible velocity field.

VI VI. Conclusions

In the present work, we developed a numerical method suitable for accurately and efficiently investigating the behavior of population dynamics and genetics under flow. This approach allows for the study of a large number of individuals by, first, implementing the diffusion and advection processes, particle by particle, and afterwards, for each box composing the 2D lattice, performing the birth and competition steps.

In order to test and validate our method, we considered a one-dimensional system. We implemented the FKPP equation, analyzing the algorithm convergence. After that, we applied this method to the heterozygosity and Kimura formula and we found a very good agreement between the theoretical and simulated results. The method we propose does not require any dynamic management of particle positions and has no limitations on the number of individuals for mesh points. Both features imply major simplifications in computer coding, especially for a large number of individuals and for parallel computation. It is worth remarking that for a large number of individuals, one can increase the computational efficiency of our method by directly sampling the binomial distribution in each mesh point along the lines discussed in binomial .

For the 2D system, we retraced the procedural scheme of the one-dimensional system and we investigated the larger system under an advection field composed of two sinks and two sources. Our main result was to find, for the 2D system, a net growth of particles born in proximity of a source, as compared to the individuals at different initial positions.

Many interesting studies can follow up on our work. One of these would be to implement a realistic oceanographic advection field and to understand the population and genetic evolution. Another topic to investigate could be the study of the effect of stochastic fluctuations in antagonist population dynamics and the exploration of the effect of external velocity on the genetic nucleation theory.

VII Acknowledgments

The authors would like to thank David Nelson for useful discussions. The work has been performed under the Project HPC-EUROPA3 (Project No. INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme HPC-LEAP; in particular, the authors gratefully acknowledge the computer resources and technical support provided by SurfSARA.