Stationary localised patterns without Turing instability

Since the pioneering work of Turing, it has been known that diffusion can destablise a homogeneous solution that is stable in the underlying model in the absence of diffusion. The destabilisation of the homogeneous solutions leads to the generation of patterns. In recent years, techniques have been developed to analyse so‐called localised spatial structures. These are solutions in which the spatial structure occurs in a localised region. Unlike Turing patterns, they do not spread out across the whole domain. We investigate the existence of localised structures that occur in two predator‐prey models. The functionalities chosen have been widely used in the literature. The existence of localised spatial structures has not investigated previously. Indeed, it is easy to show that these models cannot exhibit the Turing instability. This has perhaps led earlier researchers to conclude that interesting spatial solutions can therefore not occur for these models. The novelty of our paper is that we show the existence of stationary localised patterns in systems which do not undergo the Turing instability. The mathematical tools used are a combination of Linear and weakly nonlinear analysis with supporting numerical methods. By combining these methods, we are able to identify conditions for a wide range of increasing exotic behaviour. This includes the Belyckov‐Devaney transition, a codimension two spatial instability point and the formation of localised patterns. The combination of spectral computations and numerical simulations reveals the crucial role played by the Hopf bifurcation in mediating the stability of localised spatial solutions. Finally, numerical solutions in two spatial dimensions confirms the onset of intricate spatio‐temporal patterns within the parameter regions identified within one spatial dimension.


INTRODUCTION
The study of predator-prey models has a long history in mathematical ecology, with the first investigations appearing in the late 1920s. 1,2 At their most basic, such models take the following form.
Rate of change of the prey population (u): Rate of change of the predator population (v): Here, the function G 1 (u) models the growth of the prey in the absence of the predator, the function P(u, v) models the increase in mortality of the prey due to predation, the function G 2 (v) models the growth of the predator in the absence of prey and the function G 3 (u, v) models the growth of the predator due to predation upon the prey. (Often, the predator growth function G 3 is proportional to the predation function P). Such models remain a motivator of much research since there are many twists and turns in determining the functionality of the terms on the right-hand side that are appropriate for a given predator and prey. For a historic overview of research in this area, we refer to Berryman. 3 A natural extension of the basic predator-prey framework (1) is to allow both species to move at random throughout an otherwise homogeneous environment. With zero-flux boundary conditions, steady-state solutions of system (1) are also steady-state solutions of the equivalent reaction-diffusion model. It has been known since the pioneering work of Turing 4 that diffusion can destabilise a (spatially uniform) solution that is stable in model (1). This instability, known as the Turing instability, leads to the formation of patterns. There has been much interest in using the Turing instability as a mechanism to generate the phenomenon of self-organisation in many areas of mathematical biology. 5 There is a vast literature on pattern formation through the Turing instability in predator-prey systems. It is not possible to cover this extensive field of research here. A good overview for the beginner is provided by the review article of Banerjee. 6 We now draw the reader's attention to some significant papers suggested by anonymous reviewers of this paper. Djilali and Bentout 7 investigate pattern formation in a predator-prey model involving prey social behaviour. The homogeneous coexistence steady-state is shown to lose stability at a Hopf bifurcation. The existence of a non-homogeneous steady-state solution is established rigorously using Leray-Schauder degree theory.
Djilali 8 investigates a predator-prey model with a prey herd shape effect and cross diffusion. The model exhibits rich range of dynamic behaviour including a Hopf bifurcation, the Turing instability, and the Turing-Hopf instability. The spatiotemporal behaviour in the vicinity of this last point is investigated using normal form theory. An interesting feature of this model is that the Turing instability is impossible in the absence of cross-diffusion. Another paper which investigates the effect of cross-diffusion upon pattern formation in predator-prey systems is due to Sun et al. 9 Souna et al 10 investigate a predator-prey model in which the prey moves in a herd. When the predators attack a fraction of prey, they try to escape by splitting off from the herd. The predator-prey functionality is consequently a convex combination of two functions: one for predation upon the herd and one for predation upon the escapees. Although this system exhibits a Hopf bifurcation, the authors show that the Turing instability is impossible.
Souna et al 11 investigate a predator-prey model in which the prey population exhibits social behaviour whilst the predator population is subject to a quadratic harvesting term. The authors show that the non-trivial equilbrium can undergo a Hopf bifurcation, a Turing bifurcation, and a Turing-Hopf bifurcation. The authors investigate the dynamic behaviour near the Hopf bifurcation by determining its normal form.
However, the Turing instability is not the only mechanism that can lead to spatial patterns. In recent years, there has been a growing interest within the dynamical systems community in the formation of localised patters in reaction-diffusion equations. 12,13 These localised patterns are a result of a bistability between the homogeneous steady state and the spatially periodic patterned state. [13][14][15] Such patterns have been found in models from a variety of fields including biology, 13,14 ecology, 16,17 fluid dynamics, 18 and optics. 19 In this paper, we consider two predator-prey systems which can not lose stability through the Turing mechanism. 20,21 We show that they do lose stability at a more general type of spatial instability. This leads to the formation of spatially localised solutions. We investigate numerically the existence of such patterns in the parameter space. Our results reinforce the growing recognition that the Turning instability is not a sine qua non for the formation of spatial patterns.
The specific models that we investigate are the scaled Ivlev model.
and the scaled Holling II model These models are investigated in either one or two dimensions with homogeneous Neumann boundary conditions. These systems have been chosen since neither of them does the homogeneous steady-state solution lose stability at a Turing bifurcation. The rest of the paper is organised as follows. In Section 2, we carry out preliminary analysis using a generalised predator-prey model that employs a system of reaction-diffusion equations in one or two spatial dimensions. Specifically, through the employment of linear stability analysis, the necessary conditions under which a homogeneous solution losses stability at a Hopf bifurcation are determined. This spatial instability is not the classical Turing instability. Rather, it is a more general spatial instability, 22 although like the Turing instability, it is generically a codimension one bifurcation. Through the use of weakly nonlinear analysis, we obtain an analytical expression for a curve in parameter space along which the codimension one spatial instability bifurcation exhibits a degeneracy, i.e., a codimension two spatial instability. This is where the spatial instability curve bifurcates from being subcritical to supercritical.
In Sections 3 and 4, we investigate the Ivlev predator-prey model and the Holling II predator-prey model, respectively. In both sections, we establish the existence of a region in which localised solutions are found. By using specific parameter values, the bifurcation diagram in a two parameter space is produced. Through simulations and numerical continuation, stable localised patterns are generated for one spatial dimension (both models) and two spatial dimension (the Ivlevl model) for both one and two spatial dimensions. Concluding remarks are given in Section 5. We emphasise that the novelty of this work is to investigate localised pattern formation in two well-studied systems that do not undergo a Turing instability.
Finally, we note that the techniques we use have been developed for models using integer valued derivatives. In recent years, there has been growing interest in the development of numerical methods for problems employing fractional time derivatives, such as Kumar and Zeidan. 23 The exploration of localised solutions within such a modelling framework may be fruitful.

MATHEMATICAL PRELIMINARIES
In this section, we set up the mathematical tools that are used in Sections 3 and 4 to investigate localised structures in two predator-prey models.
We start in Section 2.1 by discussing the formulation of the predator-prey models. We do this in a general setting, i.e., without specifying the precise formulation of the reaction functions. Following this, we discuss homogeneous steady-state solutions of the in Section 2.2. In Section 2.3, we outline the linear stability analysis of a homogeneous steady-state solution. This includes a discussion of two important bifurcations and the Belyckov-Devaney curve.
The linear analysis is extended to a weakly nonlinear analysis in Section 2.4. Modern symbolic calculation packages allow the ready application of normal form theory to problems of the type considered here. Finally, in Section 2.5, we outline how reformulating the steady-state problems via the spatial dynamics approach allows the use of the continuation software AUTO to investigate non-homogeneous solutions.

Generalised model equations
Both of our models consist of two coupled nonlinear reaction-diffusion equations, in either one or two spatial dimensions. In the model, equations u and v are population densities, t is time, and the spatial co-ordinates are x and possibly .
for u(x, , t) and v(x, , t) ∈ R and x, ∈ (−∞, ∞) with asymptotically flat boundary conditions In this system, the Laplacian operator in 2D space is given by whilst (u, v) and g(u, v) are generic functions that are specified in later sections. The parameter is the ratio of the diffusion constant for species u divided by that of species v. It is assumed that < 1. Consequently, the species v is more mobile than that of species u.
For numerical calculations, we simulate the model on a finite domain x, ∈ [−L, L] with L ≫ 1 and replace the asymptotic boundary conditions with the Neumann boundary conditions

Homogeneous steady states
A natural starting place to analyse model (2) is to investigate homogeneous steady state solutions by setting both spatial and time derivatives to zero, i.e., by solving the system of equations: If it exists, a non-zero homogeneous coexistence steady state solution is written: This solution is only of interest when the values of u * and v * are positive. This requirement will impose conditions on some of the model parameters.

Linear stability analysis
In Section 2.3.1, we review the standard linear stability analysis for a homogeneous steady-state solution. In Section 2.3.2, we discuss bifurcations of the homogeneous steady-state solution.

Linear stability analysis of the homogeneous steady-state solution
In this section, we outline the linear stability analysis of a homogeneous steady-state solution, following Cross and Hohenberg. 22 Doing this, it is assumed that the solution for model (2) is near to a homogeneous steady-state solution (5).
In this equation, represents the linear operator of (2) with J being the Jacobian matrix evaluated at (5), D is a diagonal matrix that accounts for the diffusion coefficients associated with the second order spatial derivatives, and NL(U, V)| S + is a scalar quantity that reflects the nonlinear terms in (2) after using the transformation (6).
The system (7) is now linearised by neglecting the nonlinear terms. The stability of the homogeneous solution is investigated using the ansatz In the ansatz, the constants a and b, with ||(a, b)|| ≪ 1, are regarded as the real-valued amplitudes of the periodic solution for (U, V). The parameters k and are the wave number and the growth factor or temporal eigenvalue, respectively, and i is the usual imaginary unit.
Substituting the ansatz (8) into the linear form of (7) gives a dispersion equation that can be written in the form where T k and D k are the trace and determinant of the matrix . The expression for in indicates that the spectrum of eigenvalues is dependent upon the wave number k. The homogeneous steady-state solution is stable provided that the real part of all the eigenvalues is negative; it is unstable if there is at least one eigenvalue with positive real part.

Bifurcations of the homogeneous steady-state solution
In this section, we discuss two kinds of bifurcations that lead to symmetry-breaking solutions: the homogeneous Hopf bifurcation and the spatial instability (Turing instability) bifurcation. The Hopf bifurcation breaks the temporal symmetry of (2) to give oscillations that are stationary in space but periodic in time. The Hopf bifurcation occurs when Combining the conditions (10) with (9), we can obtain an expression for the Hopf bifurcation when k = 0. Note that S + in (5) is stable to the spatially homogeneous mode k = 0. The spatial instability bifurcation breaks the spatial symmetry of (2) to give solutions that are periodic in space but uniform in time. This bifurcation occurs when the following critical conditions are satisfied 22 : The value of the wave number (k) at which these conditions holds (k = k c ) is known as the critical value of the wave number.
Combining the conditions (11) with (9) gives two expressions. The first expression defines a bifurcation curve for the spatial instability of a homogeneous steady-state solution of system (2). The second expression is for the critical wave number value (k c ).
It is not widely recognised that the conditions 11 provide a more general mechanism for a homogeneous steady-state solution to undergo a spatial instability than the classical Turing bifurcation. A system that cannot undergo a Turing bifurcation may still undergo the generalised spatial instability. 22 Another way of looking at the generalised spatial instability is to consider the steady-state version of system (2). The system of two second-order differential equations in the spatial co-ordinate can be rewritten as a system of four first-order differential equations in space. In the spatial dynamics approach, the spatial co-ordinate is viewed as a 'time-like' variable. In this approach, the spatial instability conditions (11) are explained as a Hamiltonian-Hopf bifurcation-provided that the spatial co-ordinate extends over an infinite domain. 14 From the perspective of spatial dynamics, the pattern solutions are independent of time. The ansatz solution for (7) is now written in the form: where Φ is the vector obtained from the Ker (7) and A is the constant amplitude of the pattern. Another interesting transition can be identified from the linearised analysis of the spatial dynamics problem, namely, a curve along which the spatial eigenvalues change from being purely real to complex conjugates. This transition is called the Belyckov-Devaney (BD) curve. [24][25][26] This curve is not represent a bifurcation. Rather, the solution of the localised patterns changes from having a monotonic decaying tail to an oscillatory one.

Weakly nonlinear analysis
The parameters that govern the system's solution behaviour are contained within the functions and g. We identify two of these as distinguished parameters. Once identified, it is then possible to carry out a weakly nonlinear analysis using normal form theory to identify points that are codimension two degeneracies. Here, we sketch some of the more prominent details and provide relevant explanations.
Our focus is on obtaining analytical expression for a particular codimension two points at which the spatial instability curve changes from being a subcritical to a supercritical bifurcation. This point is identified by undertaking a normal form analysis, as detailed in Haragus and Iooss 25 and applied in Al Saadi et al 27 and Verschueren and Champneys. 28 The first step is to determine the amplitude equation arising from any spatial instability. To do this, we assume in (12) that A = A(t) and seek the amplitude and solution of (7) in the form where R (i) (A) and P (i) (A(t), x) for i = 1, 2, 3, … are functions that represent the decreasing orders of magnitude in A.
The required calculations are standard, but lengthy, and accordingly, we present only the essential steps in this procedure. The system (7) is now solved for each order in A. For example, the first-order coupled equation between R 1 (A) and P 1 (A, x) is given by where c.c. means the complex conjugate of A P (1) The first-order approximation is simply the linear approximation of (7). Consequently, the amplitude is constant, and hence, t R (1) (A) = 0. The first-order solution is therefore given by (12).
To investigate pattern solutions, we must extend this analysis in A to the third order. After the required calculations (see Al Saadi et al 27 for more details), the following equation is obtained: The values of the constants C i , i = 1, 3, are found from solvability conditions which can be represented in term of all the problem parameters. In Equation 15, is an unfolding parameter which allows a small variation around the homogeneous steady state. Once an expression for the parameter C 3 is achieved, we can determine the conditions for which it changes sign: This gives the point at which the spatial instability bifurcation changes nature, from subcritical to supercritical.

Nonhomogeneous solution
In Sections 2.2-2.4, we have focused our discussion on the understanding of the homogeneous solutions. In this section, we discuss the investigation of time-independent solutions of system (2). To simplify our account, we focus on the 1D scenario. For this case, the system (2) reduces to the following coupled system of two nonlinear spatial second-order differential equations (ODEs) in one spatial coordinate x.
To investigate the spatial dynamics of (16), we rewrite it as a system of four first order non-homogeneous nonlinear differential equations of the form In this system, , and B(x, Z) is a 4 × 1 vector-valued function of four variables that is dependent upon the functions and g. An important observation is that the system (2) is spatially invariant under reflection. This leads to the invariance of the dynamical system (17). It follows that system (2) is spatially reversible. 25,26 The advantage of this reformulation is that it allows the use of continuation methods. This allows the parameter values over which particular types of localised structures occur to be determined. Furthermore, they allow a relationship to be established between the solutions of the stationary equation (16) and the dynamical system (17). However, the temporal stability of spatial solutions found using 17 can only be determined by studying the full system of PDEs using numerical simulations.
We use the continuation package AUTO to compute the localised patterns. 29 As result of the symmetry in the spatial system, we only need to find the solution on the half domain [0, L], where L ≫ 1 with homogeneous Neumann boundary condition. Generally, when displaying a bifurcation diagram in a single parameter, we use the L 2 norm: When required, we obtain a numerical solution of the PDE model (2) using an explicit finite difference method with time and space increments dt = 10 −3 and dx = 5 × 10 −2 , respectively. To determine the stability of a non-homogeneous solution, we apply the finite difference approximation method to solve the model obtained by linearising around the solution, i.e., system (2) about the nonhomogeneous steady state. Additionally, we used Matlab to find the spectrum of eigenvalues for the linear operator by discretising the non-homogeneous solution of (2).
In Sections 3 and 4, we implement the methodology outlined above to investigate the existence of localised structures in two predator-prey models of the form (2). These models have been chosen because neither of them exhibits the classical Turing instability.

THE IVLEV MODEL
We start in Section 3.1 by reviewing the Ivlev model and stating the homogeneous steady-state solutions. In Section 3.2, we systematically apply the methods outlined in Sections 2.3-2.4 to the Ivlev model (19). This culminates in the identification of a co-dimension two bifurcation at which the spatial instability changes from being subcritical to supercritical-a crucial ingredient in the identification of localised spatial solutions.

Model formulation
Holling, 30 based on field data, suggested that the predator-prey interaction should be given by a function that is monotonically increasing and uniformly bounded. Ivlev 31 [32][33][34][35][36][37][38] and has proven to be of great interest to both mathematicians and ecologists. Ecological applications include plankton population dynamics, 39 the interactions between a host and a parasite, 40 and integrated pest control management. 41 Mathematical interest in using the Ivlev model can be found in such work as Metz and Diekmann, 42 Sherratt et al 43 38 We investigate the existence of localised patterns in the Ivlevl predator-prey reaction-diffusion model. The scaled model equations for the one-dimensional model are as follows.
Rate of change of prey density (−L < x < L): Boundary conditions: (For the 2D model, we have x, ∈ (−L, L), subject to homogeneous Neumann boundary conditions along the boundaries, i.e., x = ±L and = ±L.) In Equation 19a, the term u (1 − u) indicates that, in the absence of predation, the prey undergoes logistic growth. The term v (1 − e − u ) represents the reduction in prey numbers due to predation. (This is the aforementioned Ivlev predation model.) In Equation 19d, the tern v (1 − e − u ) represents the growth in predators due to their consumption of prey. (It is assumed that the predators only increase their numbers through predation upon the single species of prey.) The term − v represents natural mortality of the predator.
In system (19), is the efficiency of predator capture of prey, is the death, rate of the predator, and is a measure of the conversion rate of prey captured by the predator.
A straightforward calculation shows that (19) has three homogeneous equilibrium states. These correspond to extinction of both prey and predator (0, 0); a predator free equilbrium (1, 0) and a coexistence state S + = (u * , v * ), where The positive coexistent steady state is meaningful when > 1 and > 0 > .

Bifurcation conditions and the origins of localised solutions
Having found the homogeneous steady-state solutions, we apply the linear stability analysis outlined in Section 2. For the Ivlev model (19), the general dispersion relation (9) takes the form where and .
Applying (10) to (21), the critical value ( H ) at which the coexistence homogeneous steady state solution undergoes a Hopf bifurcation is found to be Note that the value of H only depends upon the parameter . Applying the conditions for the generalised spatial instability (11), we find that the condition for a generalised spatial instability is with corresponding critical value for the wave number (k c ) given by The Belyckov-Devaney transition curve is found to satisfy the equation: Figure 1A shows dispersion curves for ℜ( ) as a function of the wavenumber k for three values of the parameter . If the value of is sufficiently low (blue dashed curve), the real part of the eigenvalue is always negative, and there is no spatial instability. If the value for is sufficiently high, then there is a range of values for the wavenumber (k) for which the eigenvalues are positive (red solid curve) and the homogeneous solution is unstable to spatial perturbations. At the critical value for , the dispersion curve touches the line ℜ( ) = 0 tangentially at the critical value for the wave number (k c ) (black dot curve).
In what follows, we fix the parameter values of and to be 1.1 and 0.06, respectively. We now apply weakly nonlinear analysis, as outlined in Section 2.4, to the Ivlev model (19). In particular, we obtain an algebraic expression for the coefficient C 3 . Figure 1B displays the value of this coefficient as a function of the single variable . (It should be noted that the condition for the spatial instability gives = ( ).) As the value of increases from zero, the coefficient C 3 is positive, and the spatial instability is subcritical. For sufficiently large values for , the coefficient is negative, and the spatial instability is now supercritical. The green dot indicates the critical value of , at = 6.034, at which the coefficient is equal to zero.
This point is codimension two, because a value of implies a particular value for . Its significance is that, due to the change in sign of C 3 , bio-stability between the homogeneous steady state and the pattern state emerges which is crucial element for the birth of localised structures in reaction-diffusion systems as stated in Brena Medina and Champneys. 14 Figure 2 is a bifurcation diagram for the Ivlev model (19) in the ( , ) plane. To be more particular, this is the bifurcation diagram associated with the coexistence homogeneous steady-state solution. Shown on the diagram are several boundaries which identify particular solution characteristics or behaviour. These boundaries have been obtained using results from Section 3.2.

The bifurcation diagram and localised structures solution patterns
The vertical red dashed line at = H is where the homogeneous coexistence solution losses stability at a Hopf bifurcation. (Recall that Equation 22 shows that H is independent of .) The homogeneous solution is stable when < H and unstable when > H .
The pink curve is the spatial instability boundary (23). This curve splits the bifurcation diagram into two spatial stability regions. From (23a), we deduce that solutions to the right are unstable whilst solutions to the left are stable.
The presence of a Hopf bifurcation adds intricacies to the behaviour of the solutions of the Ivlev model (19). Solutions to the right of the Hopf bifurcation line have temporal periodic mode. The codimension two bifurcation identified through Figure 1B; i.e., the point at which the spatial instability changes from being subcritical to supercritical is indicated by the green dot. This point is located at (6.034, 1.773). The spatial instability bifurcation is supercritical for > 1.773 when increases and is subcritical for < 1.773 when decreases. The localised structures, the boundary of which is denoted by the black curve, emanates from the codimension two point.
The black curve is obtained using the numerical continuation package AUTO. The Belyckov-Devaney curve (24) divides the region of localised structures into two. The sub-region to the right of the BD curve is called the snaking region. In this region, the spatial eigenvalues of the localised solutions are complex conjugates. As a consequence, the tail of the localised solutions decays in the far field in an oscillatory manner. The sub-region to the left of the BD curve is called the single hole region. In this region, the spatial eigenvalues of the localised solutions are purely real. As a consequence, the tail of the localised solutions decays monotonically in the far field.
The leading spatial eigenvalues of the purely spatial part of the solution of (19) are summarised in Table 1.
In the next section ,we focus on investigating non-homogeneous solutions by numerical means.

Numerical continuation
In this section, we apply the procedure outlined in Section 2.5 and transform the steady-state equations of (19) into a system of four first-order nonlinear ODEs. For the Ivlev mode, the formulism of (17) becomes (25) Figure 3 shows a one-parameter continuation of the localised solutions of (19) using (25) in (17). The homoclinic snaking diagram displays two branches of localised patterns. The purple branch represents solutions with an odd number of 'holes', whilst the brown branch represents solutions with an even number of 'holes'.
Over the parameter range 5.53 ≤ ≤ 5.745, the structure of the solution branches is evidently convoluted due to the sequence of saddle node bifurcations that wobble backwards and forwards. This sequence of saddle noddle bifurcation has led to this structure being called the snaking or pinning region. 19,48,49 These solution curves are born unstable at = 5.745. Consequently, at each left (right ) saddle node bifurcation, there is a temporal stability change whereby the solution branches of the snaking gains (loses) stability. The dashed red vertical line shown on Figure 3 is the Hopf bifurcation of the homogeneous solution, identified by Equation 22 and also shown on Figure 2. This bifurcation line cuts the snaking diagram into two, creating another instability for the snaking branches. The consequence of the Hopf bifurcation is that all snaking branches that are stable following the left fold lose their stability at the Hopf bifurcation line. (In other models, localised patterns can lose their stability before the Hopf bifurcation curve. 50 ) As described in Section 2.5, we obtain the stability of each localised solution by computing the spectrum of the discretised nonhomogenous linear operator in (19). The left-hand side of Figure 4 depicts examples of localised hole patterns from the odd, the bottom two profiles, and even, the top two profiles, branches using selected values from both the stable and unstable side of the Hopf bifurcation line. The corresponding plots on the right-hand side display the leading eigenvalues of each localised pattern.
Tracing the fold of one of these saddle node's in Figure 3 in ( , )-plane shows that the black curve in Figure 2 which envelops the localised structures region.
Results of the continuation from the single hole region, displayed in Figure 2, are depicted in Figure 5. Figure 5A shows the maximum amplitude of u (MAXu) against the continuation parameter . This continuation displays a single fold which is connected to two branches at the saddle node bifurcation point. These branches are stable (unstable) and are identified by a thick (thin) line. Immediately to the right of the continuation figure are two solutions. The solution from the stable branch has a larger amplitude, whereas the solution from the unstable branch has a smaller amplitude. The leading eigenvalues for each solution are shown. Note that the instability here is the result of the fold bifurcation. In Figure 5B, the continuation parameter is rather than . On this case, the stable branch losses stability when the continuation parameter crosses the Hopf bifurcation line.
To further study the dynamics and solution behaviour of system (19), we solved the model numerically-see Section 2.5 for further details. For the localised hole solutions, these numerical simulations allow us to explore the nature of the Hopf bifurcation, i.e., whether they are sub-or super-critical. During these investigations, we observed that the solution components u and v have their 'holes' at the same spatial co-ordinate. Accordingly, we only present simulation results for the u component.
To demonstrate the stability behaviour of localised patterns, a careful selection of the values for is important. From Figure 3, values of are chosen from either side of the Hopf bifurcation line ( H = 5.552). Figure 6 shows both solution heat maps (left plots) and the associated solution time series at a single point in the spatial domain (right plots). It is evident, from both the solution heat map and the time series, that when = 5.545, a value before the Hopf line, that the solution is locally asymptotically stable. Conversely, when = 5.65, a value directly after the Hopf bifurcation line, the solution to evolves into a stable oscillation; see Figure 6B. These simulations show that the Hopf bifurcation line is supercritical. This implies that there is a stable limit cycle just beyond the Hopf bifurcation line.

Two-dimensional spatial simulation results
To investigate the existence of localised structures in two spatial dimensions, the system (19) was simulated on a finite domain, (x, ) ∈ [0, L] × [0, L], with zero flux boundary conditions. A finite difference method was used with a time step dt = 0.001 and a spatial step size dx = dy = 0.1. The simulations were run until the steady state solution was reached. As in the 1D case, the parameters = 1.1 and = 0.06 were fixed, whilst the parameters and were varied.
It was found that both the prey and predator have the same spatial structure. Therefore, in the following, we only discuss the prey. Figure 7 shows solution heat maps depicting two stable localised patterns for the prey distribution. The parameter values ( , ) = (5.5, 0.25) represent typical values from the snaking region, whilst the parameter values ( , ) = (5.25, 0.1) represent typical values from the isolated hole region. Consequently, Figure 7A,B is heat maps for the prey density in the snaking and isolated hole regions, respectively.
A deduction from these figures is that the snaking region produces more holes. From a biological perspective, the localised holes (cold or blue spot) indicates regions where the prey density has been significantly reduced by the predator. In these regions, the predator density is higher than that of the prey. This is in contrast to the behaviour exhibited across the rest of the domain where the prey density is approximately double that of the predator. We now consider the unstable region that lies within the snaking-Hopf region in Figure 2. In this region, a careful selection of the parameter values for & can produce a solution that exhibits spatio-temporal behaviour that is affected by both the Hopf bifurcation and the localised holes region. Figure 8 shows a solution heat map showing the type of spatio-temporal behaviour that may occur within this region. Taking ( , ) = (5.768, 0.602) as typical values from the snaking-Hopf region, it was found that cold spot patterns occurred but are disrupted by the unstable homogeneous solution-see Figure 8A. Figure 8B displays the prey component time series at the point (x, ) = (25,25). The oscillatory time series show that the localised solution is unstable. The dynamics displayed in this figure show the interaction of the Hopf bifurcation with a solution with localised holes. The 'static' localised solution is unstable and is replaced by a 'dynamic' localised solution: The number and location of the holes do not change, but the value of the solutions is changing. These oscillations are also occurring in the background region shown in yellow. This type of behaviour has been described as 'competing' dynamics. 51,52

HOLLING II MODEL
In this section, we investigate a predator-prey model in which the Ivlev predator-prey functional response is replaced by the Holling II response. Our interest in this new system is whether it also exhibits a region of localised structures and, if so, to determine if any of the solutions within this region manifests spatio-temporal behaviour. In doing so, we make compare and contrast the behaviour of the Ivlev predator-prey system with the Holling II system. To prevent repetition, the reader is referred back to Section 2 and to related material from Section 3.

Model formulation
The Holling type II functional response 30,53 is a ratio dependent model having the form u∕ ( + u), where is a positive constant. This functionality originated in Holling's study of the predation of small mammals on European pine sawflies. The Holling II response has been used in predator-prey models and has been extensively studied; see, for example, previous studies 21, [54][55][56][57][58] and reference therein.
We investigate the existence of localised patterns in the Holling II predator-prey reaction-diffusion model. The model equations are given in system (26a).
Rate of change of scaled density of prey: Rate of change of scaled density of predator: with Neumman boundary conditions. As before, all parameters are positive. The parameter denotes the scaled prey-density at which the per-capita predation rate is half its maximum value. The parameter is the maximum per-capita growth rate of predators. The parameter is the per-capita death rate of predators. The parameter is the ratio of the diffusion constant for prey to that of predators. In the forgoing analysis, the parameter values = 0.7 and = 0.06 are used.
A straightforward calculations shows that system (26a) has three homogeneous equilibrium states. These correspond to extinction of both prey and predator S 1 = (0, 0); a predator free equilibrium, S 2 = (1, 0); and a coexistence state The positive coexistent steady state is meaningful when > (1 + ).
We focus our analysis on the coexistent steady state, S 3 .

Bifurcation conditions and the origins of localised solutions
Having found the homogeneous steady-state solutions, we apply the linear stability analysis outlined in Section 2 to the coexistence steady-state S 3 . Applying (10), the critical value ( H ) at which the coexistence homogeneous steady state solution undergoes a Hopf bifurcation is found to be Applying the conditions for the generalised spatial instability (11), we find that the condition for a generalised spatial instability is with corresponding critical value for the wave number (k c ) given by The Belyckov-Devaney transition curve is found to satisfy the equation

The bifurcation diagram and localised structures solution patterns
After applying linear and weakly nonlinear stability analysis, as outlined in Sections 2.3 and 2.4, the bifurcation diagram Figure 9 is obtained. This figure demonstrates that there are different types of regions which represent stable and unstable solution patterns. Referring to (3), the blue region indicates spatial instability. As in the Ivlev model, the region where localised structures are found, i.e., the yellow region, is born at the codimension two point. The boundary of the yellow region is found using AUTO. Another similarity to that of the Ivlev model is that the region of localised structures is divided into two (isolated hole and snaking) sub-regions by the BD curve, i.e., Equation 4. Furthermore, the Hopf bifurcation (red dashed) curve divides the ( , )-plane into two sub-regions, in the process producing two contrasting temporal solutions. That is, the left (right) region is where stable (unstable) temporal solutions exit. It is to be noted that both the Ivlev and the Holling II responses models present similar qualitatively, but not quantitative, bifurcation diagrams. Figure 10 displays continuation bifurcation diagrams from the snaking region, Figure 10A, and from the isolated hole region, Figure 10B. The equivalent figures for the Ivlev model are Figures 3 and 5.

Numerical continuation
We first discuss Figure 10A. On the left-hand side of the figure, we see a sequence of limit-point bifurcations for values of a little bit below ≈ 0.2. There is a corresponding sequence of limit-point bifurcations for values of at a little higher than ≈ 0.21. Pick one of the limit-point bifurcations on the left-hand side and select the upper branch that emanates from the limit-point. Follow these solutions until they reach the corresponding limit-point on the right-hand side. Along this traverse of the snaking diagram, the solution is unstable; from here, the solution is stable as decreases until the Hopf bifurcation line = H is crossed, whereafter the solution is again unstable. Figure 10B shows the continuation diagram from the isolated hole region. Note that in this case, the localised solution is stabilised as the value of is increased through the Hopf bifurcation point at = H . This figure is similar to the Figure 5B in Ivlev model.
The subfigures to the immediate right of the bifurcation diagrams show two solution profiles and their associated leading eigenvalues. The values of used in these subfigures have been chosen to have one solution before and one solution after the Hopf bifurcation line.
In Figure 10A,B, the location of the spatial structures on the corresponding bifurcation diagram are shown by the black and green coloured dots. The former indicates a stable solution, whereas the latter indicates an unstable solution.
It can be observed that the solution behaviours displayed in Figure 10 are similar to that of the Ivlev system-see

CONCLUDING REMARKS
In Section 2, we reviewed the general methodology that is used to investigate both the homogeneous and non-homogeneous solutions of a general system of nonlinear PDEs. Our main contribution in this section is to bring to the attention of a wider readership the fact that the Turing bifurcation is a special case of a more general mechanism to generate spatial patterns. In particular, a system that fails to undergo a Turning instability can still exhibit spatially localised patterns. The work of Cross and Hohenberg 22 is not as well known as it deserves. In Section 2.5, we outlined the computational advantages that accrue when the steady-state problem is rewritten a system of four non-linear ordinary differential equation in which the spatial coordinator is now viewed as a 'time-like' coordinate. Continuation methods can now be applied, which allows the parameter regions over which certain localised structures occur to be readily determined.
We have illustrated the use of our methods by applying them to two predator-prey type models with zero Neumann boundary conditions. These models, the Ivlev and Holling type II systems, considered in Sections 3 and 4, respectively, were picked because it has been established by previous authors that they do not undergo the classical Turing bifurcation. 20,21 The novelty of our analysis is that we establish that both systems exhibit a range of localised spatial patterns despite the fact that neither system undergoes a Turing instability. In both cases, we determined theoretically that localised patterns exist before carrying out a numerical investigation. Pattern formation was explored in a limit in which the prey population diffuses much more slowly than the predator population.
For each of these systems, a universal two-parameter space bifurcation diagram was generated. Both of these diagrams are similar in nature: Both included a codimension two point indicating where the spatial instability bifurcation changes from being subcritical to supercritical. The region of localised spatial structures originates at the codimension two points for each of the systems. In both models, the Belyckov-Devaney transition curve divides the region of localised solutions into two. These sub-regions are the isolated hole region and the snaking region.
Within the isolated hole region, the localised structures exhibit a single hole pattern with a tail that decays monotonically: The spatial eigenvalues are purely real. In the snaking region created, the localised solutions have an oscillatory decaying tail: The leading spatial eigenvalues are complex conjugates.
In both models, the Hopf bifurcation curve for the coexistent steady state cuts through the region of localised solutions. This curve adds additional complexity via the creation of instabilities in the localised solutions. That is, temporal periodic localised solutions are found at parameter values near to the Hopf bifurcation curve, with a stable limit cycle occurring just beyond the Hopf bifurcation.
The similarity in the bifurcation behaviour of the two models suggests that a deeper theoretical analysis is warranted. It appears that the precise functionality of the predator-prey response does not determine the underlying bifurcational structures. Rather, it is the generic shape of the functional response that is important.
Finally, we have shown that the onset of intricate spatio-temporal patterns in the 2D Ivlev system can be predicted from the theory developed for the 1D model.