Landscape-induced spatial oscillations in population dynamics

We study the effect that disturbances in the ecological landscape exert on the spatial distribution of a population that evolves according to the nonlocal FKPP equation. Using both numerical and analytical techniques, we characterize, as a function of the interaction kernel, the three types of stationary profiles that can develop near abrupt spatial variations in the environmental conditions vital for population growth: sustained oscillations, decaying oscillations and exponential relaxation towards a flat profile. Through the mapping between the features of the induced wrinkles and the shape of the interaction kernel, we discuss how heterogeneities can reveal information that would be hidden in a flat landscape.

The numerical solution of the one-dimensional generalized FKPP equation presented in Eq. (3) of the main text was obtained by means of a finite-difference forward-time-centered-space (FTCS) scheme, implemented in C language.
For the spatial grid, we consider equally spaced mesh points x i = i∆x, with integer i and grid space ∆x. At each mesh point x i , the density at time t j , ρ i j = ρ(x i ,t j ), evolves according to an explicit scheme with fixed time step ∆t, such that t j = j∆t with integer j. For the spatial discretization of the second derivative, we used the centered form , and for the integration of the nonlocal competition term we used the trapezoidal rule. Lastly, a fourth-order Runge-Kutta approximation for the new values ρ j+1 i was used 1 , which improves the stability domain in comparison with the Euler method. Typical values used for spatial and temporal discretization, respectively, are ∆x = 0.05 and ∆t ∈ [10 −2 , 10 −4 ] (depending on D), leading to an estimated relative error smaller than 10 −4 .
For the initial preparation of the system, we applied a small random perturbation around the homogeneous solution ρ 0 = a/b, drawing a random number ξ i uniform in (−ε, ε), with ε ρ 0 , for each mesh point x i , such that, ρ i 0 = ρ 0 + ξ i . The boundary conditions for each configuration are described below (in Section Boundary conditions) and in Section Stationarity condition, we show the stop criterion used to determine stationarity.
The two-dimensional simulations in Fig. 7 of the main text were performed using the pseudospectral algorithm of Ref.

Boundary conditions
Different boundary conditions were adopted along the main text: (i) For the homogeneous landscape ( Fig. 1a) and for (finite) refuge case (Fig. 3), we used periodic boundary conditions. (ii) For the semi-infinite habitat (e.g., in Fig. S1, with growth rate a in the region x ≥ 0, and with growth rate a − A, otherwise), the integration was performed in a grid with −L ≤ x ≤ L, using L = 100, much larger than oscillation length-scales, under the constraints ρ(x ≤ −L) = (a − A)/b and ρ(x ≥ L) = ρ 0 = a/b. was performed in a grid with 0 ≤ x ≤ L, using L = 100, under the conditions ρ(x < 0) = 0 and ρ(x ≥ L) = ρ 0 = a/b. This choice is justified by the fact that in the limit A → ∞, the density outside the refuge vanishes, as shown in Fig. S1. Figure S2a shows the population density ρ(x,t) vs. x at different instants of time t, obtained from the integration of Eq. (3) for a semi-infinite refuge with strong harmful conditions outside (A → ∞). As time passes, the profile progressively attains a stationary form, ρ(x, ∞). This limiting value can be estimated by noting that the relative difference (discrepancy) between the profiles at different instants decays exponentially with time. In Fig. S2b, the discrepancy |1 − ρ(x i ,t)/ρ(x i , ∞)|, for selected mesh points x i is displayed, showing exponential convergence. The final simulation time was chosen such that the discrepancy is smaller than 10 −4 for all the mesh points in the interval of interest.
In Fig. S3, we show the mode growth rate, λ (k), for three values of α. And in Fig. S4a, we show the phase diagram. In Fig. S4b, we characterize the profiles as a function of parameter 2 − α, for D = 10 −3 . All these results qualitatively resemble those produced in the main text for kernel γ q (Fig. 2 and 5 of the main text).