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

## Abstract

Understanding how patterns and travelling waves form in chemical and biological reaction–diffusion models is an area which has been widely researched, yet is still experiencing fast development. Surprisingly enough, we still do not have a clear understanding about all possible types of dynamical regimes in classical reaction–diffusion models, such as Lotka–Volterra competition models with spatial dependence. In this study, we demonstrate some new types of wave propagation and pattern formation in a classical three species cyclic competition model with spatial diffusion, which have been so far missed in the literature. These new patterns are characterized by a high regularity in space, but are different from patterns previously known to exist in reaction–diffusion models, and may have important applications in improving our understanding of biological pattern formation and invasion theory. Finding these new patterns is made technically possible by using an automatic adaptive finite element method driven by a novel *a posteriori* error estimate which is proved to provide a reliable bound for the error of the numerical method. We demonstrate how this numerical framework allows us to easily explore the dynamical patterns in both two and three spatial dimensions.

### 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

*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

*ψ*_{i},

*i*=1,2,3.

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

*ϵ*

_{2},

*ϵ*

_{3}≤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 *ϵ*_{2},*ϵ*_{3}≤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*^{0}-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

*γ*<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

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

*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*^{1}(*Ω*) to denote the Hilbertian Sobolev space of functions with square-integrable first derivatives, namely

*L*

^{2}(

*Ω*) denotes the usual Lebesgue space of real-valued functions which are square-integrable over

*Ω*, and (

*L*

^{2}(

*Ω*))

^{l}, $l\in \mathbb{N}$, the space of vector fields with components in

*L*

^{2}(

*Ω*). 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*

^{1}(0,

*T*;(

*H*

^{1}(

*Ω*))

^{m}) such that

**∇**

**denotes the Jacobian matrix of**

*u***, and the**

*u**L*

^{2}inner product between matrix-valued functions

*A*(

**) and**

*x**B*(

**) is defined to be $(A,B):={\int}_{\mathit{\Omega}}A:B\mathrm{d}\mathit{x}$, where**

*x**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*^{1}(0,*T*;(*V* _{h})^{m}) such that

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*^{n+1}∈(*V* _{h})^{m} such that

*u*^{n}is the finite element approximation to the solution

**at time moment**

*u**t*

^{n}. 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

### 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*_{2}=1,*a*_{3}=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 *ϵ*_{2}=0.1 and *ϵ*_{3}=0.6. In this and all subsequent figures, we use the following colour coding: dark blue indicates regions dominated by species *u*_{1} (i.e. *u*_{1}≥*u*_{2},*u*_{3}), light yellow indicates regions dominated by species *u*_{2} and mid-green indicates regions dominated by species *u*_{3}. 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*_{1},*u*_{2},*u*_{3}∈[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.

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 2*a*).

Next, we investigated the dependence of these patterns on the diffusion coefficients *ϵ*_{2} and *ϵ*_{3}, 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 *ϵ*_{2}=0.05 appears to prevent the formation of regular droplets in the wedge of invasion (figure 3*a*). The spreading sharp wedge of invasion, however, continues to persist in this case. As *ϵ*_{2} is increased, we find that the regular pattern in the wake of the front persists until *ϵ*_{2}=0.35; for higher *ϵ*_{2}, the spatial structure in the spreading wedge becomes irregular (figure 3*b*). When the mobility of the third competitor is low (e.g. for *ϵ*_{3}=0.35), only a single central droplet structure emerges in the spreading wedge (figure 3*c*). When *ϵ*_{3} 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 3*d*), which are well known for this system.

Figure 4 demonstrates an example of the species spreading with a stripe-like structure, and corresponds to *ϵ*_{2}=0.1 and *ϵ*_{3}=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.

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).

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*_{3} and *u*_{1} are bistable. Our numerical study revealed a new pattern of invasion, shown in figure 6, for *ϵ*_{2}=0.55 and *ϵ*_{3}=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 *ϵ*_{2}=0.55 and *ϵ*_{3}=0.5.

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.

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*^{2}(0,*t*;*H*^{1}(*Ω*)) 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*^{1}(*Ω*)→*V* _{h} denote the *L*^{2}(*Ω*)-orthogonal projection operator given, for any *v*∈*H*^{1}(*Ω*), 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},

*v*

_{h}∈

*V*

_{h},

*g*

_{h}(

*v*

_{h}):=−

*εΔ*

_{h}

*v*

_{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

*r*=0,1, where

*H*

^{0}=

*L*

^{2}(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

*h*

_{T}and

*h*

_{e}denote the diameter of a mesh element and face, respectively. The notation $[\phantom{\rule{-0.15em}{0ex}}[\mathrm{\nabla}{u}_{h}]\phantom{\rule{-0.15em}{0ex}}]$ is used to denote the

*jump*of ∇

*u*

_{h}over a mesh face, with the definition that, for a point

**∈**

*x**e*, $[\phantom{\rule{-0.15em}{0ex}}[\mathrm{\nabla}{u}_{h}]\phantom{\rule{-0.15em}{0ex}}](\mathit{x})=\mathit{n}\cdot ({(\mathrm{\nabla}{u}_{h})}^{+}(\mathit{x})-{(\mathrm{\nabla}{u}_{h})}^{-}(\mathit{x}))$, where (∇

*u*

_{h})

^{+}(

**) and (∇**

*x**u*

_{h})

^{−}(

**) are the values at**

*x***of ∇**

*x**u*

_{h}on each of the elements

*T*

^{+}and

*T*

^{−}meeting at the face

*e*, respectively, and

**denotes the unit vector normal to**

*n**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*

^{2}(

*Ω*). 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*

^{2}

*(0,t;H*

^{1}

*(Ω)) 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

*u*

_{h}satisfies

*u*

_{h,t}−

*εΔ*

_{h}

*u*

_{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

*v*=

*ρ*∈

*H*

^{1}(

*Ω*),

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

*γ*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

*α*,

*β*>0. Moreover, we shall use the result of lemma 5.1 in [24]: that, for any

*v*∈

*H*

^{1}(

*Ω*), 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

*C*

_{1}:=

*C*

_{f}

*C*

_{Ω}2

^{γ+1}((

*γ*+1)/(

*γ*+2)) and

*C*

_{2}:=

*C*

_{1}+

*C*

_{f}

*C*

_{Ω}(

*γ*+2). Applying this and (5.6), (5.4) becomes

*ρ*(0)=0, and the functional

*δ*is given by

*a*

^{γ}

*b*≤(

*a*

^{2}+

*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

*ρ*(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*≤

*τ**,

*t*≤

*τ**,

*t*(and the integral on the left is non-decreasing), it follows that

*τ**<

*T*due to the continuity of

*ρ*in time.

Consequently, we have

*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

*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.

### 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).

## Acknowledgements

This research used the ALICE High Performance Computing Facility at the University of Leicester. The authors thank Ruslan Davidchack (University of Leicester) for his encouragement and support of this project.

### References

- 1
Fisher RA . 1937 The wave of advance of advantageous genes.**Ann. Eugen.**, 355–369. (doi:10.1111/j.1469-1809.1937.tb02153.x) Crossref, Google Scholar**7** - 2
Kolmogorov A, Petrovsky I, Piskunov N . 1937 Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem.**Bull. Moscow State Univ. Ser. A: Math. Mech**, 1–25. Google Scholar**1** - 3
Turing AM . 1952 The chemical basis of morphogenesis.**Phil. Trans. R. Soc. Lond. B**, 37–72. (doi:10.1098/rstb.1952.0012) Link, Web of Science, Google Scholar**237** - 4
Cooper SB, Maini PK . 2012 The mathematics of nature at the Alan Turing centenary.**Interface Focus**, 393–396. (doi:10.1098/rsfs.2012.0018) Link, Web of Science, Google Scholar**2** - 5
Satnoianu RA, Menzinger M, Maini PK . 2000 Turing instabilities in general systems.**J. Math. Biol.**, 493–512. (doi:10.1007/s002850000056) Crossref, PubMed, Web of Science, Google Scholar**41** - 6
Volpert V, Petrovskii S . 2009 Reaction–diffusion waves in biology.**Phys. Life Rev.**, 267–310. (doi:10.1016/j.plrev.2009.10.002) Crossref, PubMed, Web of Science, Google Scholar**6** - 7
Sarfaraz W, Madzvamuse A . 2017 Classification of parameter spaces for a reaction-diffusion model on stationary domains.**Chaos Solitons Fractals**, 33–51. (doi:10.1016/j.chaos.2017.05.032) Crossref, Web of Science, Google Scholar**103** - 8
Tsyganov M, Biktashev V . 2014 Classification of wave regimes in excitable systems with linear cross diffusion.**Phys. Rev. E**, 062912. (doi:10.1103/PhysRevE.90.062912) Crossref, Web of Science, Google Scholar**90** - 9
Aydogmus O . 2015 Patterns and transitions to instability in an intraspecific competition model with nonlocal diffusion and interaction.**Math. Model. Nat. Phenom.**, 17–29. (doi:10.1051/mmnp/201510603) Crossref, Web of Science, Google Scholar**10** - 10
Crampin EJ, Gaffney EA, Maini PK . 2002 Mode-doubling and tripling in reaction-diffusion patterns on growing domains: a piecewise linear model.**J. Math. Biol.**, 107–128. (doi:10.1007/s002850100112) Crossref, PubMed, Web of Science, Google Scholar**44** - 11
Page KM, Maini PK, Monk NA . 2005 Complex pattern formation in reaction–diffusion systems with spatially varying parameters.**Physica D**, 95–115. (doi:10.1016/j.physd.2005.01.022) Crossref, Web of Science, Google Scholar**202** - 12
Petrovskii S, Kawasaki K, Takasu F, Shigesada N . 2001 Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species.**Jpn. J. Ind. Appl. Math.**, 459–481. (doi:10.1007/BF03168586) Crossref, Web of Science, Google Scholar**18** - 13
Mimura M, Tohma M . 2015 Dynamic coexistence in a three-species competition–diffusion system.**Ecol. Complexity**, 215–232. (doi:10.1016/j.ecocom.2014.05.004) Crossref, Web of Science, Google Scholar**21** - 14
Murakawa H, Ninomiya H . 2011 Fast reaction limit of a three-component reaction-diffusion system.**J. Math. Anal. Appl.**, 150–170. (doi:10.1016/j.jmaa.2010.12.040) Crossref, Web of Science, Google Scholar**379** - 15
Adamson MW, Morozov AY . 2012 Revising the role of species mobility in maintaining biodiversity in communities with cyclic competition.**Bull. Math. Biol.**, 2004–2031. (doi:10.1007/s11538-012-9743-z) Crossref, PubMed, Web of Science, Google Scholar**74** - 16
Chen C-C, Hung L-C, Mimura M, Ueyama D . 2012 Exact travelling wave solutions of three-species competition-diffusion systems.**Discrete Contin. Dyn. Syst. Ser. B**, 2653–2669. (doi:10.3934/dcdsb.2012.17.2653) Crossref, Web of Science, Google Scholar**17** - 17
Contento L, Mimura M, Tohma M . 2015 Two-dimensional traveling waves arising from planar front interaction in a three-species competition-diffusion system.**Jpn. J. Ind. Appl. Math.**, 707–747. (doi:10.1007/s13160-015-0186-4) Crossref, Web of Science, Google Scholar**32** - 18
Petrovskii SV, Morozov AY, Venturino E . 2002 Allee effect makes possible patchy invasion in a predator–prey system.**Ecol. Lett.**, 345–352. (doi:10.1046/j.1461-0248.2002.00324.x) Crossref, Web of Science, Google Scholar**5** - 19
Wilson W, McCauley E, De Roos A . 1995 Effect of dimensionality on Lotka-Volterra predator-prey dynamics: individual based simulation results.**Bull. Math. Biol.**, 507–526. (doi:10.1007/BF02460780) Crossref, Web of Science, Google Scholar**57** - 20
Morozov A, Li B-L . 2007 On the importance of dimensionality of space in models of space-mediated population persistence.**Theor. Popul. Biol.**, 278–289. (doi:10.1016/j.tpb.2006.12.005) Crossref, PubMed, Web of Science, Google Scholar**71** - 21
Ainsworth M, Oden JT . 2000**A posteriori error estimation in finite element analysis**.Pure and Applied Mathematics (New York) . New York, NY: Wiley-Interscience. Google Scholar - 22
Verfürth R . 2013*A posteriori error estimation techniques for finite element methods*. Numerical Mathematics and Scientific Computation. Oxford, UK: Oxford University Press. Google Scholar - 23
Makridakis C, Nochetto RH . 2003 Elliptic reconstruction and a posteriori error estimates for parabolic problems.**SIAM J. Numer. Anal.**, 1585–1594. (doi:10.1137/S0036142902406314) Crossref, Web of Science, Google Scholar**41** - 24
Cangiani A, Georgoulis EH, Jensen M . 2013 Discontinuous Galerkin methods for mass transfer through semipermeable membranes.**SIAM J. Numer. Anal.**, 2911–2934. (doi:10.1137/120890429) Crossref, Web of Science, Google Scholar**51** - 25
Cangiani A, Georgoulis EH, Kyza I, Metcalfe S . 2016 Adaptivity and blow-up detection for nonlinear evolution problems.**SIAM J. Sci. Comput.**, A3833–A3856. (doi:10.1137/16M106073X) Crossref, Web of Science, Google Scholar**38** - 26
May RM, Leonard WJ . 1975 Nonlinear aspects of competition between three species.**SIAM J. Appl. Math.**, 243–253. (doi:10.1137/0129022) Crossref, Web of Science, Google Scholar**29** - 27
Kishimoto K . 1982/83 The diffusive Lotka-Volterra system with three species can have a stable nonconstant equilibrium solution.**J. Math. Biol.**, 103–112. (doi:10.1007/BF00275163) Crossref, Web of Science, Google Scholar**16** - 28
Mimura M, Fife PC . 1986 A 3-component system of competition and diffusion.**Hiroshima Math. J.**, 189–207. Crossref, Google Scholar**16** - 29
Sinervo B, Lively CM . 1996 The rock-paper-scissors game and the evolution of alternative male strategies.**Nature**, 240. (doi:10.1038/380240a0) Crossref, Web of Science, Google Scholar**380** - 30
Buss LW, Jackson JBC . 1979 Competitive networks: nontransitive competitive relationships in cryptic coral reef environments.**Am. Nat.**, 223–234. (doi:10.1086/283381) Crossref, Web of Science, Google Scholar**113** - 31
Paquin CE, Adams J . 1983 Relative fitness can decrease in evolving asexual populations of*S. cerevisiae*.**Nature**, 368. Crossref, PubMed, Web of Science, Google Scholar**306** - 32
Kirkup BC, Riley MA . 2004 Antibiotic-mediated antagonism leads to a bacterial game of rock-paper-scissor in vivo.**Nature**, 412. (doi:10.1038/nature02429) Crossref, PubMed, Web of Science, Google Scholar**428** - 33
Hauert C, De Monte S, Hofbauer J, Sigmund K . 2002 Replicator dynamics for optional public good games.**J. Theor. Biol.**, 187–194. (doi:10.1006/jtbi.2002.3067) Crossref, PubMed, Web of Science, Google Scholar**218** - 34
Tainaka K-I . 1993 Paradoxical effect in a three-candidate voter model.**Phys. Lett. A**, 303–306. (doi:10.1016/0375-9601(93)90923-N) Crossref, Web of Science, Google Scholar**176** - 35
Hosono Y . 1998 The minimal speed of traveling fronts for a diffusive Lotka-Volterra competition model.**Bull. Math. Biol.**, 435–448. Crossref, Web of Science, Google Scholar**60** - 36
Alzahrani EO, Davidson FA, Dodds N . 2010 Travelling waves in near-degenerate bistable competition models.**Math. Model. Nat. Phenom.**, 13–35. (doi:10.1051/mmnp/20105502) Crossref, Web of Science, Google Scholar**5** - 37
Cantrell RS, Cosner C . 2004**Spatial ecology via reaction-diffusion equations**. New York, NY: John Wiley & Sons. Crossref, Google Scholar - 38
Adams RA, Fournier JJF . 2003**Sobolev spaces**, 2nd edn.Pure and Applied Mathematics, vol. 140 . Amsterdam, The Netherlands/New York, NY: Elsevier/Academic Press. Google Scholar - 39
Wloka J . 1987**Partial differential equations**. Cambridge, UK: Cambridge University Press. [Translated from the German by C. B. Thomas and M. J. Thomas.] Crossref, Google Scholar - 40
Ciarlet PG . 2002**The finite element method for elliptic problems**.Classics in Applied Mathematics, vol. 40 . Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)]. Google Scholar - 41
Hundsdorfer W, Verwer J . 2003**Numerical solution of time-dependent advection-diffusion-reaction equations**.Springer Series in Computational Mathematics, vol. 33 . Berlin, Germany: Springer-Verlag. Google Scholar - 42
Madzvamuse A . 2006 Time-stepping schemes for moving grid finite elements applied to reaction-diffusion systems on fixed and growing domains.**J. Comput. Phys.**, 239–263. (doi:10.1016/j.jcp.2005.09.012) Crossref, Web of Science, Google Scholar**214** - 43
Bangerth W, Hartmann R, Kanschat G . 2007 deal.II—a general-purpose object-oriented finite element library.**ACM Trans. Math. Softw.**, Art. 24, p. 27. (doi:10.1145/1268776.1268779) Crossref, Web of Science, Google Scholar**33** - 44
Akrivis G, Makridakis C, Nochetto RH . 2006 A posteriori error estimates for the Crank-Nicolson method for parabolic equations.**Math. Comp.**, 511–531. Crossref, Web of Science, Google Scholar**75** - 45
Bänsch E, Karakatsani F, Makridakis C . 2012 A posteriori error control for fully discrete Crank-Nicolson schemes.**SIAM J. Numer. Anal.**, 2845–2872. (doi:10.1137/110839424) Crossref, Web of Science, Google Scholar**50** - 46
Dagbovie AS, Sherratt JA . 2014 Absolute stability and dynamical stabilisation in predator-prey systems.**J. Math. Biol.**, 1403–1421. (doi:10.1007/s00285-013-0672-8) Crossref, PubMed, Web of Science, Google Scholar**68** - 47
Montana C . 1992 The colonization of bare areas in two-phase mosaics of an arid ecosystem.**J. Ecol.**, 315–327. (doi:10.2307/2261014) Crossref, Web of Science, Google Scholar**80** - 48
Lefever R, Lejeune O . 1997 On the origin of tiger bush.**Bull. Math. Biol.**, 263–294. (doi:10.1007/BF02462004) Crossref, Web of Science, Google Scholar**59** - 49
Rietkerk M, Van de Koppel J . 2008 Regular pattern formation in real ecosystems.**TREE**, 169–175. (doi:10.1016/j.tree.2007.10.013) PubMed, Web of Science, Google Scholar**23** - 50
Klausmeier CA . 1999 Regular and irregular patterns in semiarid vegetation.**Science**, 1826–1828. (doi:10.1126/science.284.5421.1826) Crossref, PubMed, Web of Science, Google Scholar**284** - 51
Dagbovie AS, Sherratt JA . 2014 Pattern selection and hysteresis in the Rietkerk model for banded vegetation in semi-arid environments.**J. R. Soc. Interface**, 20140465. (doi:10.1098/rsif.2014.0465) Link, Web of Science, Google Scholar**11** - 52
Meinhardt H . 1992 Pattern formation in biology: a comparison of models and experiments.**Rep. Prog. Phys**, 797–849. (doi:10.1088/0034-4885/55/6/003) Crossref, Web of Science, Google Scholar**55** - 53
Fowler DR, Meinhardt H, Prusinkiewicz P . 1992 Modeling seashells.**SIGGRAPH Comput. Graph.**, 379–387. (doi:10.1145/142920.134096) Crossref, Google Scholar**26** - 54
Shigesada N, Kawasaki K . 1997**Biological invasions: theory and practice**. Oxford, UK: Oxford University Press. Crossref, Google Scholar - 55
O’Dowd DJ, Green PT, Lake PS . 2003 Invasional ‘meltdown’ on an oceanic island.**Ecol. Lett.**, 812–817. (doi:10.1046/j.1461-0248.2003.00512.x) Crossref, Web of Science, Google Scholar**6** - 56
Simberloff D, Von Holle B . 1999 Positive interactions of nonindigenous species: invasional meltdown?**Biol. Invasions**, 21–32. (doi:10.1023/A:1010086329619) Crossref, Google Scholar**1**