Invasion patterns in competitive systems

TL;DR

This paper explores how Fokker-Planck diffusion and a nonlinear environmental noise model influence invasion patterns in a spatial Lotka-Volterra system, revealing habitat segregation and coexistence phenomena.

Contribution

It introduces a nonlinear noise model for environmental variability and applies Fokker-Planck diffusion to study invasion dynamics, offering new insights into habitat segregation and coexistence.

Findings

Fokker-Planck diffusion causes habitat segregation of resident and invader.

Nonlinear noise can break segregation, enabling coexistence.

Heterogeneous distributions influence invasion success.

Abstract

Stochastic reaction-diffusion equations are a popular modelling approach for studying interacting populations in a heterogeneous environment under the influence of environmental fluctuations. Although the theoretical basis of alternative models such as Fokker-Planck diffusion is not less convincing, movement of populations is commonly modelled using the diffusion law due to Fick. It is an interesting feature of Fokker-Planck diffusion that for spatially varying diffusion coefficients the stationary solution is not a homogeneous distribution; in contrast to Fickian diffusion. Instead, concentration accumulates in regions of low diffusivity and tends to lower levels for areas of high diffusivity. Thus, we may interpret the stationary distribution of the Fokker-Planck diffusion as a reflection of different levels of habitat quality. Moreover, the most common model for environmental fluctuations, linear multiplicative noise, is based on the assumption that individuals respond independently to stochastic environmental fluctuations. For large population densities the assumption of independence is debatable and the model further implies that noise intensities can increase to arbitrarily high levels. Therefore, instead of the commonly used linear multiplicative noise model, we implement environmental variability by an alternative nonlinear noise term which never exceeds a certain maximum noise intensity. With Fokker-Planck diffusion and the nonlinear noise model replacing the classical approaches we investigate a simple invasive system based on the Lotka-Volterra competition model. We observe that the heterogeneous stationary distribution generated by Fokker-Planck diffusion generally facilitates the formation of segregated habitats of resident and invader. However, this segregation can be broken by nonlinear noise leading to coexistence of resident and invader across the whole spatial domain.

Tables

Table 1. Table 1: Parameter values for the simulations presented in Figures  3 - 9 .

$m$	$n$	$d_1$	$d_2$	a	${\omega }_1$	${\omega }_2$	${\gamma }_1$	${\gamma }_2$	Fig.
1	1	5	25	19	0.1	0.2	1.0	1.0	3
1	1	25	5	19	0.1	0.2	1.0	1.0	4
1	1	25	12.5	19	0.1	0.2	1.0	1.0	5
1	1	25	12.5	19	0.1	0.4	1.0	1.0	6
1	2	25	12.5	19	0.1	0.4	1.0	1.0	7
1	2	25	12.5	0	0.1	0.4	1.0	0.1	8
1	2	25	12.5	19	0.1	0.4	1.0	0.1	9

Equations

$$\frac{\partial X_i(\vec{r},t)}{\partial t}=\underset{\text{reaction}}{\underbrace{f_i(\mathbf{X}(\vec{r},t))}}+\underset{\text{diffusion}}{\underbrace{\nabla \cdot \left(-\mathbf{\mu }(\mathbf{x},t)X_i(\vec{r},t)+\nabla [D(\vec{r},t)X_i(\vec{r},t)]\right)}}+\underset{\text{stochastic}}{\underbrace{g_i(\mathbf{X}(\vec{r},t))\xi (\vec{r},t)}},$$

$$d{\mathbf{X}}_t=\mathbf{\mu }({\mathbf{X}}_t,t)dt+\mathbf{\Sigma }({\mathbf{X}}_t,t)d{\mathbf{W}}_t$$

$$\frac{\partial }{\partial t}p(\mathbf{x},t)=\nabla \cdot \left(-\mathbf{\mu }(\mathbf{x},t)p(\mathbf{x},t)+\nabla [D(\mathbf{x},t)p(\mathbf{x},t)]\right)$$

$$D_{i,j}(\mathbf{x},t)=⟨{\mathbf{\Sigma }}_{i,\cdot },{\mathbf{\Sigma }}_{j,\cdot }⟩.$$

$$\mathbf{\mu }(\mathbf{x},t):=\alpha \nabla D(\mathbf{x},t),\alpha \in \mathbb{R},$$

$$\frac{\partial }{\partial t}p(\mathbf{x},t)=\nabla \cdot \left[D(\mathbf{x},t)\nabla p(\mathbf{x},t)\right]$$

$$\frac{\partial }{\partial t}p(\mathbf{x},t)=\Delta [D(\mathbf{x},t)p(\mathbf{x},t)].$$

$$g_i(\mathbf{X})=\frac{{\omega }_i{X}_i^m}{{\gamma }_i+{\alpha }_i{X}_i^n}.$$

$$\frac{\partial X_1}{\partial t}=$$

$$\frac{\partial X_2}{\partial t}=$$

$$D^{\ast}(x,y)=D_{0}+\left\{\begin{array}[]{lc}a\left(\sin(\sqrt{x^{2}+y^{2}})\right)^{k}&\text{if }\sqrt{x^{2}+y^{2}}<3\pi\,,\\ a\left(\sin(3\pi)\right)^{k}&\text{else\,.}\end{array}\right.$$

$$D_0=1,a=19,k=8$$

$$c_{12}=c_{21}=1.2$$

$${\alpha }_1={\alpha }_2=\mathrm{0.1.}$$

$$h_t=0.02,h_x=h_y=15$$

Full text

1 Introduction

It is needless to give a broad review of the classical publications on spatial and spatiotemporal pattern formation in non-equilibrium nonlinear systems. However, on occasion of Ulrike’s significant birthday, one should remember her seminal contributions to the theory of pattern formation in electrochemical systems from the eighties of last century Ebe:83 ; Feu:84b ; Feu:84a , of course without forgetting all her impressive later work until today. In that mentioned early period of her academic career, one of us (H.M.) had the chance to learn from and to work with Ulrike in Werner Ebeling’s research group at the Sektion Physik of Humboldt-Universität zu Berlin Mal:86 . Later on, a theoretical bridging of electrochemical and ecological diffusive systems was found Mal:91 ; Mal:89b . Ulrike became interested in ecological and environmental dynamics as well, and from time to time we manage to meet and to chat and sometimes even to work, cf. Ben:14 . Also in the academic career of I.S., Ulrike Feudel has played an important role—Ulrike kindly served as an external examiner for his PhD defence. He has fond memories of being quizzed on the Lotka-Volterra competition model which also plays a central role in the study presented here.

Interactions and movements of populations $\mathbf{X}(\vec{r},t)=\{X_{i}(\vec{r},t);~{}i=1,2,\ldots,N\}$ in a heterogeneous and variable environment are often modelled with stochastic reaction-diffusion equations:

[TABLE]

Here, the reaction terms $f_{i}(\mathbf{X}(\vec{r},t))$ describe the interactions between individuals of a population with individuals of the same or a different population. This enables us to represent processes as diverse as transmission of infectious diseases, predator-prey interactions or competition for resources.

The diffusion term is derived from an underlying stochastic model of the movement of individuals. Consider the stochastic differential equation for the position $\mathbf{X}_{t}\in\mathbb{R}^{d}$ of a particle that moves stochastically in $d$-dimensional space:

[TABLE]

The drift coefficients $\bm{\mu}:\mathbb{R}^{d}\times\mathbb{R}^{+}\to\mathbb{R}^{d}$ are account for deterministic movement of the particle, $d\mathbf{W}_{t}$ is a $d$-dimensional Wiener process and the matrix-valued intensity coefficients are $\bm{\Sigma}:\mathbb{R}^{d}\times\mathbb{R}^{+}\to\mathbb{R}^{d\times d}$. An alternative description to the SDE (2) which enables us to calculate the stochastic location $\mathbf{X}_{t}$ of a particle is the probability density $p(\mathbf{x},t)$ for finding the particle at position $\mathbf{x}$ at time $t$. This analogous representation of the system is especially suitable for considering a large population of particles because in this case, the probability density $p(\mathbf{x},t)$ can be interpreted as the fraction of particles that are expected to be found at a location $\mathbf{x}$ at time $t$. The probability density $p(\mathbf{x},t)$ can be shown to satisfy a deterministic partial differential equation (PDE), the Kolmogorov forward or Fokker-Planck equation Kol:31 , (Oku:01, , chapter 5)

[TABLE]

where $\nabla$ is the gradient, $\nabla\cdot$ denotes the divergence operator and $D(\mathbf{x},t)$ is related to the intensities $\bm{\Sigma}(\mathbf{x},t)$ by the standard scalar product $\langle\cdot,\cdot\rangle$:

[TABLE]

Now, choosing

[TABLE]

we obtain several alternative laws for the dynamics of the probability distribution $p(\mathbf{x},t)$ i.e. the collective movement of the population. First, it is important to note that in this model different choices of the parameter $\alpha$ only have an effect for spatially varying diffusion coefficients $D(\mathbf{x},t)$. Second, it is easy to see that the stationary solution of (3) depends on $\alpha$. For $\alpha=1$ an easy calculation shows that (3) reduces to the well-known diffusion law due to Fick Fic:55

[TABLE]

which is the diffusion term most commonly used in modelling applications. From (6) it is clear that the homogeneous distribution $p(\mathbf{x},t)=1$ is a stationary solution. In contrast, for $\alpha=0$ we obtain a law known as Fokker-Planck diffusion Fok:14 ; Pla:17 :

[TABLE]

Here it is obvious that the stationary distribution cannot be homogeneous for spatially varying diffusion coefficients $D(\mathbf{x},t)$. The reason can be understood by regarding (2). By choosing $\alpha=0$ the movement of individual particles is purely stochastic but nevertheless for spatially inhomogeneous coefficients, the movement is biased towards directions of larger intensities $\bm{\Sigma}(\mathbf{x},t)$. Similarly, we see that in the situation of Fick’s law ($\alpha=1$), the drift term $-\nabla D(\mathbf{x},t)$ opposes the gradient of $D(\mathbf{x},t)$ and in this way seems to exactly balance the bias of the movement towards larger diffusivities.

For many physicochemical systems, Fick’s law is the model of choice due to the fact that the particles move in a purely passive way which is consistent with a flux opposed to the concentration gradient. But for biological populations whose individuals are able to actively influence their direction of movement, there is no reason to restrict ourselves to models based on Fick’s law. Instead, it seems more appropriate to start from a general model for stochastic movement such as (3). As explained in more detail in Ben:16b the drift term (5) is not only a phenomenological description but can be interpreted as the ability of an individual to “sense” environmental conditions over a distance increasing with $0<\alpha\leq 1$ and “choose” its direction of movement accordingly Fro:13 . Potapov et al. Pot:14a considered $\alpha=2$. However, regardless of possible interpretations of the stochastic movement of individuals underlying a particular diffusion term, the most important qualitative feature of the alternative models to Fick’s law presented here is the fact that the stationary distributions of the populations are in most situations inhomogeneous. Indeed, Bengfort et al. Ben:16b showed that Fokker-Planck diffusion leads to pattern formation for situations where this would not be expected for Fickian diffusion.

Whereas Bengfort et al. considered a wide range of deterministic models, most recently we have also investigated the combined influence of Fokker-Planck diffusion and stochastic environmental fluctuations Ben:16c . As illustrated in (1), by adding a stochastic term, fluctuations in environmental conditions such as temperature, nutrient availability etc. over time and/or in space, can be incorporated without having to represent each source of environmental variability individually. Usually, multiplicative Gaussian noise which is uncorrelated both in space as well as in time is used, i.e. the standard normally distributed random variable $\xi_{i}(\mathbf{x},t)$ is multiplied by the population density $X_{i}$ so that we have $g_{i}(\mathbf{X}(\vec{r},t))=X_{i}(\vec{r},t)$. As explained in more detail in Siek:16a , this particular choice of $g_{i}$ can be related to an individual-based model, the branching process in a random environment (BPRE). In a BPRE, the amount of offspring produced by each individual is modulated by a stochastic process that represents environmental fluctuations. For large population numbers, the BPRE can be approximated by a stochastic differential equation with one term that accounts for demographic stochasticity as well as a multiplicative noise term $X_{i}(\vec{r},t)\xi_{i}(\mathbf{x},t)$ representing the influence of environmental fluctuations. Thus, the commonly used multiplicative noise model for environmental stochasticity can be derived from a so-called diffusion approximation of a BPRE for large populations where demographic stochasticity is neglected.

That the effect of environmental fluctuations scales with the population number is essentially due to the fact that environmental stochasticity is assumed to affect each individual independently. However, the more the population density increases, the more likely it seems that individuals are located so close to each other that instead of responding independently to stochastic perturbations they are similarly affected. From these considerations it is expected that instead of increasing linearly with the population number, the intensity of the environmental fluctuations saturates, which suggests a model of the following form

[TABLE]

which in slightly more general form has recently been proposed by Siekmann and Malchow Siek:16a . For $m=n$ the noise intensity $g_{ii}$ monotonically tends to a maximal noise intensity $\omega_{ii}/\alpha_{ii}$. The half-saturation constant $K:=(\gamma_{ii}/\alpha_{ii})^{1/n}$ is the population density $X_{i}$ at which half of the maximal noise intensity is reached. For $m<n$ the noise intensity decreases for large population numbers $X_{i}$. Regardless of the mechanistic interpretation given above, the most important qualitative difference of (8) to the multiplicative model is the fact that the noise intensity is bounded.

The purpose of this study is to examine in more detail the combined effect of Fokker-Planck diffusion (7) and nonlinear noise (8). This is motivated by our most recent work where we found that by varying the standard approaches for modelling movement of populations Ben:16b or environmental fluctuations Siek:16a , respectively, a wide range of interesting effects could be observed, even in well-known classical models such as the Lotka-Volterra competition model. Here, we use the Lotka-Volterra model for representing an invasive species and investigate which effect Fokker-Planck diffusion and nonlinear noise have on the success of the invader. This study continues Bengfort et al. Ben:16c where we combined Fokker-Planck diffusion and the classical multiplicative environmental noise model with the Lotka-Volterra model presented here.

2 The stochastic competition-diffusion model

The dynamics of a resident species $X_{1}$ and an invader $X_{2}$ is described by

[TABLE]

The spatial dependency of the resident’s diffusivity is chosen as

[TABLE]

Here, the parameter $k$ is an even number witch controls the steepness of $D^{\ast}$. Throughout this paper we will use the parameters

[TABLE]

This functional form of the diffusivity $D^{\ast}(x,y)$, see Figure 1 for a plot, is meant to represent the resident’s varying levels of preference for different areas of the spatial domain. The coefficient $D^{\ast}(x,y)$ can be regarded as being inversely proportional to the resident’s preference for a particular location $(x,y)$. Namely, the lower $D^{\ast}(x,y)$, the lower the tendency to leave $(x,y)$ which can be interpreted as a high level of preference. How these preferences for different parts of the habitat affect the stationary distribution of the resident is fundamentally different for Fickian diffusion and Fokker-Planck diffusion. For Fickian diffusion the spatially heterogeneous diffusion coefficient $D^{\ast}(x,y)$ only affects the transient dynamics of the resident because the stationary distribution is always homogeneous, regardless of the particular functional form of $D^{\ast}(x,y)$. In contrast, for Fokker-Planck diffusion, the stationary solution is approximately inversely proportional to $D^{\ast}(x,y)$ which provides us with a simple model for a fragmented habitat that mimics the resident’s levels of preference. In the absence of the invader $X_{2}$, the resident $X_{1}$ tends to the distribution shown in the right panel of Figure 1—note that for the initial condition we have always set the resident’s population to zero within a square with a side length of 50 length units in order to mimic the onset of a biological invasion.

We ensure that resident and invader do not differ in competitive strength by letting the competition parameters coincide

[TABLE]

Because both $c_{12}$ and $c_{21}$ exceed unity, the system is in the bistable parameter range i.e. in the absence of diffusion or noise the competitor with the larger initial density will drive its opponent to extinction. For spatially extended systems, Malchow et al. showed that survival depends on the ratio of the diffusion coefficients of invader and resident—in general, the competitor with the larger diffusion coefficient prevails Mal:11 . Because the spatially varying diffusion coefficient $d_{1}\cdot D^{\ast}(x,y)$ of the resident is larger than the constant diffusion coefficient $d_{2}$ of the invader in some areas and smaller in others, it is expected from Malchow et al. Mal:11 that the spatial domain becomes segregated. Whereas in some regions invasion is successful due to relatively low diffusivity of the resident, other areas act as barriers for invasion where the resident’s diffusivity is relatively high. We will vary $d_{1}$ and $d_{2}$ in order to explore the effect of different ratios of the local diffusion coefficients.

We will investigate this situation under the influence of nonlinear noise (8), thereby extending our previous study where we applied the standard model of multiplicative noise Ben:16c . Thus, throughout the manuscript the noise intensities $g_{1}$ and $g_{2}$ are of the form (8) for which we always choose

[TABLE]

The parameters $\omega_{i}$ and $\gamma_{i}$ will be varied in order to study different dependencies of the noise intensity on the population densities $X_{i}$.

The initial condition for a spatial grid of $200\times 200$ grid cells is indicated in Figures 1 and 2, the invader is set to zero and the resident is initialised with the spatially heterogeneous stationary distribution in the whole spatial domain except for a patch of $50\times 50$ grid cells in the upper left. Here, the resident is set to zero whereas the invader is set to its carrying capacity 1. An exception is one simulation (see Figures 8 and 9) where we compare Fick’s law with Fokker-Planck diffusion where the resident is initialised with a spatially homogeneous distribution.

We numerically solve (9), (10) using an alternating direction implicit (ADI) scheme for efficiently implementing the Crank-Nicholson method Cra:47 as described previously Ben:16c . The stochastic terms for which we use the Stratonovich interpretations are numerically integrated with the derivative-free Milstein method Klo:99 ; Mil:95 as explained in Siek:16a . The temporal and spatial step sizes (in non-dimensional units) are usually

[TABLE]

except for Figures 8 and 9 where a smaller temporal step width of $h_{t}=0.002$ was required.

3 Numerical simulations and results

With the parameters described in the previous section we now investigate the success of invasion depending on various choices of the scaling parameters $d_{1}$ and $d_{2}$ of the resident’s and invader’s diffusivities as well as different parameters $\omega_{i}$ and $\gamma_{i}$ for the noise. As mentioned above, an initial patch of invaders of $50\times 50$ grid cells in size, i.e., 6.25% of the total model area, is placed in a corner of the spatial domain, cf. Figures 1 and 2. Parameters for the simulations presented here are summarised in Table 1.

The influence of the ratio of $d_{1}$ and $d_{2}$ is demonstrated in Figures 3-5. For Figures 3 and 4, the only difference between both parameter sets is that $d_{1}$ and $d_{2}$ are swapped. Consistent with Malchow et al. Mal:11 mentioned in the previous section if $d_{1}$ is low compared to $d_{2}$, the invasion is successful (Figure 3) whereas it fails for $d_{1}$ much larger than $d_{2}$ (Figure 4), although the invader initially seems to be able to enter the realm of the resident ($t=100$). Interestingly, Figure 3 shows that even in the case of successful invasion, the invader is not able to overcome all barriers created by large diffusivities of the resident so that the resident prevails in the centre of the spatial domain. The result is a segregation of the spatial domain in one habitat dominated by the invader and one habitat dominated by the resident.

For Figure 5, the diffusivity $d_{2}$ has been increased. In contrast to Figure 4 the invader not only manages to proceed into the territory of the resident but succeeds in establishing itself along a ring where the stationary population of the resident is at a low level. Very slowly, the invader manages to displace the resident from this ring without being able to occupy any other region in the spatial domain. Thus, like in Figure 3 we again end up with spatially segregated habitats but here the invader is confined between two regions occupied by the resident.

Increasing the noise intensity of the invader increases its ability to cross invasion barriers caused by large values of the diffusion coefficient $D^{\ast}(x,y)$ of the resident. Consistently, in Figure 6 we observe a similar situation as in our initial simulation (Figure 3) with low diffusivity of the resident in comparison with the invader. Again, the invader is able to conquer a habitat that extends to the boundary of the spatial domain instead of remaining constrained to a ring.

In Figure 7 we see that this result is hardly unchanged if we replace the monotonously increasing noise by terms that initially increase and then decrease for large densities. This is achieved by choosing a larger exponent in the denominator of (8) i.e. $m<n$. The resulting pattern strongly resembles Figure 6, also note that the dynamics even develops on a similar time scale.

For this form of the noise term even coexistence of invader and resident within the whole spatial domain rather than the emergence of separate habitats is possible. The fact that the coexistence pattern can be observed both for Fickian as well as Fokker-Planck diffusion shows that this is an effect mostly mediated by the noise term. In comparison with the previous simulation (Figure 7) we shift the half-saturation constant of the invader to lower population densities by decreasing $\gamma_{2}$. This has the effect of a steeper increase of the noise term (8). In Figure 8 we show that for Fickian diffusion both populations quickly mix. The initial invader patch develops into a mixed front consisting of both resident and invader that replaces the area formerly occupied only by the resident.

The behaviour is similar for Fokker-Planck diffusion but here it can be noted that the preference of the resident for certain regions of the spatial domain is apparent—within rings less favoured by the resident where its diffusivity has maximum values the density of the invader is noticeably higher (Figure 9).

4 Conclusions

We have studied a biological invasion based on a spatio-temporal Lotka-Volterra competition model under the influence of stochastic environmental fluctuations. We found surprisingly rich dynamics after replacing the two standard models for movement of populations and environmental variability. Instead of Fickian diffusion, dispersal of the resident was modelled by the so-called Fokker-Planck law of diffusion. The most important difference of Fokker-Planck diffusion to classical Fickian diffusion is that for spatially varying diffusion coefficients, the stationary distribution is heterogeneous which can be interpreted as the result of a population’s preference for different spatial regions within their habitat. Moreover, the Fokker-Planck law of diffusion has a mechanistic basis that is not less well-suited for modelling movement of biological populations than Fickian diffusion. Our results show that Fokker-Planck diffusion facilitates invasion because the invader may manage to establish itself in spatial areas that are less favoured by the resident.

Similarly, we have considered an alternative for the classical model for environmental fluctuations based on noise intensities that increase linearly with the population densities. One interpretation based on branching processes in a random environment (BPRE) is that each individual is affected independently by environmental stochasticity, an assumption that is clearly debatable for large population densities. Replacing the linearly increasing noise intensities by a model for which noise intensities saturate or even decrease with increasing population densities accounts for the fact that individuals rather than responding independently are similarly affected by environmental perturbations. The most important qualitative difference of this model is the fact that the noise intensity does not exceed a certain upper bound even for large population densities. Similar to the classical model, increasing the noise intensity of an invader usually facilitates invasion by increasing the speed of invasion.

Whereas for many parameter sets the model predicts a segregation of the spatial domain into separate habitats for resident and invader, the new noise model produces an even more interesting effect. Mediated by the noise, resident and competitor are even able to coexist in a mixed habitat that extends across the whole spatial domain. It could be shown that this effect is indeed primarily due to noise because it occurs for both Fokker-Planck as well as Fickian diffusion.