Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch article

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

A. Cangiani

A. Cangiani

Department of Mathematics, University of Leicester, University Road, Leicester LE1 7RH, UK

Google Scholar

Find this author on PubMed

,
E. H. Georgoulis

E. H. Georgoulis

Department of Mathematics, University of Leicester, University Road, Leicester LE1 7RH, UK

Department of Mathematics, School of Mathematical and Physical Sciences, National Technical University of Athens, Zografou 15780, Greece

Google Scholar

Find this author on PubMed

,
A. Yu. Morozov

A. Yu. Morozov

Department of Mathematics, University of Leicester, University Road, Leicester LE1 7RH, UK

[email protected]

Google Scholar

Find this author on PubMed

and
O. J. Sutton

O. J. Sutton

Department of Mathematics and Statistics, University of Reading, Whiteknights, PO Box 220, Reading RG6 6AX, UK

Google Scholar

Find this author on PubMed

    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 [47]. 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 [1217]. 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 [2123]. 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 ΩRd, where d=2 or 3, reads

    (u1)t=D1Δu1+a1u1(1α1,1u1α1,2u2α1,3u3),(u2)t=D2Δu2+a2u2(1α2,1u1α2,2u2α2,3u3)and(u3)t=D3Δu3+a3u3(1α3,1u1α3,2u2α3,3u3),
    where uiui(t,x) is the density of species i at time t and location x; ai is the intrinsic growth rate of species ui, and each coefficient αi,j represents the limiting effect that the presence of species uj has on species ui (the term αi,i describes the self-limitation of the population). The diffusion coefficients Di 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

    ui(x,0)=ψi(x)in Ω,
    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 u~i=αi,iui, and then rescale both time and space to obtain the simplified three-species competition system,

    (u1)t=Δu1+u1(1u1α1,2u2α1,3u3),(u2)t=ε2Δu2+a2u2(1α2,1u1u2α2,3u3)and(u3)t=ε3Δu3+a3u3(1α3,1u1α3,2u2u3).}2.1
    Note that we can assume without loss of generality that ϵ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 [1217,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 C0-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 mN 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

    utΔ(εu)=f(u)in Ω,uin=0on Ωi=1,2,,mandu(x,0)=ψ(x)in Ω.
    The nonlinear interaction term f:RmRm is assumed to satisfy the general Lipschitz-style growth condition that there exists 0≤γ<2 when d=2, or 0γ43 when d=3, such that, for any v,wRm, we have
    f(v)f(w)Cf(1+vγ+wγ)vw,3.1
    with denoting the Euclidean norm on Rm.

    Let the matrix ARm×m contain the interaction parameters, with Ai,j=αi,j. The interaction term of the model (2.1) may be seen to satisfy this condition with γ=1 by writing it as

    f(u)=(a^u)(1Au),
    where 1Rm denotes the vector with all entries equal to 1, ° denotes the Hadamard product between tensors and for rRm we use 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 H1(Ω) to denote the Hilbertian Sobolev space of functions with square-integrable first derivatives, namely

    H1(Ω)={vL2(Ω):v(L2(Ω))d}.
    Here, L2(Ω) denotes the usual Lebesgue space of real-valued functions which are square-integrable over Ω, and (L2(Ω))l, lN, the space of vector fields with components in L2(Ω). 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 uH1(0,T;(H1(Ω))m) such that

    (ut,v)+(εu,v)(f(u),v)=0v(H1(Ω))m.3.2
    Here, u denotes the Jacobian matrix of u, and the L2 inner product between matrix-valued functions A(x) and B(x) is defined to be (A,B):=ΩA:Bdx, 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 Th 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 uhH1(0,T;(V h)m) such that

    (uh,t,χ)+(εuh,χ)(f(uh),χ)=0χ(Vh)m.3.3
    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≤nN−1, find un+1∈(V h)m such that

    (un+1unτ,χ)+(ε2(un+1+un),χ)=(32f(un)12f(un1),χ)χ(Vh)m,
    where the function un is the finite element approximation to the solution u at time moment tn. 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.

    Inline Graphic

    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 a2=1,a3=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 u1 (i.e. u1u2,u3), light yellow indicates regions dominated by species u2 and mid-green indicates regions dominated by species u3. 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 u1,u2,u3∈[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.

    Figure 1.

    Figure 1. The evolution of the ‘triangular droplet-like’ patterns which are observed when ϵ2=0.1 and ϵ3=0.6 using the ‘triangular’ initial conditions. (a) t=50, (b) t=150, (c) t=360 and (d) t=900.

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

    Figure 2.

    Figure 2. The evolution of the ‘triangular droplet-like’ patterns with ϵ2=0.1 and ϵ3=0.6 using the initial conditions shown in (a). (a) t=0, (b) t=40, (c) t=120, (d) t=240, (e) t=400 and (f) t=600.

    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 3a). 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 3b). 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 3c). 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 3d), which are well known for this system.

    Figure 3.

    Figure 3. The formation of droplet patterns for various values of the diffusion parameters ε2 and ε3. Solutions are plotted for t=900, unless otherwise stated, having evolved from the same initial conditions as shown in figure 1. (a) ϵ2=0.05, ϵ3=0.55, (b) ϵ2=0.35, ϵ3=0.55 (shown at t=560), (c) ϵ2=0.1, ϵ3=0.35 and (d) ϵ2=ϵ3=1.

    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.

    Figure 4.

    Figure 4. The evolution of the ‘band’ patterns which are observed when ε= (1,0.1,0.9) and the same initial conditions as in figure 1. (a) t=40, (b) t=160, (c) t=320, (d) t=640, (e) t=960 and (f) t=1400.

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

    Figure 5.

    Figure 5. Schematic of the movement of a droplet-like unit.

    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 C1,2>C2,3. The curvature of the wave reduces the spread of the propagation of the front of species 1 in the droplet, thus C1,2>V 1,2(R)=C2,3 and the spread of the droplet becomes synchronized with C2,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 C2,3>C1,2, thus C2,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 u3 and u1 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.

    Figure 6.

    Figure 6. The evolution of the ‘glider’ patterns which are observed when ϵ2=0.55 and ϵ3=0.5, α3,1=1.3 and the initial conditions are the same as in figure 1. The three colours indicate the regions of the domain in which each of the three competing species dominates. (a) t=60, (b) t=120, (c) t=160 and (d) t=290.

    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.

    Figure 7.

    Figure 7. Evolution in three spatial dimensions using the parameters ϵ2=0.1 and ϵ3=0.6 (as in figure 1), computed on the domain Ω=[0,600]3. Column (a) shows the domains dominated by each species, while those in (b) show u2 and u3 only in the subdomain x>300.

    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 L2(0,t;H1(Ω)) 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 fh:H1(Ω)→V h denote the L2(Ω)-orthogonal projection operator given, for any vH1(Ω), by (fh(v),χ)=(f(v),χ), for all χV h. Next, we define the discrete Laplacian operator Δh:V hV h to be such that, for any vhV h,

    (εΔhvh,χ)=(εvh,χ)χVh.
    Following [23], we also define the elliptic reconstruction operator R:VhH1(Ω) to satisfy, for any vhV h,
    (εRvh,v)=(gh(vh),v)vH1(Ω),
    where gh(vh):=−εΔhvhfh(vh)+f(vh). As uh thus coincides with the finite element solution to the elliptic problem with true solution Ruh, we assume that we have at our disposal a computable a posteriori error bound of the form
    RuhuhHrCelipEr(uh),5.1
    for r=0,1, where H0=L2 (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
    Er(uh)2=TThhT42rεΔuhεΔhuhfh(uh)+f(uh)L2(T)2+eEhhe32rε[[uh]]L2(e)2,5.2
    where Th denotes the set of mesh elements, Eh denotes the set of mesh faces, and hT and he denote the diameter of a mesh element and face, respectively. The notation [[uh]] is used to denote the jump of ∇uh over a mesh face, with the definition that, for a point xe, [[uh]](x)=n((uh)+(x)(uh)(x)), where (∇uh)+(x) and (∇uh)(x) are the values at x of ∇uh 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 to denote the norm on L2(Ω). 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 uh 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−uh satisfies the L(0,t;L2(Ω)) error estimate

    Cmaxs[0,t]e(s)2maxs[0,t]E0(uh(s))+e(1+4MCf)t0tE0(uh,t(s))2+4MCfE0(uh(s))2ds+e(1+4MCf)t0tC1E0(uh(s))γE1(uh(s))2ds
    and the L2(0,t;H1(Ω)) error estimate
    C0tεe2ds0tεE1(uh(s))2ds+e(1+4MCf)t0tE0(uh,t(s))2+4MCfE0(uh(s))2ds+e(1+4MCf)t0tC1E0(uh(s))γE1(uh(s))2ds,
    where M:=supt[0,T](1+2uh(t)L(Ω)γ). In each case, the positive constant C is independent of h, u and uh.

    Proof.

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

    (uh,tεΔhuhfh(uh),χ)=0χVh,
    implying that uh satisfies uh,tεΔhuhfh(uh)=0 in strong form and therefore, from the definition of the elliptic reconstruction,
    (uh,t,v)+(εRuh,v)+(f(uh),v)=0vH1(Ω).5.3

    To derive the required bounds, we first decompose the error e=uuh into ρ:=uRuh and θ:=Ruhuh. 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

    (ρt,v)+(ερ,v)=(θt,v)+(f(u)f(uh),v)vH1(Ω),
    and therefore, picking v=ρH1(Ω),
    12ddt(ρ2)+ερ2θt2+ρ2+(f(u)f(uh),ρ).5.4

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

    (f(u)f(uh),ρ)Ωf(u)f(uh)ρdxCfΩ(1+uγ+uhγ)uuhρdx.
    The restriction on γ implies that a+bγ2max{1,γ}1(aγ+bγ)2(aγ+bγ), from which, with a=uuh and b=uh, it follows that
    (f(u)f(uh),ρ)2CfΩMuuhρ+uuhγ+1ρdx.5.5
    The first term on the right-hand side of (5.5) can be bounded as
    ΩuuhρdxΩρ2+θρdx32ρ2+12θ2.5.6
    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 Ωαγ+1βdx((γ+1)/(γ+2))Ωαγ+2dx+(1/(γ+2))Ωβγ+2dx, for α,β>0. Moreover, we shall use the result of lemma 5.1 in [24]: that, for any vH1(Ω), we have vLγ+2(Ω)γ+2CΩvγv2, where CΩ is a constant depending only on the domain Ω. We have
    CfΩuuhγ+1ρdxCfγ+1γ+2Ωuuhγ+2dx+Cfγ+2Ωργ+2dxCf(γ+1)2γ+1γ+2θLγ+2(Ω)γ+2+Cf1+(γ+1)2γ+1γ+2ρLγ+2(Ω)γ+2C1θγθ2+C2ργρ2,
    where C1:=CfCΩ2γ+1((γ+1)/(γ+2)) and C2:=C1+CfCΩ(γ+2). Applying this and (5.6), (5.4) becomes
    12ddt(ρ2)+ερ2θt2+4MCfθ2+(1+4MCf)ρ2+C1θγθ2+C2ργρ2.
    Integrating through time, we find that
    ρ(t)2+0tερ2dtδ(θ)2+(1+4MCf)0tρ2dt+C20tργρ2dt,
    where we observe that, by construction, ρ(0)=0, and the functional δ is given by
    δ(θ)2:=0Tθt2+4MCfθ2+C1θγθ2dt.
    The inequality aγb≤(a2+b)1+γ/2 (a consequence of Young’s inequality) then implies that
    ρ(t)2+0tερ2dsδ(θ)2+(1+4MCf)0tρ2ds+C2sups[0,t]ρ(s)γ0tρ2dtδ(θ)2+(1+4MCf)0tρ2ds+C2(sups[0,t]ρ(s)2+0tρ2ds)1+γ/2.5.7

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

    δ(θ)C2γ(4e(1+4MCf)T)(2+γ)/2γ,
    implying that
    C2(4δ(θ)2e(1+4MCf)T)1+γ/2δ(θ)2.
    Consider the set
    I={τ[0,T]:sups[0,τ]ρ(s)2+0τερ2ds4δ(θ)2e(1+4MCf)T}.
    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 τ=maxI, we suppose that τ*<T. Then, for tτ*,
    C2(sups[0,t]ρ(s)2+0tερ2ds)1+γ2C2(4δ(θ)2e(1+4MCf)T)1+γ2δ(θ)2.
    Substituting this into (5.7), we find that, for tτ*,
    ρ(t)2+0tερ2ds2δ(θ)2+(1+4MCf)0tρ2ds,
    from which Gronwall’s inequality implies
    ρ(t)2+0tερ2ds2δ(θ)2e(1+4MCf)T.
    As this is true for all t (and the integral on the left is non-decreasing), it follows that
    sups[0,t]ρ(s)2+0tερ2ds2δ(θ)2e(1+4MCf)T,
    which contradicts the assumption that τ*<T due to the continuity of ρ in time.

    Consequently, we have

    ρ(t)2+0tερ2ds2e(1+4MCf)s0tθt2+4MCfθ2+C1θγθ2ds,
    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-norm of the finite element solution uh 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

    ϵE1(un)2+E0(unun1τ)2,
    which may naturally be broken into contributions from each element by observing the structure of Er 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.

    Figure 8.

    Figure 8. Some examples of meshes produced by the adaptive algorithm, demonstrating the reduction in the number of elements required compared with a uniform mesh with the same resolution around the layers. (a) Adapted mesh for the solution shown in figure 2d containing 18 418 elements (73% saving over an equivalent uniform mesh). (b) Adapted mesh for the solution shown in figure 3d, containing 8935 elements (83% saving over an equivalent uniform mesh). (c) Adapted mesh for the solution shown in figure 6d, containing 5953 elements (91% saving over an equivalent uniform mesh). (d) Adapted mesh for the solution shown in figure 7, containing 1 015 660 elements (51% saving over an equivalent uniform mesh).

    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 [1217]. 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 [4749]. 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 [1317] 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.

    Footnotes

    Published by the Royal Society. All rights reserved.

    References