Revealing new dynamical patterns in a reaction–diffusion model with cyclic competition via a novel computational framework

1. Introduction

The formation of patterns and travelling waves in chemical, physical, biological and game-theoretical models described by systems of advection–reaction–diffusion equations are widely researched phenomena discussed in an immense number of publications. The pioneering works in this area were published as early as the late 1930s with the works [1,2] on wave propagation and Turing’s demonstration of pattern-forming mechanisms in the early 1950s [3]. A nice introduction to the current state of the art can be found in [4–7]. To this day, the research area continues to experience fast development: new types of patterns and waves have recently been reported in a variety of increasingly complex models. For example, a number of exotic patterns, such as envelope and multi-envelope quasi-solitons [8], have been observed in systems with cross-diffusion (i.e. where the diffusivity matrix is non-diagonal). Similarly, novel patterns have also been found in reaction–diffusion systems with non-local terms, i.e. in the case where the reaction terms depend not only on the local species densities but also on the species distribution in some neighbourhood, described by an integral term [6,9]. Furthermore, models of interactions on growing domains or in the case when the spatial parameters are variable can result in the emergence of new types of patterns with different properties from those on fixed domains or with homogeneous parameters [10,11].

Surprisingly enough, despite the current trend towards investigating the behaviour of ever-more sophisticated reaction–diffusion systems (including density-dependent diffusion coefficients, cross-diffusion, time delay, non-local terms, etc.), we still lack a fundamental understanding of the possible types of dynamical regimes in well-known reaction–diffusion models, including some mainstays of introductory textbooks on mathematical biology.

A notable example is the Lotka–Volterra competition model. This model is well known in the literature with numerous applications in mathematical biology and areas of game theory such as voting models [12–17]. In particular, this system is often considered as a paradigm for biodiversity modelling. New types of patterns have been recently demonstrated in this system with spatially homogeneous diffusion coefficients [13,15,17]. This includes, for example, patchy invasion (the spread of a species via the formation and propagation of chaotic patches without a smooth population front), which was originally believed to occur only in predator–prey- or inhibitor–activator-type models [18]. Despite this, the Lotka–Volterra system still remains poorly understood, especially when the interacting species diffuse at different rates. Here, we show the existence of several new dynamical patterns, related to the spread of travelling waves, which have been missed in the literature so far, and which may have important biological applications. In particular, we demonstrate spreading patterns exhibiting complex regular spatial structure which have not been observed so far in reaction–diffusion models with a non-transitive competition, such as the Lotka–Volterra cyclic model.

Another important gap in our knowledge about patterns and waves in reaction–diffusion models is that most existing results have been obtained for either one or two dimensions in space. The simple reason for this is that exploring interactions in three dimensions presents a challenge from both the analytical and the numerical point of view. This is rather a nuisance because many applications of these models, modelling microbiological interactions in water bodies, various medical applications as well as chemical interactions within a substantial volume, are inherently three dimensional. In particular, an important question is how the dimension of the system’s spatial domain influences species persistence and biodiversity [19,20].

A thorough investigation of the Lotka–Volterra model with complex spatio-temporal patterns in two and three space dimensions is made possible here by applying a novel adaptive numerical method, based on a finite element method (FEM) coupled with a reliable a posteriori error estimator. Our choice of method is motivated by the observation that solutions to this model (and many others) feature large patches of relative homogeneity in which a single species dominates, separated by narrow travelling fronts where interactions occur. These narrow fronts can only be resolved by a suitably fine computational mesh, although using such a fine mesh across the large domains and long time scales required for the full system dynamics to develop can be prohibitively expensive, a difficulty which is amplified as the number of spatial dimensions increases. As this resolution is not required in large areas where a single species dominates, the Lotka–Volterra system is a perfect target for mesh adaptivity, tightly focusing the computational effort in the areas of interaction and sparsely deploying it elsewhere.

The algorithm we employ is based on mathematically rigorous error estimates which are computable at each time step because they only depend on ‘known’ quantities such as the numerical solution and problem data. A posteriori error estimation for both stationary and dynamic linear equations is now relatively well understood [21–23]. This approach is extended here by developing energy-norm a posteriori error bounds for semilinear systems of parabolic equations for which the nonlinear reaction terms satisfy suitable growth conditions, ensuring that the framework we present also includes many other similar reaction–diffusion-type systems typically encountered as models of biological and ecological phenomena. For this class of equations, discretized by FEMs, we prove new a posteriori upper bounds for the error of the numerical method. The analysis is based on the elliptic reconstruction technique of [23] and exploits, in an a posteriori fashion, an argument from [24] to bound the nonlinear terms; see also [25] for a similar argument applied to a single reaction–diffusion equation whose solution blows up in finite time. Based on the a posteriori error estimator, we devise a local error indicator which is used to drive the computational mesh adaptation algorithm used in all of the numerical simulations presented here. The standard second-order FEM is judiciously combined with the second-order Crank–Nicolson/Adams–Bashforth implicit–explicit (IMEX) discretization in time, allowing us to treat the linear part of the differential operator implicitly, while the nonlinear reaction term is treated explicitly. The resulting scheme, therefore, allows us to efficiently explore the spatio-temporal patterns arising in the Lotka–Volterra cyclic competition model in both two and three spatial dimensions. We demonstrate its effectiveness by revealing several novel spatio-temporal patterns in the model in two and three spatial dimensions, and explore their properties.

The paper is organized as follows. Section 2 introduces the cyclic competition model with diffusion, while §3 provides the details of the novel numerical method used in our simulations. Some results from our simulations are presented in §4, demonstrating the new dynamical patterns in the cases of two and three spatial dimensions. In §5, we present the proof of a reliable a posteriori error bound for a general class of semilinear parabolic problems, and discuss how it is used to drive mesh adaptivity. A concluding discussion of the significance of our results from the biological and computational viewpoints is given in §6.

2. Cyclic competition in a reaction–diffusion model

The spatio-temporal interactions of three competing species are described using a reaction–diffusion scheme based on the three-species Lotka–Volterra competition model, as proposed by May & Leonard [26] and Mimura & Tohma [13]. The unscaled system in Ω×[0,T], with $\mathit{\Omega }\subset {\mathbb{R}}^d$, where d=2 or 3, reads

where u_i≡u_i(t,x) is the density of species i at time t and location x; a_i is the intrinsic growth rate of species u_i, and each coefficient α_i,j represents the limiting effect that the presence of species u_j has on species u_i (the term α_i,i describes the self-limitation of the population). The diffusion coefficients D_i describe the dispersal rate/mobility of each species. All model parameters are assumed to be positive.

The model is completed by homogeneous Neumann boundary conditions on ∂Ω and the initial conditions

for various different choices of functions ψ_i, i=1,2,3.

As the model contains a large number of parameters, it is convenient to scale each species’ density by introducing ${\tilde{\mathit{u}}}_i={\alpha }_{i,i}{\mathit{u}}_i$, and then rescale both time and space to obtain the simplified three-species competition system,

Note that we can assume without loss of generality that ϵ₂,ϵ₃≤1 because we can always choose to scale the system to the largest diffusion coefficient, and consider the corresponding species as the first one.

The Lotka–Volterra competition model with diffusion has been studied in a number of papers, where it was shown that it can possess complex patterns of dynamics [12–17,27,28]. The outcome of interactions between the species is strongly affected by the hierarchy of the competition structure, as dictated by examining each pairwise interaction. One of the most interesting scenarios includes a cyclic competition structure, in which (roughly speaking) species 1 outcompetes species 2, species 2 outcompetes species 3 and species 3 outcompetes species 1. Such cyclic dominance is analogous to the popular game of ‘rock–paper–scissors’ and may be observed throughout Nature. Well-known examples include competition between side-blotched lizards [29], coral reef invertebrates [30], and various yeast [31] and bacterial strains [32]. The same model also arises in non-biological situations such as many-player Prisoner’s Dilemma games [33] or some types of voter models [34].

Formalizing the above characterization of cyclic competition can be tricky, although here we will follow the definition given by Adamson & Morozov [15]. This is based on considering the outcomes of pairwise interactions in an unbounded one-dimensional spatial domain (i.e. in the absence of a third species), starting from initial conditions such that the species densities at positive and negative infinity are equal to the carrying capacity for one species and zero for the other. Cyclic competition is then said to occur if the direction of the resulting travelling waves follows the cyclic order 1>2>3>1. For example, the domain occupied by species 2 at its carrying capacity level should eventually be replaced by a spreading wave of species 1. This generic definition of cyclic dominance allows two main types of local dynamics [15,17]. In classical cyclic competition, the phase portrait of each pairwise interaction should involve only one stable steady state corresponding to the presence of the stronger competitor at its carrying capacity. In this case, adding a spatial dimension to the local interaction does not reverse the outcome of the competition because the corresponding travelling wave will be directed from the domain occupied by the stronger competitor to that of the weaker competitor [35]. The mathematical conditions for this to occur are: α_i,i+1≤1 and α_i+1,i>1. Under conditional cyclic competition, on the other hand, some local pairwise interactions can be bistable: both of the axial steady states (corresponding to the carrying capacities of one species and zero density for the other) are locally stable and the final outcome of the local competition will depend on the initial conditions. Mathematically, assuming bistability occurs for interactions between species 1 and 3, this means that α_1,3>1 and α_3,1>1. Adding a spatial dimension into the model with conditional cyclic competition should preserve the displacement order 1>2>3>1 as in the classical cyclic competition. However, the conditional cyclic competition involves some constraints on the diffusion coefficients [36]: using arbitrary diffusion coefficients will not guarantee the cyclic dominance. In this paper, we will explore patterns corresponding to both the classical and the conditional cyclic competition scenarios.

Interestingly, for the model of cyclic competition without space, the long-term coexistence of all species is not possible [26], while coexistence can occur when space is considered [15,17]. Furthermore, the dimensionality of the spatial domain plays a role in the long-term persistence of all species [15], and various scenarios of coexistence in the model (2.1) have been reported in the case of two spatial dimensions. Along with well-known patterns of dynamics such as spiral waves and target waves and their combinations, new regimes have been recently demonstrated including the spread of irregular patches, zip-shaped waves and wedge-shaped travelling fronts [15,17]. It has also been observed that the species mobilities ϵ₂,ϵ₃≤1 can play a key role in ensuring species persistence and provide the system with a rich variety of dynamical regimes [15]. In this work, we will demonstrate novel patterns of dynamics for the cyclic Lotka–Volterra reaction–diffusion model with non-equal diffusion coefficients, which have not been reported before.

Finally, the solutions of system (2.1) largely depend upon the habitat and initial conditions. In all cases, we consider the spatial domain Ω to be a d-dimensional cube, d=2,3, and we explore several different types of initial conditions which imitate different scenarios of biological invasion. In the two-dimensional case, we use the following conditions: (i) from a point near the centre of the habitat, we extend three boundary lines at an angle 2π/3 from each other, thus dividing the habitat into three ‘triangular’ subdomains, and we start with one and only one species present in each of these subdomains at its carrying capacity density (we shall refer to this as the ‘triangular’ initial condition for brevity); (ii) in a habitat initially occupied by a single species at its carrying capacity, we introduce the other two species locally in circular subdomains such that the centres of the two circles do not coincide. In the three-dimensional case, we explore using initial conditions formed by dividing the whole domain into eight equal boxes and consider that the boxes are initially occupied by only one species at the density corresponding to its carrying capacity. In all cases, the initial conditions are smoothed by replacing the zero density/carrying capacity interfaces with a smooth narrow transition layer.

3. Numerical scheme

The numerical scheme we use to approximate solutions to the model (2.1) is based on a second-order C⁰-conforming FEM in space, coupled with a second-order Crank–Nicolson/Adams–Bashforth IMEX time discretization, initialized by a single step of an explicit Euler method. The spatial mesh is automatically locally coarsened and refined on certain time steps in response to an estimate of the local contribution to the discretization error, evaluated using the error indicator derived in §5. Such an adaptive algorithm can make the numerical scheme significantly more efficient because solutions of this problem contain a large range of length scales, from the small moving features in the chaotic interaction regions to the large areas where a single species dominates. Consequently, using a uniform mesh sufficiently fine to resolve all solution features would be prohibitively expensive in terms of memory consumption and computation time, even on high-performance hardware. This issue is exacerbated in three spatial dimensions, and the problem quickly becomes completely intractable with a uniform fine mesh. Instead, the adaptive method, outlined in algorithm 1, refines the mesh to add extra resolution where the estimated error is high, and coarsens where the estimated error is low and less resolution is acceptable.

We begin by rewriting the model (2.1) in the following more general form. Let $m\in \mathbb{N}$ be the number of variables (species) and let ε be a smooth m×m diffusion tensor which is positive definite on Ω. Then, we consider the general reaction–diffusion system given by

The nonlinear interaction term $f:{\mathbb{R}}^m\to {\mathbb{R}}^m$ is assumed to satisfy the general Lipschitz-style growth condition that there exists 0≤γ<2 when d=2, or $0\le \gamma \le \frac{4}{3}$ when d=3, such that, for any $\mathit{v},\mathit{w}\in {\mathbb{R}}^m$, we have with $\mid \cdot \mid$ denoting the Euclidean norm on ${\mathbb{R}}^m$.

Let the matrix $A\in {\mathbb{R}}^{m\times m}$ contain the interaction parameters, with A_i,j=α_i,j. The interaction term of the model (2.1) may be seen to satisfy this condition with γ=1 by writing it as

where $\mathbf{\text{1}}\in {\mathbb{R}}^m$ denotes the vector with all entries equal to 1, ° denotes the Hadamard product between tensors and for $r\in {\mathbb{R}}^m$ we use $\hat{r}$ to denote the m×m matrix with the elements of the vector r on its main diagonal. Note that the given system has a unique solution which follows from the general theory of PDEs [37].

We shall use the notation H¹(Ω) to denote the Hilbertian Sobolev space of functions with square-integrable first derivatives, namely

Here, L²(Ω) denotes the usual Lebesgue space of real-valued functions which are square-integrable over Ω, and (L²(Ω))^l, $l\in \mathbb{N}$, the space of vector fields with components in L²(Ω). We refer the reader to, for example, [38,39] for more on Sobolev and Bochner function spaces.

Multiplying by a suitable test function and integrating by parts, the problem may be written in weak form as: find u∈H¹(0,T;(H¹(Ω))^m) such that

Here, ∇u denotes the Jacobian matrix of u, and the L² inner product between matrix-valued functions A(x) and B(x) is defined to be $(A,B):={\int }_{\mathit{\Omega }}A:B\mathrm{d}\mathit{x}$, where A:B denotes the Frobenius product.

The spatially discrete FEM for solving this problem can then be written as follows. Let V _h be a conforming finite element space that is a finite dimensional space of functions which are continuous over Ω and piecewise polynomial of degree k with respect to a mesh ${\mathcal{T}}_h$ covering Ω; see, for example, the monograph [40]. Further, let (V _h)^m denote the space of vector fields with components in V _h. The spatially discrete FEM reads: find u_h∈H¹(0,T;(V _h)^m) such that

In particular, we select piecewise quadratic finite elements for our computations.

To produce a practical method, the time derivative in the spatially discrete problem presented above must also be discretized. Here we use IMEX schemes, allowing us to treat the linear part of the differential operator implicitly, while the nonlinear reaction term is treated explicitly; see [41,42] and references therein for more details. Also, the finite element space on each time step will be different, in general, although we refrain for expressing this dependence explicitly in the notation for accessibility. In particular, adopting the second-order Crank–Nicolson/Adams–Bashforth IMEX time discretization, we obtain the fully discrete problem: for 1≤n≤N−1, find uⁿ⁺¹∈(V _h)^m such that

where the function uⁿ is the finite element approximation to the solution u at time moment tⁿ. This choice of time-stepping scheme avoids the need to solve a system of nonlinear equations at each time step while still providing second-order accuracy in time.

The above numerical method was implemented using the open source deal.II C++ finite element library [43]. An essential feature of this library is its support for high-order spatial discretizations and adaptive meshes in both two and three spatial dimensions.

4. New patterns of spread in the cyclic competition model

We start our investigation of (2.1) with the case d=2. Our goal is to report some previously unknown transient regimes in which the area where all three species coexist spreads into areas where only one species is present. We first consider the classical cyclic competition scenario with model parameters given by a₂=1,a₃=1;α_1,1=1; α_1,2=1; α_1,3=2; α_2,1=2; α_2,2=1; α_2,3=1; α_3,1=1; α_3,2=2; α_3,3=1, for which our numerical simulations revealed new dynamical patterns of travelling population waves with regular spatial structures. The geometric shape of the regular structures in the wake of the front largely depends on the diffusion coefficients.

For small values of the diffusion coefficient of the second species, the resulting patterns have a triangular droplet-like shape which is shown in figure 1. This corresponds to the diffusion coefficients ϵ₂=0.1 and ϵ₃=0.6. In this and all subsequent figures, we use the following colour coding: dark blue indicates regions dominated by species u₁ (i.e. u₁≥u₂,u₃), light yellow indicates regions dominated by species u₂ and mid-green indicates regions dominated by species u₃. To keep the figures as clear as possible in the presence of the narrow interfaces between the species, which are very small relative to the size of the domain, the colouring used for each species does not vary with its magnitude. We note, however, that u₁,u₂,u₃∈[0,1] for all points in time and space. While this property is not strictly preserved by the numerical method, the impact on the numerical results is small.

• Download figure
• Open in new tab
• Download PowerPoint

We use the symmetrical triangular initial conditions described in §2. During the initial stages of the simulation, a spiral tip is formed. This subsequently transforms into a droplet-like domain with a sharp wedge and a complex structure inside. The sharp wedge moves at a constant speed towards the lower left corner of the domain. At the tip of the wedge, all three species are present. A regular structure emerges at the front of the wedge, and both move with the same speed. This regular structure is formed of small droplet-like units which are periodic in the direction orthogonal to the bisector of the wedge. The number of small droplets in the wedge is constantly growing through time. However, moving progressively from the wedge to the centre of the spreading domain in which the species coexist, the spatial structure becomes more irregular, thus the domain of regularity in the species spread is transient. At much larger times, after the spreading wedge as well as the other spreading boundaries eventually leave the habitat Ω, the dynamics of the patches becomes chaotic in both space and time (we do not show this pattern for the sake of brevity).

We explored the robustness of the pattern shown in figure 1 with respect to the choice of initial conditions. We found that similar patterns can be observed for the various initial conditions mentioned in §2 and, in particular, when the initial angles of the sectors of species distributions are different and when the species are initially located in six sectors instead of three (results not shown). In other cases, such as when the three species are initially contained within squares or discs, we found that spreading wedges generating regular droplet-like structures are once again produced. The key ingredient for droplet-like structures to occur seems to be that the three species should all meet at at least one point in the domain. However, with highly symmetric initial arrangements, these spreading regions of regularity can eventually annihilate one another by colliding. Figure 2 shows an example of this occurring for the initial conditions with circular subdomains (see figure 2a).

• Download figure
• Open in new tab
• Download PowerPoint

Next, we investigated the dependence of these patterns on the diffusion coefficients ϵ₂ and ϵ₃, using the same initial conditions as for figure 1. The results of this investigation are summarized in figure 3. Decreasing the diffusion coefficient of the second competitor to ϵ₂=0.05 appears to prevent the formation of regular droplets in the wedge of invasion (figure 3a). The spreading sharp wedge of invasion, however, continues to persist in this case. As ϵ₂ is increased, we find that the regular pattern in the wake of the front persists until ϵ₂=0.35; for higher ϵ₂, the spatial structure in the spreading wedge becomes irregular (figure 3b). When the mobility of the third competitor is low (e.g. for ϵ₃=0.35), only a single central droplet structure emerges in the spreading wedge (figure 3c). When ϵ₃ is close to 1, the angle of the spreading wedge abruptly increases and a new pattern of dynamics appears, described below. Finally, for all diffusion coefficients very close to 1, the pattern of spread consists of spiral waves (figure 3d), which are well known for this system.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 4 demonstrates an example of the species spreading with a stripe-like structure, and corresponds to ϵ₂=0.1 and ϵ₃=0.9 with the initial condition shown as for figure 1, but for a larger spatial domain Ω. In this case, we observe that the invasion of the domain of coexistence occurs via the formation of a train of parallel bands which move towards the left-hand side boundary. Each band has the same width and is slightly round at the front. Behind the bands, the species distribution is highly irregular. The spread of the domain of coexistence in the opposite direction (i.e. towards the right boundary) occurs in a different way, via the propagation of irregular patches, which we refer to as a wave of chaos. Thus, there is a clear anisotropy in terms of patterns of the population spread depending on the direction, which in turn is determined by the initial condition. The domain occupied by the irregular patches eventually grows in size until it invades the whole domain once the bands travel out of the domain. One can also see that regular droplet-like patterns travel around the edge of the chaotic domain.

• Download figure
• Open in new tab
• Download PowerPoint

The transition from the pattern of spread via droplet-like units (figure 1) to the one containing bands (figure 4) can be understood by exploring the schematic diagram shown in figure 5. The spread of the droplets in the wedge can be described by considering pairwise interactions of species, most of which actually occur via plane wave interactions. We can neglect the presence of the third species in each case because the density of each species rapidly drops when entering the domain dominated by another species (except the points where all three species meet, as at the tip of the wedge). In figure 5, we show the direction of the spread of plane waves of the cyclic displacement of species; here C_i,j denotes the speed of the plane wave replacing species j by its stronger competitor i. One can also see a round interface between species 1 and species 2. The corresponding wave speed is denoted by V _1,2=V _1,2(R), where R is the radius of the curvature. The values of C_i,j and V _1,2 can be determined by considering the one-dimensional case (in the case of V _1,2, one should explore the system in polar coordinates).

• Download figure
• Open in new tab
• Download PowerPoint

Our simulations show that, for the parameters from figure 1, in the one-dimensional case the propagation of a travelling pulse composed of all three species is impossible, whereas for pairwise switch waves we have C_1,2>C_2,3. The curvature of the wave reduces the spread of the propagation of the front of species 1 in the droplet, thus C_1,2>V _1,2(R)=C_2,3 and the spread of the droplet becomes synchronized with C_2,3 and this gives the condition for R. However, in the case where the diffusion of species 3 increases (as in figure 4), our one-dimensional simulations show that C_2,3>C_1,2, thus C_2,3>V _1,2(R) for any R and the propagation of a travelling pulse composed of all three species now becomes possible. As a result, the droplet-shape structure breaks down and a one-dimensional band composed of three species is eventually formed. In the case of figure 4, our simulations demonstrate that the propagation of one-dimensional pulses is stable with respect to the two-dimensional case, i.e. small two-dimensional perturbations would not destabilize them; however, this does not provide an explanation of why the formation of bands does not occur in the opposite direction (i.e. towards the right), thus more investigation is needed to fully understand this problem.

From figure 1, one can also observe another new type of transient pattern which we call the dynamical droplets. Such structures are formed along the boundary of the domain in which the species coexist. This is essentially a transient regime, however, because, although their life time can be long, the dynamical droplet structures will eventually collide and annihilate one another, resulting in chaotic dynamics. The tip of the dynamical droplet pattern is a point where high densities of all three species meet each other. One can see in the figure that the growing dynamical droplet structure is generated by a certain small area which can be considered as a generator. This is similar to a target wave emanating from a pacemaker; however, the generation of dynamical droplet-shaped moving patches is a transient phenomenon that eventually stops after collision with irregular patches.

We also considered the other scenario of cyclic competition: conditional cyclic competition. This is realized for α_3,1=1.3 (the other interaction parameters being the same as before), so that local interactions between u₃ and u₁ are bistable. Our numerical study revealed a new pattern of invasion, shown in figure 6, for ϵ₂=0.55 and ϵ₃=0.5. The initial spiral tip emanates a droplet composed of species 2 and 3. The droplet takes the form of a glider which moves towards the lower left corner. Some glider-shaped structures quickly disappear, whereas others can persist by splitting and merging with other patches. The tip which generates glider-shaped structures eventually reaches the boundary and disappears. Irregular patches produced in this dynamical regime can persist for a long time; however, eventually they disappear and only one species survives. Thus, this new regime of invasion can only guarantee species coexistence at the moving front of invasion. Our investigation of the parameter space of diffusion coefficients reveals that this pattern is robust and is observed within a 10% variation of ϵ₂=0.55 and ϵ₃=0.5.

• Download figure
• Open in new tab
• Download PowerPoint

Finally, we extend our analysis to the three-dimensional case. We focus on exploring the three-dimensional analogue of the regular droplet-like structures observed in the two-dimensional model under the classical cyclic competition shown in figure 1. Our results when the species are initially distributed within separate cubic regions are shown in figure 7.

• Download figure
• Open in new tab
• Download PowerPoint

We observe also droplet-like regular structures here, although such patterns appear to be prismatic analogues of the two-dimensional structures. We also found that the long-time dynamics is essentially irregular and chaotic: the area occupied by irregular patches of coexisting species will eventually invade the entire domain Ω (the result is not shown here for brevity). Note that similar results are obtained when the species initially dominate within overlapping yet mutually exclusive spherical regions, although for brevity these are not reported here.

5. A posteriori error analysis

The adaptive numerical scheme described in §3, and used to compute the numerical results presented here, relies on a local error indicator functional which is used to decide which mesh elements should be refined or coarsened. This error indicator is derived from a fully computable estimate for the numerical error measured in the L²(0,t;H¹(Ω)) norm in the form of a residual-type a posteriori error estimate, developed in this section.

We develop the a posteriori error analysis for the spatially semidiscrete method (3.3); this choice is commented on at the end of the section. To keep the presentation clear, the analysis we present here focuses on the single-equation situation (i.e. when m=1), although all the results follow for when m>1 in a completely analogous fashion. In this case, ε>0 is simply a constant scalar.

Let f_h:H¹(Ω)→V _h denote the L²(Ω)-orthogonal projection operator given, for any v∈H¹(Ω), by (f_h(v),χ)=(f(v),χ), for all χ∈V _h. Next, we define the discrete Laplacian operator Δ_h:V _h→V _h to be such that, for any v_h∈V _h,

Following [23], we also define the elliptic reconstruction operator $\mathcal{R}:V_h\to H^1(\mathit{\Omega })$ to satisfy, for any v_h∈V _h, where g_h(v_h):=−εΔ_hv_h−f_h(v_h)+f(v_h). As u_h thus coincides with the finite element solution to the elliptic problem with true solution $\mathcal{R}{u}_h$, we assume that we have at our disposal a computable a posteriori error bound of the form for r=0,1, where H⁰=L² (see [21,22] for several examples of such results). For instance, a residual-type a posteriori bound on the error would be of the form where ${\mathcal{T}}_h$ denotes the set of mesh elements, ${\mathcal{E}}_h$ denotes the set of mesh faces, and h_T and h_e denote the diameter of a mesh element and face, respectively. The notation $[[\mathrm{\nabla }{u}_h]]$ is used to denote the jump of ∇u_h over a mesh face, with the definition that, for a point x∈e, $[[\mathrm{\nabla }{u}_h]](\mathit{x})=\mathit{n}\cdot ({(\mathrm{\nabla }{u}_h)}^{+}(\mathit{x})-{(\mathrm{\nabla }{u}_h)}^{-}(\mathit{x}))$, where (∇u_h)⁺(x) and (∇u_h)⁻(x) are the values at x of ∇u_h on each of the elements T⁺ and T⁻ meeting at the face e, respectively, and n denotes the unit vector normal to e pointing from T⁺ to T⁻. In the interests of brevity, throughout the remainder of this section, we shall use $\parallel \cdot \parallel$ to denote the norm on L²(Ω). We are now in a position to derive an a posteriori bound on the error between the exact and the finite element solution.

Theorem 5.1

Let u be the solution of (3.2) and u_h be the solution of (3.3), respectively, under the assumption that the nonlinear term f satisfies the Lipschitz-style growth condition (3.1). The error e:=u−u_h satisfies the $L^{\mathrm{\infty }}(0,t;L^2(\mathit{\Omega }))$ error estimate

and the L²(0,t;H¹(Ω)) error estimate where $M:=\underset{t\in [0,T]}{sup}(1+2{\parallel u_h(t)\parallel }_{L^{\mathrm{\infty }}(\mathit{\Omega })}^{\gamma })$. In each case, the positive constant C is independent of h, u and u_h.

Proof.

Using the definitions of f_h and the discrete Laplacian operator, we can rewrite (3.3) as

implying that u_h satisfies u_h,t−εΔ_hu_h−f_h(u_h)=0 in strong form and therefore, from the definition of the elliptic reconstruction,

To derive the required bounds, we first decompose the error e=u−u_h into $\rho :=u-\mathcal{R}{u}_h$ and $\theta :=\mathcal{R}{u}_h-u_h$. As (5.1) provides a bound on the reconstruction error θ, we focus principally on deriving a bound for ρ. Subtracting (5.3) from the original weak form (3.2), we find that

and therefore, picking v=ρ∈H¹(Ω),

To treat the nonlinear term in (5.4), we use assumption (3.1) on the growth of f to find

The restriction on γ implies that ${\mid a+b\mid }^{\gamma }\le 2^{max\{1,\gamma \}-1}({\mid a\mid }^{\gamma }+{\mid b\mid }^{\gamma })\le 2({\mid a\mid }^{\gamma }+{\mid b\mid }^{\gamma })$, from which, with a=u−u_h and b=u_h, it follows that The first term on the right-hand side of (5.5) can be bounded as In view of bounding the second term on the right-hand side of (5.5), we first note that, as a consequence of Hölder’s inequality, we have ${\int }_{\mathit{\Omega }}{\alpha }^{\gamma +1}\beta \hspace{0.17em}\mathrm{d}\mathit{x}\le ((\gamma +1)/(\gamma +2)){\int }_{\mathit{\Omega }}{\alpha }^{\gamma +2}\hspace{0.17em}\mathrm{d}\mathit{x}+(1/(\gamma +2)){\int }_{\mathit{\Omega }}{\beta }^{\gamma +2}\hspace{0.17em}\mathrm{d}\mathit{x},$ for α,β>0. Moreover, we shall use the result of lemma 5.1 in [24]: that, for any v∈H¹(Ω), we have ${\parallel v\parallel }_{L^{\gamma +2}(\mathit{\Omega })}^{\gamma +2}\le C_{\mathit{\Omega }}{\parallel v\parallel }^{\gamma }{\parallel \mathrm{\nabla }v\parallel }^2$, where C_Ω is a constant depending only on the domain Ω. We have where C₁:=C_fC_Ω2^γ+1((γ+1)/(γ+2)) and C₂:=C₁+C_fC_Ω(γ+2). Applying this and (5.6), (5.4) becomes Integrating through time, we find that where we observe that, by construction, ρ(0)=0, and the functional δ is given by The inequality a^γb≤(a²+b)^1+γ/2 (a consequence of Young’s inequality) then implies that

To bound the final terms on the right-hand side using δ, suppose that the maximum size $h_{max}$ of the mesh used to partition the domain Ω is small enough that, for $h<h_{max}$, the reconstruction error θ satisfies

implying that Consider the set Upon observing that, by construction, we have ρ(0)=0, this set is clearly not empty because it contains τ=0. Moreover, the continuity of ρ in time implies that I must be closed, and thus the maximum of the set is well defined. Thus denoting ${\tau }^{\ast }=maxI$, we suppose that τ*<T. Then, for t≤τ*, Substituting this into (5.7), we find that, for t≤τ*, from which Gronwall’s inequality implies As this is true for all t (and the integral on the left is non-decreasing), it follows that which contradicts the assumption that τ*<T due to the continuity of ρ in time.

Consequently, we have

for any s∈[0,T]. Invoking the triangle inequality and the bound (5.1) for θ then produces the result. ▪

We note that the above error estimates are computable. Indeed, assuming the existence of the constant M bounding the $L^{\mathrm{\infty }}$-norm of the finite element solution u_h is not unreasonable as we can assume to have computed it, and thus have access to its maximal and minimal values. Hence both bounds of theorem 5.1 are computable because they depend only on the discrete solution and problem data.

The error indicator which we use to mark mesh elements for either refinement or coarsening is derived from the a posteriori error bound of theorem 5.1 and is of the form

which may naturally be broken into contributions from each element by observing the structure of E_r in (5.2). We remark that a posteriori bounds for the time discretization by the Crank–Nicolson method are also available [44,45]. Given the nature of the simulations, however, whereby the length scales present do not change over time (but only in position), the incorporation of a full-space time a posteriori analysis in the spirit of [45] was not deemed necessary in this case. Crucially, however, the modified Crank–Nicolson method of [45] was used in the present context of temporal mesh modification.

Figure 8 shows some examples of the computational meshes used to obtain the results of §4, reporting the number of elements saved compared with an equivalent uniform mesh in each case (which may be used as a rough estimate of the computational effort required to compute the solution). What this clearly demonstrates is the effectiveness of the resulting adaptive scheme, because the number of elements required is reduced by over 50% in each case and typically dramatically more. Moreover, examining the areas in which the algorithm has opted to refine or coarsen the mesh indicates that computational effort (in the form of high mesh resolution) is tightly targeted at the areas where it is required, demonstrating that this indicator works well in driving the adaptive algorithm.

• Download figure
• Open in new tab
• Download PowerPoint

6. Discussion and conclusion

Although pattern formation and travelling waves in reaction–diffusion models are widely researched, it seems that—surprisingly—we still lack a clear understanding of all types of patterns in some well-known systems, such as the Lotka–Volterra three-species competition model [12–17]. In this paper, we study such a model with a cyclic competition interaction structure and reveal several novel patterns of population waves in this system which have not been reported before, but may have important applications in biology, chemistry or game-theoretical models. To achieve our goal, we used a novel, adaptive numerical method driven by an a posteriori error estimate which allowed us to efficiently run simulations in both two and three spatial dimensions.

For a cyclic competition system without space, it is well known that the eventual result is a single exclusive species, the densities of the other species being dropped below an extinction threshold [15,26]. Adding a spatial dependence, however, allows for long-term species coexistence. A major question is therefore to understand the dynamical patterns which allow the domain in which the species coexist to spread into areas occupied by a single species. The richness we observe in the system, in terms of its many varied dynamics regimes, stems from the fact that we allow the three species to have different mobilities. Note that this is a natural case which is observed, for example, in the cyclic system of coral reef invertebrates [30], and is modelled here by allowing different diffusion coefficients for each species.

In previous studies of model (2.1), it was reported that the area of mutual coexistence can spread through space via spiral patterns, interacting patches or different types of travelling wedges [13,15,17]. Here we demonstrate new spatially regular patterns of spread. Interestingly, in their work Contento et al. [17] hypothesized the existence of droplet-like structures in a spreading wedge which is close to that shown in figure 1, although they did not find a practical realization of such a pattern and assumed that it would be unstable [17]. Here we found a stable pattern consisting of droplets in a spreading wedge. It is worth observing, however, that in our case the mechanism by which droplets form is slightly different. For example, in our case the spread is based on pairwise interactions between species and, unlike in the cited work, the corresponding one-dimensional case does not allow the spread of a travelling pulse involving all three species. Moreover, the authors of [17] hypothesized that their pattern would exist under conditional cyclic competition, while the droplet-shaped structures in figure 1 are found in classical (i.e. unconditional) cyclic competition. It therefore remains to be determined whether the patterns predicted by Contento et al. are actually possible in the case of conditional cyclic competition involving local bistability. Interestingly, the patterns demonstrated here emerge via a non-Turing scenario: our (linear) stability analysis of the spatially homogenous equilibrium (which is an unstable focus in the system without diffusion) reveals a monotonic decrease in the real part of eigenvalues from positive to negative values with an increase in the wavenumber.

The pattern of spread shown in figure 4 is of particular interest, not only because of the apparent regularity in the direction of movement and irregularity in the opposite direction. A novel feature of the pattern seems to be the coexistence of two different waves moving towards the left-hand boundary: a wave of regularity composed of almost parallel bands in the middle and two wedge-shaped waves of chaos on each side of the bands. Our simulations show a long-term coexistence of both types of waves. Using this pattern, one can describe a complex spread of species involving regular and irregular population patches.

These newly demonstrated patterns of travelling waves with spatially regular structure can be interpreted following the idea of the dynamical stabilization considered in [6,46]. Indeed, the numerical results we present suggest that the regular travelling structures which develop are stable only in a reference frame moving with the speed of the front of the wave, whereas stationary spatially heterogeneous patterns cannot be realized in this system. Note that, compared with the originally reported dynamical stabilization, our patterns have more complex structure. Finally, the regular wave patterns shown here are not globally stable (figure 4) in the sense that, depending on the initial conditions, both the waves of regularity and the waves of chaos can be simultaneously realized in the same system.

The spatially regular geometric shapes found by this study to exist in the wake of spreading waves may have applications in the life sciences. It is well known that the distribution of vegetation in semiarid or other regions shows regular band-shaped patterns which slowly move over time [47–49]. The common point of view on the origin of these vegetation patterns is the interaction between the soil and plants controlled by the level of moisture via various mechanisms such as Turing pattern formation or periodic travelling waves [50,51]. Here we suggest an alternative scenario for the formation of such bands, due to the interaction of competitive plant species which, for instance, does not require the existence of a steady gradient in the system.

The pattern containing regular droplet-shaped structures shown in figure 1, and in figure 3, can potentially be realized on growing domains [10,11] such as in the pigmentation and relief-like patterns found on mollusc shells, which remain a long-standing question [52,53]. Previously, it was suggested that regular patterns in mollusc shells are the result of inhibitor–activator-type interactions via a Turing mechanism. Here we show that similar patterns can be produced by a cyclic competition type of interaction via a non-Turing scenario. Finally, the transient glider-type patterns shown in figure 6 in the case of conditional cyclic competition provide an example of a new scenario of patchy spread of invasive species. This new pattern can be used to improve our understanding and modelling of biological invasions because empirical observations often report that the spread of a species into the habitat occupied by another species occurs via the propagation of irregular patches [54]. This also supports the recent ecological theory of invasional meltdown, when an invasive species facilitates the invasion of some other invasive species [55]. Note that unlike the original concept of invasional meltdown, suggesting mutual facilitation of invasion of species via mutualistic interactions, in our case we consider the case of antagonistic competitors [56].

Our results have also allowed us to improve our understanding of the role of dimensionality on the pattern formation in the considered type of models. This can be seen by comparing figures 1 and 6 alongside the corresponding one-dimensional simulations (not shown here for the sake of simplicity). In one spatial dimension, a wave of mutual spread of three species is impossible for the given parameters: only pairwise switch waves are observed. With two spatial dimensions, the droplet-shape pattern can emerge even though it is simply the result of the pairwise interaction of plane waves, as shown in figure 5 (the round-shaped interface can be formally considered as a plain wave in polar coordinates). Thus, increasing from one to two spatial dimensions allows for species coexistence through a structure which was previously impossible. Interestingly, adding a third dimension continues to allow the persistence of the droplet-shape structure, although we argue here that the pattern remains primarily two dimensional because the three-dimensional droplets are observed to have a prismatic structure, and can still be described via pairwise interactions of locally plain prismatic waves.

Bearing in mind the large domains and long time scales required for the full solution dynamics to evolve, it is clear that the computational cost of these simulations would be prohibitively expensive using a uniform mesh, an issue which is amplified in three spatial dimensions. Instead, the adaptive numerical method described in §3, based on the novel computable error indicator derived in §5, allows us to obtain accurate simulations using just a fraction of the computational effort of an equivalent non-adaptive scheme, as demonstrated by the adapted computational meshes shown in figure 8. The savings this method provides mean that highly accurate simulations of this model are within reach of researchers without needing access to high-performance computing facilities. Moreover, as the error analysis of §5 is applicable to a much wider class of semilinear reaction–diffusion problems, the adaptive method which we describe can be easily applied by researchers wishing to study other phenomena.

We should point out that our numerical investigation of the model cyclic Lotka–Voltera system is by no means exhaustive. We do not claim that combined with the previous studies of the system [13–17] we have now completed a full classification of possible patterns. Further research will be needed especially to further explore the case of conditional cyclic competition. Another interesting direction is to further explore the influence of the number of spatial dimensions on the species persistence. In other words, it is interesting to verify whether or not adding a third dimension will enhance the coexistence of all species and which possible patterns of mutual coexistence can occur. This is a biologically relevant question which is important for understanding, for example, the coexistence of competing bacterial strains or microalgae in laboratory and natural conditions.

Data accessibility

This article has no additional data.

Authors' contributions

A.C., E.H.G. and A.Yu.M. designed and coordinated the study and helped draft the manuscript. O.J.S. designed the study, ran the simulations and helped draft the manuscript. All the authors gave their final approval for publication.

Competing interests

The authors have no competing interests.

Funding

A.C. acknowledges financial support by the EPSRC (grant no. EP/L022745/1). E.H.G. acknowledges financial support by The Leverhulme Trust (grant no. RPG-2015-306).