In this section, the heterogeneity of the landscape is introduced by assuming that its profile can be written as (Psi (x) = a + psi (x)), where (psi (x)) represents the spatial variations of the environment around a reference level a.
The results that we will present were obtained through theoretical and numerical techniques. The theoretical approach is based on the mode linear stability analysis discussed in the previous section. Numerical integration of Eq. (3), starting from a homogeneous state plus a random perturbation, was performed following an explicit forward-time-centered-space scheme, with boundary conditions suitably chosen for each case (see Supplementary Information for details).
Refuge
As a paradigm of a heterogeneous environment with sharp borders, we first consider that the spatial variations around the reference level a are given by
$$begin{aligned} psi (x) = – A[1- Theta (L/2 -|x|)] ,, end{aligned}$$
(6)
where (Theta) is the Heaviside step function. With (A>0), it represents a refuge of size L with growth rate a immersed in a less viable environment with growth rate (a-A). In a laboratory situation, this can be constructed by means of a mask delimiting a region that protects organisms from some harmful agent, for instance, shielding bacteria from UV radiation26. In natural environments, this type of localized disturbance appears due to changes in the geographical and local climate conditions27, or even engineered by other species38.
In Homogeneous landscapes section, we have seen that the uniform distribution is intrinsically stable when (q ge 0). In contrast, when there are heterogeneities in (Psi (x)), spatial structures can emerge even if (qge 0), as illustrated in Fig. 3 for the case (D=0.01).
Stationary population density (rho _s) vs. x in a refuge. This heterogeneous environment is defined by Eq. (6), with (a=b=1), (A=10^{-3}) and (L=10). The vertical lines indicate the refuge boundaries. We used the kernel (gamma _q(x)), with (q=0.1) and (ell =2), and two different values of D. Symbols are results from numerical integration of Eq. (3) under periodic boundary conditions, and solid lines from the small-A approximation given by Eq. (8), in excellent agreement with the exact numerical solution. Recall that, in a homogeneous environment, no oscillations appear for (q ge 0).
In the limit of weak heterogeneity, i.e., under the condition (|psi (x)|/a ll 1), we obtain an approximate analytical solution assuming that the steady solution of Eq. (3) can be expressed in terms of a small deviation (varepsilon _s(x)) around the homogeneous state (rho _0=a/b). Then, we substitute (rho _s(x)=rho _0+varepsilon _s(x)) into the stationary form of Eq. (3), discard terms of order equal or higher than (mathcal{O}(varepsilon ^2, Avarepsilon ,A^2)), and Fourier transform, obtaining
$$begin{aligned} tilde{varepsilon }_s(k) = dfrac{ rho _0 tilde{psi }(k)}{-lambda (k)},, end{aligned}$$
(7)
where (lambda (k)) was already defined in Eq. (5) and (tilde{psi }(k)) is the Fourier transform of the small fluctuations in the landscape quality, which for the case of Eq. (6) is (tilde{psi }(k)= A[2sin (Lk/2)/k -2pi delta (k)]).
Finally, assuming that (lambda (k^star )<0), the steady density distribution is given by
$$begin{aligned} rho _s(x) ,=, rho _0 +varepsilon _s(x) ,=, rho _0 + mathcal{F}^{-1}Bigl (dfrac{ rho _0 tilde{psi }(k)}{-lambda (k)}Bigr ) ,, end{aligned}$$
(8)
where the inverse Fourier transform (mathcal{F}^{-1}) must be numerically computed in general. For small heterogeneity, Eq. (8) is in very good agreement with the exact numerical solution obtained by integration of the dynamics Eq. (3), as can be seen in Fig. 3. Notice the two different profiles, depending on the diffusion coefficient D: one gently following the landscape heterogeneity and the other strongly oscillatory.
For small D, the induced oscillations display two evident characteristics, which depend on (tilde{gamma }_q): a well-defined wavenumber and an amplitude that decays with the distance from the interface at (x=pm L/2) (highlighted by dashed vertical lines in Fig. 3). We will see in the next section how the characteristics of the oscillations reflect the details of the kernel (gamma _q).
Semi-infinite habitat
Since oscillations are induced by changes in the landscape, it is worth focusing, from now on, on one of the interfaces between a more viable region and a less viable one. Moreover, we assume a refuge much larger than the oscillations wavelength, sufficient to follow over several cycles the structure originated at the interface. To do that, we consider a semi-infinite habitat defined by
$$begin{aligned} psi (x) = -ATheta (-x) ,, end{aligned}$$
(9)
where for convenience the interface was shifted to (x=0), such that the low-quality region is at (x<0). As an additional feature, we consider that the harmful conditions are very strong, that is, (A rightarrow infty). The purpose is twofold, on the one hand, it allows to test the robustness of the results beyond the small-A approximation, on the other, it allows a simplification as follows. When (A gg a), (rho) is very small in the unfavorable region, then the nonlinear competition term in Eq. (3) can be neglected, leading to a steady distribution that decays exponentially from the interface as (rho (x<0) sim exp [sqrt{(A-a)/D},x]). Thus, in the limit (Arightarrow infty), we have (rho (x <0,t)= 0). In addition, the semi-infinite habitat is simulated by the interval [0, L], where L ((=100) in our simulations) is large enough in comparison to oscillation length-scales. Then, far away from the interface, we set (rho (x ge L,t)=rho _0). This is the setting used to produce Fig. 1b, by numerical integration of Eq. (3).
As sketched in Fig. 4, for each steady distribution attained at long times, we measure the wavelength, from which we obtain the wavenumber (bar{k}), and the decay length (bar{x}), by observing that the envelope of the oscillations decays as (exp (-x/bar{x})).
Characterization of stationary profiles. Long-time solutions approach a stationary state characterized by the wavelength (2pi /bar{k}) and decay length (bar{x}). The slope of the straight black line is indicated in the figure. This example was obtained from numerical integration of Eq. (3), assuming a semi-infinite habitat, with parameters (D= 0.003), (gamma _q(x)) with (ell =2) and (q=-0.5).
The stationary spatial structures that emerge for (x>0) can be classified into the three types depicted in Fig. 1: sustained oscillations (lilac line, with (bar{k}>0) and (bar{x}rightarrow infty)); decaying oscillations (orange line, with (bar{k}>0) and finite (bar{x})); exponential decay (gray line (bar{k}=0) and finite (bar{x})). In the case of Fig. 1b, these three types appear when D changes. We also systematically varied the shape parameter q to construct the phase diagram in the plane (q-D) presented in Fig. 5a.
Phase diagram and characteristics of the stationary profiles as a function of diffusion coefficient D and q, in the semi-infinite habitat. We used the kernel (gamma _q(x)), with (ell =2). (a) Phase diagram in the (q-D) plane, and cuts at (b) (D=10^{-3}), (c) (q=-0.5) (d) (q=0.5). The remaining parameters are (a=b=1). In diagram (a), for each point in the grid, the type of regime was determined based on the values of (2pi /bar{k}) and (bar{x}) that characterize the solutions of Eq. (3): sustained oscillations ((bar{k}>0) and (bar{x}rightarrow infty), lilac), decaying oscillations ((bar{k}>0) and finite (bar{x}), orange), and pure exponential decay ((bar{k}=0) and finite (bar{x}), gray). In (a), the dashed and dotted lines correspond to ({k_i}=0) and ({k_r}=0), respectively, where ({k_r}) and ({k_i}) are the real and imaginary parts of the zeros of (lambda(k)), with the smallest positive imaginary part. In (b)-(d), symbols correspond to measurements of numerical profiles, according to Fig. 4, and solid lines correspond to the prediction in Eq. (10) (theoretical 1). Thin dashed lines correspond to the harmonic estimate (theoretical 2) given by Eq. (14).
To perform a theoretical prediction of (bar{k}) and (bar{x}), within the linear approximation, we consider that these oscillation parameters should be related to the poles of the integrand (mathrm{e}^{ikx} tilde{psi }(k) /[-lambda (k)]) in the expression for the inverse Fourier transform that provides the solution, according to Eq. (8). As far as the external field (psi (x)) does not introduce non-trivial poles, like in the case of a Heaviside step function ((tilde{psi }(k) sim 1/k)), only the zeros of the complex extension of (lambda (k)) matter. The dominant (more slowly decaying mode) is given by the complex poles (k=pm k_r + i k_i) ( (k_r>0)) with minimal positive imaginary part that, except for amplitude and phase constants, will approximately provide patterns of the form (mathrm{e}^{-k_i x} cos (k_r x)), allowing the identifications
$$begin{aligned} bar{k} = k_rquad text { and }quad 1/bar{x}= k_i. end{aligned}$$
(10)
This theoretical prediction42 is in very good agreement with the results of numerical simulations, as shown in Fig. 5, explaining the observed regimes.
Moreover, the modes that persist beyond the interface have relatively small amplitudes, so that the system response is approximately linear in this region.
Lastly, recall that this analysis assumes mode stability ((lambda (k)<0)). When (lambda (k^star )>0), the system is intrinsically unstable, with the poles having null imaginary part (lying on the real axis). Nevertheless, the initially fastest growing mode, given by the maximum of (lambda (k)), tends to remain the dominant one in the long term41, yielding (bar{k} simeq k^star) for the sustained oscillations ((bar{x} rightarrow infty)).
In order to obtain further insights, it is useful to consider the response function (tilde{R}(k)) that, from Eq. (7), is
$$begin{aligned} tilde{R}(k) equiv frac{|tilde{varepsilon }_s(k)|^2}{|tilde{psi }(k)|^2} = frac{rho _0^2}{lambda ^{2}(k)} ,. end{aligned}$$
(11)
Despite missing some of the dynamical information contained in the phase of (lambda (k)), it can provide a more direct estimation of the observed parameters than through calculation of the poles. In order to perform this estimation, we resort to the response function of a driven damped linear oscillator43 described by the equation (varepsilon _H”(x)+2zeta k_0varepsilon _H'(x)+k_0^2varepsilon _H(x)= f(x)). We have
$$begin{aligned} tilde{R}_H(k) equiv frac{|tilde{varepsilon }_H(k)|^2}{|tilde{f}(k)|^2}= frac{1}{|lambda _H(k)|^{2}} = frac{1}{(k^2-k_0^2)^2+4zeta ^2k_0^2k^2}, , end{aligned}$$
(12)
with (-lambda _H(k) = -k^2 + i2zeta k_0 k + k_0^2), whose zeros (poles of (1/lambda _H(k))) are (k= pm k_r + i k_i = k_0 (pm sqrt{1-zeta ^2}+izeta )), where (k_0) is the natural mode and (zeta) is the damping coefficient. Note that, under a step forcing (f(x)=k_0^2 Theta (x)), which simulates our present setting, those poles carry the essential information of the damped-oscillation solution, given by (tilde{varepsilon }_H(k) = tilde{f}(k)/[-lambda _H(k)]), where (tilde{f}(k)=k_0^2(pi delta (k) -i/k)). In the underdamped case ((zeta <1)), this solution is explicitly given by
$$begin{aligned} varepsilon _H(x) = left[ 1 – frac{k_0}{kappa } e^{-x/xi } sin (kappa x +phi )right] Theta (x) , , end{aligned}$$
(13)
where (kappa = k_0sqrt{1 -zeta ^2}) ((=k_r)), (xi = 1/(zeta k_0)) ((=1/k_i)), and the phase constant (phi =tan ^{-1}(xi kappa )). The solution for the overdamped case emerges for (zeta >1), when the zeros of (lambda (k)) are pure imaginary with (k_i=k_0(zeta pm sqrt{zeta ^2-1})). The connection between the poles of (tilde{R}_H(k)) and the dynamic solution of the driven harmonic oscillator is possible because, as previously discussed, (tilde{f}) does not introduce relevant poles, and the forced solution has a form similar to the unforced one.
The harmonic model is, in fact, the minimal model sharing characteristics with our observed structures, and the correspondence between Eqs. (12) and (13) will allow to estimate the oscillation features. In the limit of small (zeta), (tilde{R}_H(k)) has a sharp peak, characterized by a large quality factor (Q equiv k^star /Delta k), where (Delta k) is the bandwidth at half-height of (tilde{R}(k)) around (k^star)43. First, we see that the position of the peak of (tilde{R}_H) approximately gives the oscillation mode (kappa), according to (k^star = k_0sqrt{1-2zeta ^2} = kappa + mathcal{O}(zeta ^2)). Second, the bandwidth is related to the decay-length through (Delta k = 2/bar{x} + mathcal{O}(zeta ^2))<a data-track="click" data-track-action="reference anchor" data-track-label="link" data-test="citation-ref" aria-label="Reference 44" title="If $$k_-{ <}k_+$$ k – 44.
Putting all together, as long as (tilde{R}(k)) resembles the bell-shaped form of (tilde{R}_H(k)), we can use the following estimates, which are correct for the harmonic case to first order in (zeta):
$$begin{aligned} bar{k} simeq underset{k}{text {arg max}},(tilde{R}) equiv k^star quad text { and }quad bar{x} simeq frac{2}{Delta k} ,. end{aligned}$$
(14)
The expression for (bar{x}) is also valid in the overdamped limit (large (zeta) in the harmonic model), in which case the maximum is located at (k^star =0).
The adequacy of the harmonic framework as an approximation to the response function of our model, (tilde{R}(k)), is illustrated in Fig. 6. In the case (D=2times 10^{-1}), the harmonic response is able to emulate (tilde{R}(k)). Then, if the harmonic approximation holds, one expects that the estimates given by Eq. (14) should work for the population dynamics case. In fact, they do work, as we will see below. Differently, when (D=2times 10^{-4}), (tilde{R}(k)) does not follow the harmonic shape, it is multipeaked and the dominant mode observed in the simulations is not given by the absolute maximum.
Comparison of (tilde{R}(k)) with the harmonic response (tilde{R}_H(k)), both normalized to their maximal values. (tilde{R}(k)) of our model, given by Eq. (11) (solid lines) and harmonic response (tilde{R}_H(k)), given by Eq. (12) (dashed lines), where the values of (k_0) and (zeta) were obtained by fitting Eq. (12) to (tilde{R}(k)). In all cases, (q=0.5), (ell =2) and two different values of D shown in the legend were considered. Notice that for (D=2times 10^{-1}), the response can be described by the harmonic approximation. For (D=2times 10^{-4}), the response is multipeaked, indicating that the harmonic approximation fails. In fact, the dominant mode observed in the simulations is not given by the absolute maximum, but by the small hump at (ksimeq 2.1), as predicted by the analysis of complex poles.
In Fig. 5, we compare the values of (bar{k}) and (bar{x}) extracted from the numerical solutions of Eq. (3) with those estimated by Eq. (14) (dashed lines) and, more accurately, with those predicted from the poles of (tilde{R}(k)) (solid lines), which perfectly follow the numerical results. The harmonic estimates are shown in the full abscissa ranges, as a reference, even in regions where the approximation is not expected to hold, because discrepancies give an idea of the departure from the harmonic response.
Figure 5c shows outcomes for a fixed (q<0) ((q=-0.5)), corresponding to a vertical cut in the diagram of Fig. 5a. Sustained oscillations (i.e., (bar{x} rightarrow 0)) can emerge for (q<0), when diffusion is weak, namely, for (D<D_c simeq 0.0025) (lilac colored region), where (D_c) is obtained from (lambda (k^star )=0). When D increases beyond this critical value, oscillations are damped with a finite characteristic length (bar{x}). For even larger values of D, oscillations completely disappear ((bar{k}rightarrow 0)). Note that the comparison between numerics and harmonic theory (symbols vs. dashed lines) is good close to the pattern transition point (D_c), where the response peak is sharp (large Q). Despite the lack of agreement for larger D, the harmonic approximation qualitatively works with a shift of the transition from attenuated oscillations to exponential decay.
Figure 5d (which corresponds to vertical cut at (q=0.5) in the diagram of Fig. 5a) shows the corresponding results for a fixed (q>0) ((q=0.5)), which is characterized by the absence of sustained patterns.
Above (Dsimeq 0.02), the response (tilde{R}(k)) is unimodal, a bell-shaped curve that resembles the harmonic response, as in the case (D=0.2) (black lines) in Fig. 6, producing a good agreement between harmonic and numerical results, despite being far from the large-Q limit. However, for smaller values of D, the profile is multi-peaked, and not even (k^star) predicts the observed mode, indicating that the harmonic approximation does not hold, as for (D=2times 10^4) (red lines) in Fig. 6. In this regime, it is crucial to analyze the response function in terms of complex poles in order to extract the dominant mode and its decay.
Figure 5b displays (bar{k}) and (bar{x}) as a function of q, for a fixed value of the diffusion coefficient ((D=10^{-3})), corresponding to a horizontal cut in Fig. 5a. Recall that, the smaller the value of q, the more confined is the interaction (thus, the larger is (bar{x})). For (q < q_c approx -0.093) there are sustained oscillations ((bar{x} rightarrow infty)). Above (q_c), oscillations decay, which is indicated by the transition of (1/bar{x}) from null to finite values. Again, near this transition, the harmonic approximation works well, but, far from the critical point, it fails, as noticed above (qsimeq 0.2), where there is a strong mismatch between the main mode given by the harmonic approximation and the numerical one. Also in this case, a small hump in the response function represents the dominant mode, as predicted by the analysis of the complex poles of (tilde{R}).
Two-dimensional landscapes
In this section, we show results of simulations for relevant 2D scenarios, verifying that the picture of induced oscillations described up to now for 1D also holds in 2D.
Long-time spatial distribution in 2D. Simulated scenarios: (a) a circular region (with radius 5 a.u., highlighted with a black dashed boundary) where the growth rate is positive, a (in a strong negative background (a-A)); (b) a circular region (with radius 2.5 a.u., highlighted with a black dashed boundary) where the growth is strongly negative (a-A) (while outside, it is positive, a); (c) four regions with negative growth rates (a-A) (in a positive background, a); (d) time-independent random landscape (where each spatial cell is assigned a growth rate uniformly distributed in [0.5a, 1.5a]). In all cases the interaction kernel is (gamma _q), with (ell =2) and (q=0.5), (D=10^{-3}), (a=b=1) and (A=10) . Colors show the deviation from the homogeneous state (rho (x,t)-rho _0) (where (rho _0 = 1) for the chosen values of the parameters). For numerical integration, a pseudo-spectral method45 was used with (Delta x = 0.2) and (Delta t = 10^{-3}).
Snapshots of simulations for different 2D landscapes are presented in Fig. 7: a refuge (a), a defect (b), multiple defects (c) and spatial randomness (d) where many spatial scales are present. It is worth remarking that, in 2D, for the kernel (gamma _q), patterns only appear in homogeneous landscapes if (q<q_c simeq 0.25) (i.e., if (lambda (k^star )>0)). Thus, in all the cases of Fig. 7 (using (q=0.5)) we would not find oscillations if the landscape were homogeneous. In Fig. 7, we see that for 2D the same picture as in 1D is found: decaying oscillations appear near landscape disturbances with a clear wavenumber and decay length. The linear response approach presented in Semi-infinite habitat section can straightforwardly be extended to 2D. Figure 7a–c shows the case in which defects either increase or decrease the population growth rate. This can be promoted by ecosystem engineers such as termites38. Figure 7d shows a case where the landscape is random (in space, but time-independent). This situation, investigated in many previous studies20,36, produces a pattern that is noisy but has a dominant wavelength, which is related to (ell). Furthermore, although there is not a clear identification of decay length from pattern observation, the linear theory would allow one to estimate the characteristic spatial correlation length from the width of the Fourier spectrum.
Source: Ecology - nature.com
