SPATIOTEMPORAL PATTERN FORMATION IN A PREY–PREDATOR MODEL WITH GENERALIST PREDATOR

. Generalist predators exploit multiple food sources and it is economical for them to reduce predation pressure on a particular prey species when their density level becomes comparatively less. As a result, a prey-predator system tends to become more stable in the presence of a generalist predator. In this article, we investigate the roles of both the diﬀusion and nonlocal prey consumption in shaping the population distributions for interacting generalist predator and its focal prey species. In this regard, we ﬁrst derive the conditions associated with Turing instability through linear analysis. Then, we perform a weakly nonlinear analysis and derive a cubic Stuart-Landau equation governing amplitude of the resulting patterns near Turing bifurcation boundary. Further, we present a wide variety of numerical simulations to corroborate our analytical ﬁndings as well as to illustrate some other complex spatiotemporal dynamics. Interestingly, our study reveals the existence of traveling wave solutions connecting two spatially homogeneous coexistence steady states in Turing domain under the inﬂuence of temporal bistability phenomenon. Also, our investigation shows that nonlocal prey consumption acts as a stabilizing force for the system dynamics.


Introduction
Mathematical modeling plays an important role in understanding various complex phenomena arising in real world ecological systems.In this regard, prey-predator systems are commonly used by researchers to model intra-and inter-specific interactions among different species in an ecological system.Traditionally, the predators are classified into two categories such as specialist and generalist.A specialist predator depends exclusively on a focal prey species for food, and as a result, it goes to extinction in the absence of that particular prey species.On the other hand, a generalist predator feeds on a variety of prey species and can eventually persist in an ecosystem in the absence of one particular prey species.In mathematical models, the distinction between both the types of predators has been generally achieved through the incorporation of appropriate functional responses.Several authors have considered a sigmoidal functional response in the prey-predator systems to characterize the predation by a generalist predator.This type of functional response efficiently captures the prey-switching phenomenon which states that a generalist predator should stop foraging for a specific prey species with low concentration and instead keep focusing on other available prey [37].However, the models with a sigmoidal functional response and without other food sources for predator are somewhat inconsistent with the generalist nature of the predator as it can not persist in the absence of the considered focal prey species [16].
Only a few studies have considered additional food sources for generalist predator along with the sigmoidal functional response.For example, Spencer and Collie [50] incorporated this aspect of a generalist predator in the modeling framework through a logistic growth function which arises due to the inclusion of both the linear growth and quadratic death terms for the predator.They largely concentrated on the emergence of three coexistence equilibria and the resulting bistability phenomenon.van Baalen et al. [54] studied the importance of switching behavior of predators on the persistence of prey-predator systems and suggested that the numerical response of the predator population should exclusively incorporate the contribution of alternative food.Further, van Leeuwen et al. [55] considered a two-prey-one-predator system with a functional response depending on both the prey densities and investigated the conditions under which the sigmoidal Holling type III functional response would be apparent.Their study indicated that the constant density of one prey population is key for the emergence of the sigmoidal functional response.Also, a study regarding the biological control of leafminers by generalist parasitoids was conducted in [29] by considering logistic growth for both the species with a Holling type II functional response.Finally, Erbach et al. [16] investigated generalist predator-prey dynamics with a Holling type III functional response for predation on the focal prey and a Beverton-Holt-like function for predator reproduction in the absence of the focal prey.This study elucidated an interesting bistability scenario between a stable coexistence steady state and a stable population cycle along with the usual bistability between two stable coexistence steady states.
Spatial diffusion is an intrinsic characteristic for almost all ecological species and it can efficiently explain the high degree of bio-diversity observed in nature by producing spatially heterogeneous population distributions.Various field and experimental studies in ecology assert the claim of spatial heterogeneity induced species diversity.For example, Gause [17] showed the important role played by spatial heterogeneity in long term survival of a species through a laboratory experiment on paramecium and didinium.Performing a series of laboratory experiments on two species of mite such as Eotetranychus sexmaculatus and Typhlodromus occidentalis, Huffaker [19] demonstrated the long term coexistence of both the species in a heterogeneous environment, whereas, they could become extinct in a perfectly homogeneous set up.Experiments conducted by Luckinbill [26,27] emphasized the importance of dispersal on persistence and coexistence for interacting populations.The spatially heterogeneous population distributions in nature are often studied theoretically by using the reaction-diffusion modeling framework.Being inspired by the seminal work of Turing on chemical morphogenesis [52], Segel and Jackson [48] first employed this modeling approach to explain the dissipative structures observed in ecological communities.Further, this modeling approach was successfully applied to elucidate the plankton patchiness, vegetation patterns in semiarid region and invasion by an exotic species [21,24,34,42,49].Both the Turing and Hopf-Turing bifurcations are the two widely accepted mechanisms behind the emergence of spatially heterogeneous community structures in ecology [7,30,43].Due to these mechanisms, a spatially extended ecological system can induce a variety of stationary and dynamic patterns [7,30,43].Stationary patterns include labyrinthine, spot (both the cold and hot spots), and mixture of stripes and spots, while periodic, quasi-periodic and chaotic patterns are some examples of dynamic patterns [7,30,43].Also, the reaction-diffusion systems can efficiently explain the invasion through continuous population fronts propagation and the patchy invasion by an exotic species observed in nature [30,36,42].
The theoretical and numerical investigations on the pattern formation for spatially extended prey-predator systems largely concentrated on the specialist predators.On the other hand, only a limited number of studies considered the generalist predators in their modeling framework and examined the subsequent prey-predator dynamics.For example, Magal et al. [29] showed that generalist parasitoids mediated control of leafminers can be achieved only if leafminers disperse quite rapidly.They identified the dilution effect arising due to rapid spreading as a reason behind this counter-intuitive outcome.Kumari [23] studied a diffusive tri-trophic food chain model for marine ecosystem, and found that top predator being generalist in nature can lead to wave of chaos phenomenon and eventually destabilize the ecosystem.Further, Chakraborty [13] reported a wide class of stationary (spots, stripes and mixture of them) and chaotic patterns for a spatial model accounting for interactions between a generalist predator and its focal prey.Recently, we considered a reaction-diffusion system for two preys and their common predator, and studied the detailed bifurcation scenario [32].In this study, we found the existence of multiple stationary states and a long quasi-stable transient dynamics.Also, Rodrigues et al. [45] explored the dynamics of the spatial version of the model presented in [16] within Turing and Hopf-Turing domains, and observed the existence of either stationary or spatially homogeneous periodic population distributions.They also suggested that spatial dispersal can effectively suppress the hysteresis effect observed in the corresponding temporal dynamics.Finally, Arancibia-Ibarra et al. [2] investigated the complex pattern formation scenario of a modified Holling-Tanner model with an alternative food source for the predator, and identified that the inclusion of alternative food source can induce more accurate predictions regarding the spatial distributions of weasels and ermine populations in the Boreal Forest and Foothills National Regions.
Generally, the reaction-diffusion systems accounting for prey-predator dynamics assume the pointwise intraand inter-specific interactions where an individual interacts with others present at the same spatial location.Nevertheless, interactions depend on the average population density over a certain area adjacent to individuals' present location in reality due to spatial dispersal [5,11,33].This assumption leads to a system of integro-partial differential equations where a convolution integral with an appropriate kernel function takes care of this type of nonlocal interactions [5,11,33].In [8,14], it was manifested through some biological experiments that toxic substances or nutrients released by bacteria cultures in Petri dishes can effectively induce nonlocal interactions.Also, nonlocal interactions arise in the case where ladybirds (Hippodamia converges) feeding on aphids (Aphis varians) [12].Recently, there has been growing interest among researchers to model prey-predator interactions incorporating nonlocality.However, only a very few articles have dealt with nonlocal generalist predator-prey models.In this direction, Autry et al. [3] studied a nonlocal three-species food chain model in the context of biological control where the top predator is assumed to feed on both the prey (an important crop) and intermediate predator (a pest).They showed that it is impossible to control a rapidly dispersing pest species, however, introduction of a highly diffusive top predator can effectively restore the control.They also identified that nonlocal behavior of the top predator acts as an important factor to induce a robust partial control.Further, Han et al. [18] investigated a spatial intraguild predation model with a nonlocal interaction term in the growth of the shared resource and found the pattern transition from stationary to chaotic with the increasing intensity of nonlocal interaction.On the other hand, we recently showed that chaotic population oscillations can be significantly suppressed by incorporating sufficiently large extent of nonlocality, and hence, nonlocal interactions can also act as a stabilizing force for a two-prey-one-predator model [32].
In this article, we revisit the prey-generalist predator system proposed by Erbach et al. [16], and systematically explore the implications of spatial dispersal and nonlocal prey consumption on the population distributions by means of Turing and Hopf-Turing bifurcations.The rest of this article is arranged in the following manner.In the next section, we describe the non-spatial system and illustrate some interesting aspects of its dynamics.In Section 3, we introduce the corresponding spatial system and present a wide variety of resulting spatiotemporal dynamics.In order to elucidate some quantitative information regarding the emerging stationary Turing patterns near the bifurcation boundary, we derive Stuart-Landau equation for amplitude of the patterns through weakly nonlinear analysis in Section 3.1.Further, we extend the spatial system by taking into account nonlocal prey consumption and explore the influence of the extent of nonlocality on the resulting population distribution in Section 4. To end this article, we present a brief discussion in Section 5.

Non-spatial model
In this section, we describe the non-spatial prey-predator model with a generalist predator studied in [16].Also, we present a wide spectrum of interesting temporal dynamics possessed by this model in a concise manner as the prior knowledge of these dynamics can effectively help us in understanding the corresponding spatiotemporal dynamics.The concerned non-spatial model reads as Here, N stands for the density of a focal prey species, whereas the density of a generalist predator of it is denoted by P .The parameters r and K represent the intrinsic growth rate and the environmental carrying capacity of the prey species, respectively.The equation (2.1) shows that the prey species follows the logistic growth in the absence of the generalist predator and predation follows a Holling type III functional response.
In this functional response, the maximum predation rate and the handling time are respectively represented by α/h and h.The predation on the prey species in turn contributes to the predator density with γ being a conversion efficiency.The sigmoidal shape of the Holling type III functional response is usually recognized to contemplate the prey-switching phenomenon which suggests that it is prudent for the predator to stop foraging for a particular prey species with low density and instead start looking for available alternative resources [37].
To talk about few examples of generalist predators exhibiting Holling type III functional response, the mirid predator M. pygmaeus shows a type III response in consuming Trialeurodes vaporariorum (Westwood) [15].Also, the mirid predators E. varians and M. basicornis exhibit a type III response in consuming the eggs of T. absoluta [56].The generalist predator reproduces with a Beverton-Holt-like function in the absence of the focal prey.The parameters β, Q and m denote the per capita reproduction rate, the strength of density-dependence and the mortality rate for the predator population, respectively.Hence, the equation (2.2) shows that the predator species follows a logistic-like growth in the absence of the focal prey.All the parameters incorporated in the model (2.1)-(2.2) are assumed to be positive from ecological point of view.However, the inequality β > m acts as a necessary condition for predator persistence in the absence of the focal prey [16].Therefore, the dynamics of the generalist predator is modeled by using a sigmoidal response function and a separate growth term combining all other available food sources.In this case, the alternative growth term can also be viewed as predator interference.
We can observe that the model (2.1)-(2.2) depends on eight parameters.In order to reduce the number of parameters, we non-dimensionalize the model (2.1)-(2.2) by setting N = N/K, P = (αKP )/r and t = rt.With this change of variables and omitting tildes, we obtain the following non-dimensionalized system where a = hK 2 , b = (αγK 2 )/r, c = β/r, d = (rQ)/(αK) and e = m/r.In this case, we have c > e from the condition β > m.The system (2.3)-(2.4) is subjected to non-negative initial conditions.The model (2.3)-(2.4)always admits the total extinction steady state E 0 = (0, 0), the prey-only steady state E 1 = (1, 0) and the predator-only steady state E 2 = 0, c−e de .A typical coexistence steady state E * = (N * , P * ) of the system can be obtained by solving the following two algebraic equations simultaneously for positive solutions: ) From the equation (2.6), we can easily observe that any feasible positive value of N * corresponds to a positive P * as N * < 1.However, it is almost impossible to find an analytical expression for a feasible solution N * as we need to deal with a quintic polynomial (that is, equation (2.5)).Therefore, we need to depend on numerical computations to obtain a feasible coexistence steady state.The paper [16] demonstrated that the model always admits at least one coexistence steady state and the number of coexistence steady states can vary from one to three depending on the parameter values.In this study, we are particularly interested on the dynamics of the system around the coexistence steady states as they generally indicate the well-functioning of an ecosystem.It is easy to obtain that the extinction steady state E 0 is an unstable node, and both the axial equilibria E 1 and E 2 are saddle points always.Now, we derive the conditions for local stability of a typical coexistence steady state E * = (N * , P * ) by linearizing the system (2.3)-(2.4)about it.The corresponding Jacobian matrix is given by where Then, the characteristic equation is given by where tr(J ) = J 11 + J 22 and det(J ) = J 11 J 22 − J 12 J 21 .The coexistence steady state E * = (N * , P * ) is locally asymptotically stable when the following Routh-Hurwitz criteria are satisfied [22]: The stable coexistence steady state E * = (N * , P * ) loses its stability due to the Hopf bifurcation when a pair of complex conjugate eigenvalues passes through the imaginary axis.If we consider b as the bifurcation parameter, then the system (2.3)-(2.4)undergoes Hopf bifurcation at b H when the following conditions are satisfied simultaneously [22]: (2.10)However, it is not possible to obtain an analytical expression for the Hopf bifurcation threshold b H as the value of N * depends on the parameter b (see Eq. (2.5)).Therefore, we need to rely on numerical computations of b H in this case also.Note that various types of local and global bifurcation results can be obtained for the model under consideration, but such detailed results are beyond the scope of this study.Nevertheless, we demonstrate some interesting bifurcation scenarios such as super-/sub-critical Hopf, saddle-node, homoclinic bifurcations numerically in the following subsection.

Numerical results
In this section, we demonstrate numerically a wide variety of interesting dynamics possessed by the system (2.3)-(2.4) in a concise manner.Firstly, we consider the parameter values a = 100, c = 1, d = 0.9 and e = 0.5 [16].For this particular parameter set along with any feasible value of the parameter b, the system (2.3)-(2.4)admits a unique coexistence steady state.For b ∈ [40,80], the dynamics of the system (2.3)-(2.4)are encapsulated in Figure 1 as a bifurcation diagram.From Figure 1, we can observe that the coexistence steady state loses its stability and the stable oscillatory coexistence arises through a super-critical Hopf bifurcation at b 1 H = 48.116.Further, the stable oscillatory coexistence disappears through another super-critical Hopf bifurcation at b 2 H = 73.435and the coexistence steady state regains its stability.Note that the parameter b is proportional to the predation rate α.Thus, the dynamics illustrated in Figure 1 can be described in terms of the following ecological perspectives: low predation rate enforces the stable stationary coexistence, while increased predation rate induces stable oscillatory coexistence.However, sufficiently large predation rate can effectively diminish the population cycle, and consequently, the stationary coexistence becomes stable again with low prey density.Such a stable coexistence at higher predation rate is a characteristic feature of prey-predator models with generalist predators.
Next, we consider the parameter values a = 100, c = 1, d = 0.1 and e = 0.5 [16].Figure 2 elucidates the dynamics of the system (2.3)-(2.4)for b ∈ [20,26].In this case, a unique stable coexistence steady state bifurcates into three coexistence steady states through a saddle-node bifurcation at b 1 SN = 20.7807.Further, these three steady states disappear through another saddle-node bifurcation at b 2 SN = 23.2997 and only one coexistence steady state with low prey density persists for higher values of the parameter b.In the case of three coexistence steady states, the state with higher prey density always remains stable, while the other state with intermediate prey density remains as a saddle point.On the other hand, the unstable steady state with lower prey density becomes stable through a sub-critical Hopf bifurcation at b H = 20.9176 and an unstable limit cycle emerges around it.The unstable limit cycle then disappears through a homoclinic bifurcation at b HM = 21.014.From this figure, we can observe that the system experiences bistable situation for b ∈ (20.9276, 23.2997) where two coexistence steady states become stable.In this case, the long-time dynamics of the system is determined by the initial population densities.In terms of ecological point of view, it indicates that interventions like harvesting or replenishing can poke the system to a coexistence state from the other one.

Spatial model
In this section, we extend the model (2.3)-(2.4)by incorporating random dispersal of both the species and investigate the resulting spatiotemporal dynamics.For the sake of simplicity, we consider the one-dimensional spatial domain Ω = [−L, L] with L(> 0) ∈ R. The corresponding spatial prey-predator system is given by subjected to non-negative initial conditions: and periodic boundary conditions: Here, N (x, t) and P (x, t) denote the densities of prey and predator, respectively, at spatial location x ∈ Ω and time instant t.Note that we have assumed the spatial domain Ω to be one-dimensional for the sake of simplicity.The parameters D N and D P represent the diffusion coefficients of both the species, and both these coefficients are taken to be positive for ecological feasibility.All other kinetic parameters are described in Section 2. Recently, Rodrigues et al. [45] studied the system (3.1)-(3.2) in two-dimensional spatial domain with no-flux boundary conditions and paid special attention to the pattern formation scenario in Turing and Hopf-Turing domains.However, their study seems to leave out few crucial dynamical aspects.In this section, we will revisit the spatiotemporal dynamics of the model with one-dimensional space and periodic boundary conditions in a more complete manner.Note that the consideration of periodic boundary conditions is primarily due to mathematical simplicity.From the ecological point of view, periodic boundary conditions can take care of both the inflow and outflow of populations through the boundary, and thus, it seems more realistic in the case where the chosen spatial domain is just a part of the entire habitat.Of course in a compromised way, we may assume that the amount and rate of inflow from one end is same as the amount and rate of outflow from other end.The spatially homogeneous steady states of the system (3.1)-(3.2) are the steady states of the corresponding non-spatial system (2.3)-(2.4),and accordingly, we use the same notations to denote them for the sake of convenience.In order to study the local stability of a typical spatially homogeneous coexistence steady state E * = (N * , P * ) and subsequent bifurcations, we perturb the system (3.1)-(3.2) with small amplitude spatiotemporal perturbation about it and then linearize the system.The corresponding Jacobian matrix is given by where k represents the wavenumber and the quantities J ij (i, j ∈ {1, 2}) are described in the previous section.
Then the characteristic equation is given by where Thus, the spatially homogeneous coexistence steady state E * = (N * , P * ) is locally asymptotically stable if the following Routh-Hurwitz criteria are satisfied for all k > 0: Also, we can easily observe that tr(J (k)) < 0 for all k > 0 when E * = (N * , P * ) is locally stable in the corresponding non-spatial system.Now, the system (3.1)-(3.2) can undergo the Turing instability when E * = (N * , P * ) is stable under any homogeneous perturbation, however, it becomes unstable with respect to spatially heterogeneous perturbations [7,38,52].Thus, we need to have det(J (k)) < 0 for some k > 0 in order to induce Turing instability.Differentiating det(J (k)) with respect to k 2 and then setting the derivative to zero, we obtain the following critical value For a feasible k 2 crit , we have the necessary condition D P J 11 + D N J 22 > 0 which implies that the quantities J 11 and J 22 must have opposite signs.As we have J 22 < 0, then we need J 11 > 0 to satisfy this necessary condition.At this critical value, det(J (k)) becomes minimum and the minimum value is given by det Further, the minimum value (3.7) becomes negative when The condition (3.8) serves as a sufficient condition to induce Turing instability.Therefore, the necessary and sufficient conditions for Turing instability are summarized as follows [7,38,52]: The condition (ii) further implies that the quantities J 12 and J 21 must also have opposite signs which is actually true for our considered system (in this case, J 12 < 0 and J 21 > 0).The above discussion also leads to the following result regarding a condition under which the system (3.1)-(3.2) does not undergo the Turing or diffusion-driven instability.

Weakly nonlinear analysis
The above linear analysis regarding the Turing instability and the emergence of Turing patterns is qualitative in nature, and it fails to capture any quantitative information of the resulting patterns.On the other hand, weakly nonlinear analysis acts as an efficient method to quantify the emerging patterns sufficiently close to the Turing bifurcation threshold [10,20,25,28,39].This method enables us to derive an evolution equation for the amplitude of the resulting pattern through the technique of multiple-scales, since the pattern evolves on a slow time-scale for parameter values close to the bifurcation threshold [10,20,25,28,39].In this section, we aim to perform a weakly nonlinear analysis to derive the amplitude equation of the pattern for the spatial system where Here, N − N * and P − P * are respectively denoted by N and P for the sake of notational convenience.Note that a 10 = J 11 , a 01 = J 12 , b 10 = J 21 and b 01 = J 22 .Let U = (N, P ).Then the above system (3.9)-(3.10)can be rewritten in the following matrix form: where Let D * P denote the threshold value for the diffusion coefficient of predator population to induce the Turing instability.We also introduce the perturbation parameter ( 1), and then expand the bifurcation parameter D P , the time t and U as follows: P + 2 D (2) P + O( 4), (3.12) with U i = (N i , P i ) for i = 1, 2, 3. Then by simple computations, we obtain where Further, we have where . Now, substituting the above expansions in equation (3.11) and collecting the terms of O( ) to O( 3 ), we obtain the following equations, respectively: where The solution of the linear homogeneous equation (3.15) with periodic boundary conditions is represented by where . Now, let us denote L † * as the adjoint operator of L * .Then, the solution of the corresponding adjoint system of (3.15) is given by where ψ = (g, 1) ∈ Ker(L † * ) with g = b10 D N k 2 * −a10 .In order to satisfy the Fredholm solvability condition for equation (3.16), that is S 2 , U 1 = 0, we need to take D (1) P = 0 and T 1 = 0 [10,20,25,28,39].Then, we obtain where Thus, the solution to equation (3.16) can be given by where Again, substituting U 1 and U 2 into equation (3.17) we obtain S 3 as 3 A + S (3) where S Now, imposing the Fredholm solvability condition for the equation (3.17), we obtain S 3 , U 1 = 0 which results in the following cubic Stuart-Landau equation for the amplitude A [10,20,25,28,39]: where 3 , ψ ψ, ψ and 3 = S It is well-known that the growth rate coefficient 1 is always greater than zero in the Turing unstable region [10,28].Thus, the pattern formation scenario of the spatial system (3.1)-(3.2) close to the bifurcation threshold depends on the sign of the Landau constant 3 which determines the stability property of the cubic Stuart-Landau equation (3.22).If 3 > 0, the Stuart-Landau equation (3.22) admits a stable equilibrium solution , and consequently, the instability becomes supercritical in nature [10,28].In this case, the longtime behavior of the solution of the spatial system (3.1)-(3.2) is given by U = ψA ∞ {cos(k * x) + sin(k * x)}.On the other hand, the equation (3.22) does not admit any stable equilibrium solution for 3 < 0 which leads to subcritical instability [10,28].In this case, higher order terms (particularly, terms up to O( 5)) must be taken into account for the weakly nonlinear analysis and we should recover the quintic Stuart-Landau equation to determine the actual long-time behavior [10,28].However, it is not possible to predict the sign of the Landau constant 3 from the analysis itself, and hence, we need to rely on numerical computations.In order to numerically integrate the spatial system (3.1)-(3.2),we used finite difference approximation of the system by employing central difference and forward Euler schemes for diffusion and reaction parts, respectively.Both the temporal and spatial grid sizes have been chosen adequately so that numerical artifacts could not arise.Furthermore, numerical simulations have been carried out over a one-dimensional spatial domain [−100, 100] with the periodic boundary conditions and the following pulse-type initial conditions:

Numerical results
where |ε N |, |ε P | 1.Now, we illustrate various stationary and dynamic patterns for some representative parameter sets from Turing, Hopf and Hopf-Turing domains.Figure 4(a) represents the stationary Turing pattern for the prey species.As the corresponding pattern for the predator species is qualitatively similar, we chose not to display it here.The spatial density distributions for both the species are illustrated in Figure 4(b) which shows that patches with high prey aggregation attract the predator individuals and, in turn, high density patches for both the species coincide.Further, Figure 5 demonstrates two different dynamic patterns such as quasi-periodic and periodic patterns for parameter sets from the Hopf domain.Note that the periodic pattern illustrated in Figure 5(b) is homogeneous in space.The corresponding phase portraits of spatially averaged densities for both the species are illustrated in Figure 6.Extensive numerical simulations suggest that any parameter set from the Hopf domain can eventually lead to any one of these two dynamic patterns.Finally, we have observed the emergence of stationary and periodic patterns (similar to the patterns illustrated in Figures 4(a) and 5(b)) for the parameter sets chosen from the Hopf-Turing domain and accordingly, we have avoided to include them in this study.In this case, we have found that the periodic pattern emerges for the parameter values close to the Turing curve and the pattern becomes stationary when the parameter values are taken sufficiently away from this curve.with low prey density.We would like to mention here that the illustrated traveling wave solution is independent of the choice of boundary conditions as we stopped our simulation before the wave could hit the boundaries [31].For larger values of D P , the traveling wave front persists with more oscillatory transition zone.Similar situation arises when we choose the parameter values from the Turing domain for E (1) * and b lies on the left of the homoclinic bifurcation threshold b HM = 21.014.However, the situation becomes more complex when we cross the threshold b HM in the Turing domain for E (1) * .In this case, the emerging traveling wave front disappears, and eventually, a stationary heterogeneous structure arises for sufficiently large value of D P (for example, see the Figure 10 for b = 22).Note that the Turing threshold for E (1) * in this case is given by D T P = 1.5322.Figure 10(a) illustrates the stationary spatial distributions for both the species with D P = 20, while Figure 10(b) demonstrates the temporal trajectories of the spatially averaged density of the prey species (< N >) for different values of D P from the Turing domain.Figure 10(b) clearly shows how spatial dispersal enriches the overall population level of a species.In this figure, the black dashed lines represent the lower and upper branches of the prey density at the coexistence steady states, respectively.The corresponding temporal evolution for the predator species is similar, and hence, we did not present them here.Further, we do not observe any traveling wave front when we move into the Turing domain where both E

Nonlocal model
In this section, we explore the effects of nonlocal prey consumption on the resulting spatiotemporal dynamics of the system (3.1)-(3.2) which assumes that individuals interact locally.However, any dispersing individual can also interact with other individuals situated at adjacent locations [5,32,33].Thus, interactions among individuals should depend on weighted average of the population densities in a suitable region around their current position [5,32,33].In particular, prey population situated at the spatial location x can effectively be consumed by their generalist predators situated in some adjacent neighborhood, and in this case, the consumption rate is regulated by the integral ∞ −∞ φ(x − y)P (y, t)dy where φ(x − y) signifies the efficiency of a generalist predator located at the spatial position y in consuming its focal prey located at the spatial position x [5].Consequently, the kernel function φ should depend on the distance between x and y, and hence, it is considered to be a function of the difference x − y.Also, predator's capacity of consuming prey is actually limited in reality and this limitation can be applied to each spatial location independently [5].Thus, the functional response at each spatial location x in this case takes the form: On the other hand, numerical response of generalist predators at the spatial location x is proportional to their density at this location and to their focal prey population accessible around them [5].As the predators situated at the spatial location x can consume their focal preys available in some adjacent neighborhood, then the quantity of consumed focal preys is determined by the integral 1+aN 2 (y,t) dy [5].Thus, the numerical response at each spatial location x in this case takes the form: Therefore, the prey-predator model with a generalist predator and nonlocal prey consumption is given by supplemented with non-negative initial and periodic boundary conditions.For simplicity, we assume the kernel function φ as the following top-hat kernel function: This type of kernel function efficiently captures the unbiased dispersal of individuals and is non-negative.Further, it is in the normalized form as ∞ −∞ φ(z)dz = 1.Apart from mathematical simplicity, another important motivation behind the consideration of periodic boundary conditions is that the spatially uniform steady states for the spatial model correspond to the steady states for the associated nonlocal model.In the case of no-flux boundary conditions, we need to redefine the chosen kernel function in the vicinity of the boundary depending upon the extent of nonlocality [40].
Note that the steady states of the non-spatial system (2.3)-(2.4)correspond to the spatially homogeneous steady states of the nonlocal model (4.1)-(4.2) and accordingly, we denote them by the same notations.Now, we use the perturbations N = N * + ε N e ηt+ιkx and P = P * + ε P e ηt+ιkx with | ε | 1 to linearize the system (4.1)-(4.2) around a typical spatially homogeneous coexistence steady state E * = (N * , P * ).The Jacobian matrix of the corresponding linearized system is given by Thus, the corresponding characteristic equation reads as where By employing the Routh-Hurwitz criteria, we arrive at the following conditions accounting for the local asymptotic stability of E * = (N * , P * ): tr(J (k, M )) < 0 and det(J (k, M )) > 0. ( We can easily observe that tr(J (k, M )) < 0 holds true for all feasible k and M if E * = (N * , P * ) is locally asymptotically stable in the corresponding non-spatial model (2.3)-(2.4).Therefore, the system (4.1)-(4.2) undergoes the Turing instability only if det(J (k, M )) becomes negative for some k > 0 and M > 0. Thus, in order to obtain the Turing bifurcation threshold we solve the following two equations simultaneously [32,33]: Solving the first equation of (4.6) for D P , we obtain The above equation (4.9) presents a transcendental equation of the wavenumber k when other parameters are fixed.Thus, a solution of of the equation (4.9) represents a critical wavenumber k T and then by substituting it in the equation (4.9) we obtain the threshold value D T P of the diffusion coefficient of predator species to induce Turing instability.However, it is not possible to evaluate D T P analytically, and hence, we need to rely on numerical computations to obtain this threshold value.Finally, the nonlocal model (4.1)-(4.2) does not undergo the Turing instability when J 11 < 0 as it would imply det(J (k, M )) > 0 for all possible k and M .

Numerical results
In this section, we illustrate the spatiotemporal dynamics of the nonlocal system (4.1)-(4.2).Especially, we investigate the effects of nonlocal prey consumption on the resulting population distributions of both the species.For this purpose, we first consider the parameter values a = 100, c = 1, d = 0.9, e = 0.5 and D N = 0.1 as fixed, and vary the remaining parameters b, D P and M .Figure 11 illustrates the effect of nonlocal prey consumption on the location of Turing curve.It shows that the critical value of D P is inversely correlated with the extent of nonlocality.Further, we would like to mention that introduction of nonlocality in prey consumption cannot be able to alter the fate of regions indicated by light blue color (in those regions diffusive instability is not possible as the corresponding J 11 terms are negative).Now, we present some interesting spatiotemporal patterns of the nonlocal system (4.1)-(4.2) by simulating the model numerically.We integrated the system (4.1)-(4.2) numerically by using central difference, forward Euler and trapezoidal schemes for diffusion, reaction and nonlocal prey consumption terms, respectively.Further, we carried out the simulations over an one-dimensional spatial domain [−100, 100] with the periodic boundary conditions and the pulse-type initial conditions specified in (3.24).Also, we chose appropriate temporal and spatial grid sizes to discard any sorts of numerical artifacts.
Our numerical simulations suggest that the nonlocal prey consumption acts as a stabilizing force for the population distributions as a dynamic pattern can eventually become a stationary one in the presence of a sufficiently large extent of nonlocality (M ).Also, we have observed that a sufficiently large M can induce stationary population distributions when D N ≥ D P which is not possible in the absence of nonlocal prey consumption (that is, for the system (3.1)-(3.2)).In order to illustrate one such instance, we consider the parameter set a = 100, b = 72, c = 1, d = 0.9, e = 0.5 and D N = D P = 0.1.For this parameter set, the resulting patterns for the prey species with different values of M are encapsulated in Figure 12.From this figure, we can observe the transition from quasi-periodic to stationary patterns through both the periodic and quasi-periodic patterns.The nature of the dynamic patterns can be easily understood from the corresponding phase portraits presented in Figure 13.As the patterns for the predator species are qualitatively similar to that of the prey species, we did not include the corresponding results to present the study in a concise manner.Further increments in the value of M retain the stationary nature of the patterns.However in this case, sufficiently large M makes the predator distribution homogeneous in space.The spatial heterogeneity in predator distribution can be regained by adequately increasing the value of D P .An interesting pattern formation scenario has been observed for b = 60 and M = 3.5 along with the other parameter values mentioned above.The resulting pattern for the prey species and the corresponding phase portrait are encapsulated in Figure 14.In this case, we can observe from Figure 14(a) that the population distribution is stationary in a region inside the spatial domain and it becomes oscillatory towards the boundaries.Also, the limit cycle presented in Figure 14(b) is away from the spatially homogeneous steady state which we did not observe for other periodic patterns.
Finally, we explore the effects of nonlocal prey consumption on the spatiotemporal dynamics for the parameter set considered for Figure 2. In this case, nonlocality does not have any drastic effect on the emerging traveling wave solutions for b ∈ (20.7807, 20.9176) where large extent of nonlocality can only make the transition zones more oscillatory.However, sufficiently large extent of nonlocality can eventually transform the traveling wave solution to a stationary pattern for b ∈ (20.9176, 23.2295), and hence, it reassures the stabilizing effect of nonlocal prey consumption.For the sake of brevity, we did not incorporate the corresponding numerical results in this study.

Discussion
In real ecosystems, predators generally utilize various prey species as their potential food sources.Thus, it becomes really important to model prey-predator interactions by considering the predator species being generalist in nature to understand properly the functioning of a particular ecosystem.However, there is a scarcity of literature which dealt with the spatiotemporal dynamics of a generalist predator-prey model.Most of the studies in this direction assumed the predator being a specialist one which can not persist in the absence of its focal prey species.In this article, we have considered the prey-predator model proposed by Erbach et al. [16], and investigated how spatial dispersal and nonlocal prey consumption shape the population distributions under different circumstances.
To begin with, we have encapsulated some interesting temporal dynamics of the considered model in Section 2. Depending upon parameter values, the non-spatial system can admit a wide variety of bifurcation scenario such as super-/sub-critical Hopf, saddle-node and homoclinic bifurcations, and bistability phenomenon   where two of the coexistence equilibria act as attractors.In the case of bistability, the initial population concentrations completely determine the fate of the ecosystem, and hence, external interventions such as environmental conditions and harvesting can eventually reduce the population levels permanently.This phenomenon is known as ecological regime shift and it is almost impossible to revert it back [9,44,46].Thus, it becomes crucial to predict this phenomenon adequately from the point of view of species conservation [9,44,46].On the other hand, it can have important consequences in biological control of some harmful pest species by their generalist predators.Further, we have incorporated the spatial diffusion in the modeling framework and investigated its influence on the resulting population distributions in Section 3. First, we have derived the necessary and sufficient conditions for Turing instability by employing the linear analysis.Then, we have derived a cubic Stuart-Landau equation through weakly nonlinear analysis to provide quantitative information regarding the stationary patterns emerging near the Turing bifurcation threshold.Our numerical simulations show the existence of stationary, periodic spatially uniform and quasi-periodic patterns, and traveling wave solutions.However, we could not find the occurrence of any chaotic population distribution.Note that stationary patterns basically demonstrate isolated static high concentration patches separated by patches with relatively low concentration.Whereas, spatially heterogeneous dynamic patterns correspond to localized extinction and regeneration of a species.It is well established now that diffusive prey-predator models with specialist predators generally exhibit chaotic patterns in Hopf and Hopf-Turing domains (especially, with equal rates of diffusivity for both the species) [4,30,32].Thus, it seems contradictory to our results, however, it might be a consequence of our consideration of predators being generalist in nature.On the other hand, Chakraborty [13] showed the existence of chaotic patterns for a generalist predator-prey model where predation is modeled via Holling type II functional response.But, this type of functional response does not reflect the prey-switching phenomenon, and hence, it seems inappropriate for generalist predators [16].Nevertheless, our model takes care of the prey-switching phenomenon appropriately by assuming Holling type III functional response for predation.Thus, our study indicates that chaotic patterns are less probable to find in the case of generalist predators.
The emergence of traveling wave solutions connecting two coexistence steady states which contribute to the temporal bistability phenomenon can have significant ecological relevance.Firstly, it indicates that the incorporation of spatial dispersal eradicates the temporal bistability phenomenon and the system asymptotically approaches the uniform state with higher density levels.Thus, our study suggests that spatial dispersal may help a system to avoid any sort of regime shift and it can be useful from the point of view of species conservation.This observation also suggests that biological control of a harmful pest species through generalist predators may not be successful.Another remarkable fact about the traveling wave solutions is that its presence in the Turing domain.In his seminal work [52], Turing suggested that the stationary spatially heterogeneous patterns are expected to form in Turing domain.However, we should not associate traveling waves to a stationary spatially heterogeneous structure.Thus, it possesses a contradiction to Turing's idea and demonstrates a limitation of the linear analysis.This type of limitation has already been documented in literature.For example, Vanag and Epstein [53] demonstrated the emergence of oscillatory spatially heterogeneous patterns in Turing domain for a bistable reaction-diffusion system with periodic boundary conditions.Similarly, Aragón et al. [1] found oscillatory and even chaotic patterns in Turing domain by investigating a toy problem with no-flux boundary conditions.
Finally, we have further extended the spatial system by incorporating nonlocal prey consumption and explored its subsequent influence on the resulting population distributions in Section 4. It is well established that the consideration of nonlocal prey consumption can not produce stationary Turing patterns in Rosenzweig-MacArthur model [6].However, in the presence of intra-specific competition among predators [5] and here for generalist predators-it can induce stationary Turing patterns.In this case, we have derived the Turing instability threshold through linear analysis.Our extensive numerical simulations suggest that nonlocal prey consumption acts as a stabilizing force as sufficiently large extent of nonlocality can eventually transform a dynamic pattern to a stationary one.This observation is consistent with our previous work [32] on nonlocal intra-specific interactions.Another interesting outcome of nonlocal prey consumption is the pattern involving both the stationary and dynamic regimes (see Fig. 14(a)).The associated phase portrait of the spatially averaged densities of both the species reveals a stable limit cycle located away from the corresponding spatially uniform coexistence steady state (see Fig. 14(b)).Predominantly, limit cycles (both stable and unstable) arise around a coexistence state in population dynamics.However, this is not the case here because of the existence of another spatially heterogeneous coexistence steady state due to the combined effect of spatial dispersal and nonlocal prey consumption.
Though we have presented a thorough investigation regarding the spatiotemporal dynamics of the considered model in this study, we have left one interesting parameter set for which Erbach et al. [16] showed the bistability between a stable equilibrium and a stable limit cycle.Thus, it would be interesting to examine the resulting patterns corresponding to this bistability phenomenon.Also, top-hat kernel has been used as a starting point to model nonlocal prey consumption in this study.Nonetheless, characteristics of nonlocal interactions are species dependent, and accordingly, other types of kernels such as triangular, parabolic, Gaussian and Laplacian kernels are considered in literature [35,41,47].Shape of these kernels are different with respect to standard deviation and kurtosis [35].For example, a kernel with large standard deviation and low kurtosis basically leads to a shape with broad peak and thin tail [35].Generally, large amounts of information regarding a particular species are compressed to obtain the shape of a kernel function, and hence, it encapsulates a mere approximation of real events [51].As a result, patterns obtained for a specific kernel function may be inconvenient to observe in reality for a certain species.Therefore, it would be another interesting research prospect to investigate the influences of different kernels on the resulting population distributions for our model.Further, we would like to check the sensitivity of the patterns with respect to different boundary conditions as was done in [40].We will explore these interesting problems in near future.

Figure 1 .
Figure 1.Bifurcation diagram for the prey population density (N ) with respect to the parameter b for the non-spatial model (2.3)-(2.4).The other parameter values are a = 100, c = 1, d = 0.9 and e = 0.5.Blue color corresponds to the stable coexistence (both stationary and oscillatory) of both the species and red color corresponds to the unstable steady state.Two super-critical Hopf bifurcation thresholds (b 1 H = 48.116and b 2 H = 73.435)are indicated by black vertical lines.

Figure 2 .
Figure 2. Bifurcation diagram for the prey population density (N ) with respect to the parameter b for the non-spatial model (2.3)-(2.4).The other parameter values are a = 100, c = 1, d = 0.1 and e = 0.5.Blue solid curves correspond to the stable coexistence steady states and red solid curves denote either unstable coexistence steady state or unstable limit cycle.Further, red dashed curve corresponds to the saddle point.Here, sub-critical Hopf bifurcation, homoclinic bifurcation and two saddle-node bifurcation thresholds are given by b H = 20.9176,b HM = 21.014,b 1 SN = 20.7807 and b 2 SN = 23.2997,respectively.
(3.1)-(3.2) using the method of multiple-scales.For this purpose, let us denote the reaction parts of the system (3.1)-(3.2) by F (N, P ) and G(N, P ) such that F (N, P ) := N (1 − N ) − N 2 P 1+aN 2 and G(N, P ) := bN 2 P 1+aN 2 + cP 1+dP − eP .By employing the Taylor series expansions of F and G up to third order around the spatially homogeneous coexistence steady state E * = (N * , P * ), we can write the system (3.1)-(3.2) as

Figure 3 .
Figure 3. Plots of temporal-Hopf thresholds (black vertical lines) and Turing curve (red curve) in (b, D P )-parameter space for the spatial model (3.1)-(3.2) with the parameter values a = 100, c = 1, d = 0.9, e = 0.5 and D N = 0.1.The light blue regions correspond to a situation where the diffusion-driven instability is not possible.

Figure 4 .
Figure 4. Plots of stationary Turing pattern for the prey species (panel (a)), and spatial density distributions for both the species (panel (b)).The used parameter values are a = 100, b = 73.5, c = 1, d = 0.9, e = 0.5, D N = 0.1 and D P = 50.The pattern and spatial distribution are presented for a certain time period and for a particular time instant, respectively, after neglecting the initial transients.

Figure 5 .
Figure 5. Plots of different dynamic patterns ((a) quasi-periodic and (b) periodic) for the prey species with (a) D P = 0.1 and (b) D P = 10.The other parameter values are a = 100, b = 72, c = 1, d = 0.9, e = 0.5 and D N = 0.1.Both the patterns are presented for a certain time period after neglecting the initial transients.

Figure 6 .
Figure 6.Phase portraits of the spatially averaged densities for both the species corresponding to the dynamic patterns illustrated in Figure5.Here, the large red dots represent the corresponding spatially homogeneous coexistence steady states.

Figure 8 .
Figure 8. Plots of temporal-Hopf threshold (black vertical line) and Turing curves (blue and red curves) in (b, D P )-parameter space for the spatial model (3.1)-(3.2) with the parameter values a = 100, c = 1, d = 0.1, e = 0.5 and D N = 0.1.The light blue regions correspond to a situation where the diffusion-driven instability is not possible.

Figure 9
demonstrates that the system (3.1)-(3.2) admits traveling wave solution connecting the steady states E (1) * = (0.1783, 19.2607) and E (3) * = (0.6377, 23.671) for the parameter values mentioned in the caption.Note that this set of parameter values lies in the Hopf-Turing domain for the steady state E (1) * unstable.In this case, only the stationary Turing patterns are possible to emerge.

Figure 10 .
Figure 10.(a) Plot of spatial distribution of the densities of both the species after they settle in the stationary regime for D P = 20.(b) Plot of the temporal evolution of spatially averaged density of the prey species (< N >) for different values of D P .The used parameter values are a = 100, b = 22, c = 1, d = 0.1, e = 0.5 and D N = 0.1.

Figure 11 .
Figure 11.Plots of relative position of the Turing curves in (b, D P )-parameter space for different values of the extent of nonlocality (M ).Red, blue and magenta curves are drawn for M = 0, M = 1 and M = 2, respectively.The other parameter values are a = 100, c = 1, d = 0.9, e = 0.5 and D N = 0.1.Both the black vertical lines correspond to the temporal-Hopf thresholds and the light blue regions correspond to a situation where the diffusion-driven instability is not possible.

Figure 12 .
Figure 12.Plots of dynamic and stationary patterns for the prey species with (a) M = 1, (b) M = 2, (c) M = 3.5 and (d) M = 4.The other parameter values are a = 100, b = 72, c = 1, d = 0.9, e = 0.5 and D N = D P = 0.1.All the patterns are illustrated for a certain time period after neglecting the initial transients.

Figure 13 .
Figure 13.Phase portraits of the spatially averaged densities for both the species corresponding to the dynamic patterns illustrated in Figures 12(a), (b) and (c), respectively.Here, the large red dots represent the corresponding spatially homogeneous coexistence steady states.

Figure 14 .
Figure 14.Plots of a dynamic pattern for the prey species (panel (a)) and the corresponding phase portrait of the spatially averaged densities for both the species (panel (b)) with the parameter values a = 100, b = 60, c = 1, d = 0.9, e = 0.5, D N = D P = 0.1 and M = 3.5.The large red dot in the phase portrait represents the corresponding spatially homogeneous coexistence steady state.