A Stable Finite-Difference Scheme for Population Growth and Diffusion on a Map

We describe a general Godunov-type splitting for numerical simulations of the Fisher–Kolmogorov–Petrovski–Piskunov growth and diffusion equation on a world map with Neumann boundary conditions. The procedure is semi-implicit, hence quite stable. Our principal application for this solver is modeling human population dispersal over geographical maps with changing paleovegetation and paleoclimate in the late Pleistocene. As a proxy for carrying capacity we use Net Primary Productivity (NPP) to predict times for human arrival in the Americas.

Fisher [10] studied the problem-via a growth-diffusion equation-of an analogous but one-dimensional situation: the propagation of an advantageous genetic mutation within an already-present population, situated along a coast line. Kolmogorov, Petrovskii and Piskunov [11] were more general; in particular, their analysis treated the two-dimensional case. Such a model (called Fisher/KPP in the following) was applied to the dispersal and growth of a population by Skellam [12], and serves as an important control for designing and validating other more complex spatiotemporal population models [5]. Coupling population dynamics with models of large-scale changes in continental topography, climate, and ecosystem productivity is necessary to understand the role of environmental constraints on patterns of genetic, phenetic, and cultural variation among human populations [5].
Here we present a stable and efficient finite-difference solver for the Fisher/KPP equation on the 2-D domains of geographical maps, and show how it can be extended to include environmental fluctuations. In a brief outline of our paper, we will: review the derivation of the Fisher/KPP equation (Section The Fisher/KPP equation); develop finite-difference schemes in 1 and 2 dimensions for variable environmental carrying capacity K (Section Numerical a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 methods and splitting) depending on both space and time Kðx; tÞ (Section Space-time dependent capacity maps); and show our application of this technique to the out-of-Africa dispersal of Homo sapiens by using net primal productivity (NPP) [1] as a proxy for K (Section World-wide hominin dispersal).

The Fisher/KPP equation
An intuitive way to get the Fisher/KPP equation [10][11][12] is as follows. A current j of particles (e.g., individuals) moving across an interface located at x is proportional to the gradient of the population density p (called "Fickian diffusion") j ¼ À c rp: The rate of change of p is then given by the mass balance equation [13], which for Fickian diffusion reads @p @t À c r 2 p ¼ r: If ρ = 0, this is the heat equation when c = D/2 and D is the diffusion coefficient. For lack of a better model, we assume c is a constant (Young & Bettinger [4]). In Subsection Parameter optimization, we optimize on D = 2c. The source term ρ is usually modeled by a logistic growth function, r ¼ lpð1 À p=KÞ, and gives the Fisher/KPP equation where K is called carrying capacity and λ is the growth rate. In Eq (1) we wrote K as a constant.
Since ρ is a local density, it can be space, x, and time, t, dependent. To get relative density information, all we need is that there is an upper limit " K, in which case 0 p " K. If Kðx; tÞ is evolving, locally the solution p Kðx; tÞ " K. As in Table 1, " K can be scaled to unity, although obviously only the relative density p= " K can be computed. A scaled version for constant K is given in Eq (2), while space and space-time dependent Kðx; tÞ versions are given in the Section Space-time dependent capacity maps. Fluctuations in climate produce environmental changes in vegetation, sea levels, opening/closing of land bridges, waxing/waning of ice sheets, and perturbations to habitable areas in general. Thus, time-dependent environments compel us to deal with space-and time-dependent Kðx; tÞ (see Subsection Time interpolation of maps), Section Space-time dependent capacity maps, Eq (9). Fisher and KPP were particularly interested in the traveling wave case, p(x, t) = f(x 0 + vt). It was shown in [14,15] that asymptotically (large t), the speed is v ¼ 2 ffiffiffiffi ffi lc p ¼ ffiffiffiffiffiffiffiffi 2rD p in our notation. By the rescalings show in Table 1, for constant K the Fisher/KPP Eq (1) is written where the only sensible solutions have 0 u 1. The initial distribution u(x, 0) = u 0 (x) must be defined for all x on the habitable regions of the map. In Murray ([14], eq. (11.17)) our growth coefficient λ is called r and c is denoted by D, whereas in Young and Bettinger [4] the growth coefficient is R and the diffusion coefficient is K. We use r to denote radius in Section Numerical methods and splitting, and k for the CFL parameter, so λ and D are used here and that next section. Elsewhere, λ = r is used interchangeably. These inputs to our code are given in units of yr −1 and km 2 /yr, respectively.

Numerical methods and splitting
In one dimension, Eq (2) can be solved using the MatLab function pdepe, or if the system is two-dimensional but rotationally symmetric, pdepe can again be used with the radial part of the Laplace operator in cylindrical coordinates, requiring only that one sets a pdepe parameter m=1. What is important is that because MatLab function pdepe is robust, it is a valuable verification tool for testing our more general solver when comparisons can be made.

The finite-difference scheme
Since the map on which we will be working is a pixelized plane, an obvious method uses finite differences. First, however, let us examine the 1-D case for Eq (2). In this situation, the second order derivative becomes a differencing operator in matrix form acting on the vector where the matrix A is If h is the time step, the Courant-Friedrichs-Lewy (CFL) parameter [16] is An explicit integrator for Eq (2) would require k < 1/4 [16,17]. In our case, because boundary conditions are so irregular on a map and there are strong spacial variations in K, we are uninterested in a method of higher order than second because stability is more important [18].
Using this notation, the lowest order approximation is Euler's method which estimates the next step u(t + h) by which should be considered a vector equation in u(t) = {u j (t), j = 1, n}. The logistic terms, which are diagonal, should be taken to mean ((1−u)u) j = (1−u j )u j for j = 1, 2, . . ., n. Euler's method is both low-accuracy and generally unusable if it is used alone over too many steps. This method is in principle conditionally stable if the step size is forced to be small enough. A step size which is too small, however, actually drives errors up, not down: e.g., C. Wm. Gear's Fig. 1 (6): Eq (7) can be solved as a linear system, because the matrix, 1 À k 2 A À h 2 ð1 À u E Þ, on the left hand side is explicit, as is the right hand side. That is, this matrix and the right hand side contain only old data, namely just information from the previous step, u(t). Euler estimate u E Eq (5) is an explicit one step computation using old data u(t). Significant advantages are: the matrix on the left hand side is tridiagonal with constants on the sub/super-diagonals, and the diagonal terms are O(1) strong. The procedure Eq (8) is only linearly stable but we find empirically that it gives good results when compared to pdepe when comparisons to this MatLab function are appropriate.

Space-time dependent capacity maps
Fluctuations in climate produce environmental changes in vegetation, sea levels, opening/closing of land bridges, waxing/waning of ice sheets, and perturbations to habitable areas in general. Thus, changing environments compel us to develop a procedure Eq (9) for the case when Kðx; tÞ depends on both space and time.
Since Eq (8) is basically the trapezoidal method (see section 5.3 in [17]), the modification for a space-time dependent Kðx; tÞ is as follows: where we have suppressed the x dependence of Kðx; tÞ for simplicity of notation. In Eq (9), the operators A x and A y are the same finite difference operators as in Eq (3) but for directions x and y, respectively. For simulations on a lattice, u ij (t) = u(x 0 + (i−1)Δx, y 0 + (j−1)Δy, t), where 1 i N x , 1 j N y and Δx = Δy, the following gives the action of the A x , A y operators: A compression scheme and code outline given in Appendix Map segmentation show that only a maximum of one row or column (i.e., max(N x , N y )) of storage is needed for u ? and u ?? .
We use the simple regularization of each version of 0 u K in Eq (9) at every x−point, namely For example: use Eq (9b) to compute u E , then set u E max(0, u E ), followed by u E min ðKðx; t þ hÞ; u E Þ. Intermediate values u ? and u ?? are constrained similarly.

Fisher/KPP on geographical maps
In order to solve Fisher/KPP on a geographical map, one segments the map into x-and ydirection pieces, imposing boundary conditions at their endpoints; the solver acts on each segment independently, changing the directions x ! y Eqs (9a) to (9b), then y ! x Eqs (9b) to (9d) [17] of integration as in Eq (8). This approach also lends itself to efficient parallelization. In the Map segmentation appendix we illustrate the scheme with an illustrative MatLab code sample.

Boundary conditions for Fisher/KPP
Two versions for boundary conditions were implemented and tested: (1) zero solution Dirichlet, and (2) zero net flux Neumann conditions. Both variations require some explanation.

Dirichlet BCs
In this case, a border of zero pixels were added to the land masses. In our Godunov method, each direction (X, respectively Y) segments acquire two zero pixel end points, and the resulting n seg + 2 one dimensional difference equations are solved, but only the n seg land mass pixel values of the solution are updated. For example, in the Appendix Map segmentation, an X direction segment passing through equatorial central Africa has n seg pixels, each of which also has a corresponding space/time dependent carrying capacity Kðx; tÞ. Two zero pixels are added, one at each end, but Kðx; tÞ for those additional points are not relevant, nor used. Single pixel land mass segments thus become 3 pixels, where only the middle point is updated. The differencing operator (3) becomes very simple with three elements/row (1, -2, 1). Simulations with high resolution, say N x = 720 and N y = 360, show that the coastal region populations look too small. This is because of smoothing between zero (water) and the second adjacent pixel away. Our ancestors fished, so we concluded Dirichlet BCs seem inappropriate.

Neumann BCs
In consequence, we use zero net flux, @u/@x = 0 or @u/@y = 0 into the sea conditions. This is more awkward and two issues must be dealt with: • Derivative information (u x , u y to O(Δx 2 )) must be available, or approximated. From LeVeque [17], we know that when n seg > 3, we may use @u=@x % 3 2 uðxÞ À 2uðx þ DxÞ À þ 1 2 uðx þ 2DxÞÞ=Dx as a one-sided (to the right, here) difference approximation. If n seg 3, however, things are clumsier.
• Furthermore, because our Godunov splitting advances two one half steps of the heat equation per time step, there is the question of whether the solver simulates these half-steps in a properly posed way. It is well known (again [17]) that the heat equation with Neumann boundary conditions can be an ill-posed problem. This is easy to see: u t ¼ 1 2 u xx with u x = 0 at the end points remains invariant to the shift u ! u + C, for any constant C. Thus, in that case, the solution cannot be unique. Uniqueness must be imposed by the logistics term in (2). That term is not invariant under the shift, so uniqueness is assured. Only single half time-steps for the heat equation steps, Eqs (9a) and (9d), are taken.
If n seg > 3, the resulting linear equations to be solved are of the form Eq (9c).
where RHS 1 = RHS n seg = 0. Elements u E,j denote the j−th elements of the Euler estimate vectors u E . Note that this system has bandwidth 5, not 3. However, the upper and lower 2nd sub/ super-diagonals contain only one element (1/3). As LeVeque [17] points out, this only slightly disturbs the tridiagonal structure. To be able to still use a tridiagonal solver, we modify Eq (10), by using explicit Heun rule estimates for u 3 and u n seg −2 and move those to the right hand side: that is, and remove the first and n seg -th row elements 1/3 in Eq (10). The u 1 and u nseg elements can be updated after the tridiagonal system solution to correct for the explicit estimates, if desired. We find no discernible differences doing so. Three special cases remain: n seg = 1, 2, 3.
• If n seg = 1, only the logistics term can be used: set A y = 0 in (11b) and (11c) and update the single element u ?? 1 . This must be kept non-negative, however.
• Surprisingly, the n seg = 2 is the most awkward. Eq (11) use the 2nd order differencing operator on the old (previous time step) solutions, and likewise A y . The resulting solutions of Eqs (9a) and (9d) contain limited derivative information, hence the Neumann boundary conditions require u ?
2 , and u 1 (t + h) = u 2 (t + h) and all these must be non-negative. Similar restrictions are also made on the Euler estimate in (11b).
• When n seg = 3 the resulting 3-vector equations for u ? u ?? , and u(t + h) all take a form A couple of simple manipulations show u 1 = u 2 = u 3 , which says that there is no derivative information if Neumann boundary conditions are imposed on this n seg = 3 case. The solution is u 2 = RHS 2 /(2S 1 + S 2 ), and thus we get the other two. They must also be kept non-negative.

World-wide hominin dispersal
We now focus on our principal application using the methods presented above: the worldwide dispersal of Homo sapiens out of Africa.

Capacity maps
Our construction of a time-dependent K uses Net Primary Productivity (NPP) as a proxy [1]. The Miami model [20] was originally formulated in 1972 to estimate NPP (in (grams of carbon in dry organic matter)/m 2 /day) from annual temperature and rainfall [21]. In order to compute our NPP maps, we obtained the temperature and precipitation data from simulations by the BRIDGE program [22] organized at the University of Bristol [23][24][25][26]. The simulation data that we use are computed on a 96 × 73 grid, which we interpolate to size 720 × 360 and convert to NPP maps [27] by applying the formulas given in [20]. World-wide NPP data are difficult to obtain, so our Miami model-like maps are somewhat rough. Via testing, we found our Godunov solver to be quite robust with respect to abrupt x −steps in the carrying capacity Kðx; tÞ.

Time interpolation of maps
We assembled 61 NPP maps, from 120 kya to 1 kya. These are in 4 kilo-year steps for the first 10; 2 ky steps for the next 21; then 1 ky steps for the remainder. The time stepper in our Godunov scheme has no information about continuity between these discrete NPP maps, so an interpolation scheme needs to be used.
If we have available some estimates for the carrying capacities K L , K H at times t L , t H , one possible estimate for KðtÞ at times t L t t H in between is a homotopy, where a sigmoid function 0 = S(t L ) S(t) S(t H ) = 1 smoothly interpolates between lower and higher time frames. There are many choices available, such as that used in [1]. We use the following variant.
Start with the classical sigmoid which is zero at z = −1 and unity at z = +1. The −1 < z < +1 interval is not what we want, but the following t 7 ! z transformation permits several variants: where ΔT = t H − t L . Notice that S(z(t L )) = 0 and S(z(t H )) = 1. The exponent ν in Eq (12) gives some freedom in choosing a particular form for z for almost any ν > 0. If ν < 1/2, d 2 S/dt 2 has more than two sign changes, so ν ! 1/2 is preferable. With choice ν = 1/2, the interpolant is nearly a straight line: see Fig 2. However, its turn-up at t = t L and turn-down at t = t H numerically resemble very quick derivative changes. So, we choose ν = 1 which also makes z(t) time scale independent. At both ends, all derivatives of S(z(t)) in t smoothly vanish. We also have the forward/backward symmetry S(z(t H − t)) = 1−S(z(t)) for t L t t H , as does [1].

Parameter optimization
Other tasks remain: optimize the r = λ and D = 2c parameters in Eq (1); choose our starting value T s , the start time for our simulation; and choose a threshold value for arrival. For example, a migration wave reaching some fraction of the local carrying capacity at a studied site would constitute arrival. Roughly, the distance from Eritrea/Djibouti to any site, arriving at time t A , looks like where a, b are unknown fractions of the distance. The terms consist of a nearly constant speed ($ ffiffiffiffiffi ffi rD p ) traveling wave, and a diffusion ffiffiffiffiffiffiffi Dt A p , respectively. This is a scaling argument using the units D * km 2 /yr and r * yr −1 . A fraction ccFRAC % 1/2 of the carrying capacity would mostly probe the first term, a ffiffiffiffiffi ffi rD p t A . Conversely, a tiny value of ccFRAC weighs too much of the diffusion tail b ffiffiffiffiffiffiffi Dt A p . The right-hand frame of Fig 1 shows that the wave consists of a steep traveling wave with extended diffusion tails. We found that arrival times are relatively insensitive to values of ccFRAC *0.2. Larger values can produce much later arrival times, as in Table 2 and its missing error bars. Hence, our choice is ccFRAC = 0.1: When the incoming population reaches 10 percent of the local carrying capacity, it has arrived. For this ccFRAC choice, the differences between NPP and FEP (Food Extraction Potential [28]) arrivals did not seem significant. The opening of the Bering Straight "bridge" occurs before any of our simulation start-times, T s , so arrival times in the Americas become almost linearly increasing for T s ! 45kya. Hence, T s = 45kya value is a reasonable lower limit. Our NPP data and masks are later than 120 kya, after the beginning of the Tarantian period. A Bering Strait crossing was possible until the Holocene period began around 10 kya. Thus migration data would permit a transit until closure. Arrival in northeastern India is known to have happened before 40kya, so T s cannot be that small. Thus to limit our optimization space, we choose T s = 45kya. For this choice of T s , we still have to choose r, D. Some estimates were known [4,5] to be approximately D % 200km 2 /yr and r % 2.0 × 10 −3 yr −1 . With the hope that we might be more precise, we optimize on the following RMS parameter, where arch MC = earliest estimated population of Meadowcroft, PA; similarly, arch FC = earliest known population of Fell's Cave, near the tip of Cape Horn, Isla Grande, Tierra del Fuego  doi:10.1371/journal.pone.0167514.g003 [29,30]. Our simulation values sim MC and sim FC depend on the diffusion/growth parameters r, D, the arrival threshold ccFRAC, and starting time T s .

Population dispersal
An initial population at t = T s = 45 kya at an Eritrean is shown in Fig 4. The integration units are scaled following Table 1 such that λ = 2.04 Á 10 −3 y −1 and D = 2c = 190 km 2 /y (consistent with those used in [4]). The solver is then run using time frame K maps described above, down to 1 kya before present. The remaining plots (Figs 5, 6, 7, 8 and 9) display the resulting population dispersal simulation on space-time dependent K maps. Using population parameters optimized in Subsection Parameter optimization gives results consistent with the literature. The gross features of the late (45-50 kya) out-of-Africa dispersal of Homo sapiens are reproduced [32], e.g. the colonization of Western Europe before 40 kya and that of South America before 14 kya [31].  table of Table 2 are also quite reasonable considering the uncertainties involved in our crude NPP maps. The 2 nd column of Table 2 also indicates archaeological uncertainties.

Conclusions
In this paper, we present a semi-implicit Godunov scheme for the Fisher/KPP equation with a variable carrying capacity K. In one dimension, the expected traveling wave [10,11] develops as shown in Fig 1. As in the one-dimensional situation, a traveling wave also develops as expected [11] in 2-D. Plots in Subsection Population dispersal clearly show these evolving waves.
In the Section Space-time dependent capacity maps we describe our procedure for handling x and time dependent Kðx; tÞ, Eq (9).
Finally, we apply our scheme to our principal objective of population dynamics: the out-of-Africa dispersal of Homo sapiens. On the Mercator projected world map, by using vegetation net primal productivity (NPP) as a proxy for carrying capacity, we get Table 2 which shows that interpolating in time between discrete KÀ maps yields stable and reasonable results. These results show ancestor arrival in NE India before 40 kya, then a crossing of the Bering Strait before 10 kya. Multiple routes into South Asia [33] are also evident. Honesty requires that we admit our size (56km) 2 pixels do not resolve the two crossing points at Bab al Mandab and Sinai adequately. Additionally, we make no attempt to model coastlines effectively. Neither simple boats nor easier passages along beaches are included in our model: only Neumann boundary conditions determine such movements here. Better modeling of coastlines and more detailed NPP maps would improve our results.
The core computations performed by our solver are independent tridiagonal solutions, which can be easily parallelized to deal with larger grids (e.g. [34], p. 116). In order to improve numerical performance, in the Appendix Map segmentation, we discuss a compressed storage scheme to integrate the Fisher/KPP equation on a projected world map. About 71% of the earth's surface is water, so this compressed storage reduces computational work by the same amount.

Provenance
For this paper, the simulations were run on either a Mac Mini, 2.3 GHz Intel Core i7, or a 2.

Map segmentation appendix
Our solver on a geographical map uses a map outline, i.e., a rectangular grid with 1's in habitable regions, and 0's in the water, see  Scheme for Population Growth and Diffusion on a Map independent rows, resp. columns, indexed by starts/ends of habitable segments. Independent columns, i = 1. . .NX in Fig 10 will have nysegs(i) of habitable segments, whose start and end positions are ystart_seg(k) and yend_seg(k), respectively, where k = 1. . .nysegs(i). Likewise, for j = 1. . .NY rows, each have nxsegs(j) also with start/ end positions. A 100 × 50 example, Fig 10, shows row 26 has 4 segments of varying size. Column 87 has 5 segments.
Each time a new carrying capacity map is loaded (up to 61 total) [1,22] a suitable mask is computed from K L and K H , then resegmentation is performed: see Subsection Time interpolation of maps. Resegmentation facilitates time changing coastlines.