Spatial dynamics of predator-prey system with hunting cooperation in predators and type I functional response

1.   Introduction

Spatial patterns formation in nature are ubiquitous, with illustrations like zebra stripe patterns on animals skin, Turing patterns in a coherent quantum field, or diffusive patterns in predator-prey models ^[1,2,3,4,5]. The spatial factors of species interplay has been recognized as a vital component in how ecological communities are created and ecological interplay occur over a broad limit of temporal and spatial scale ^[6]. Spatial population distribution is of major importance in the study of ecological systems ^[7,8,9]. Mechanisms and scenarios characterizing the spatial population distribution of ecological species in spatial habitat are a focus of special interest in population dynamics. The spatial population distribution is affected by the proliferation capacity of the species and interactions between individuals ^[10].

Spatial pattern formation of predator-prey systems have started based upon the elementary work of A. M. Turing on morphogenesis ^[11]. The spatial predator-prey systems are studied to comprehend the role of random mobility of the prey and predator, inside their residence. A fully comprehensive elucidation of the spatial impact on ecological species interplays can be observed in the book, written by Okubo et.al. ^[12].

Spatial mathematical model is an appropriate tool for investigating fundamental mechanism of complex spatiotemporal population dynamics. An appropriate mathematical structure to explain the spatial aspect of population dynamics is specified by reaction-diffusion equations. Reaction-diffusion models were initially applied to describe the ecological pattern formation by Segel and Jackson in 1972 ^[3], based on the primary work of Turing ^[11]. Over the last several decades, a lot of articles have been published on the spatial dynamics of predator-prey model based on reaction-diffusion equations and different types of patterns have emerged for these models ^{[2,3,4,5,12,13,14,15,16,17,18,19]}.

Cooperative behavior can stimulate a relation among the population density and per capita population growth rate ^[20,21]. Ecologists have accepted several mechanisms for stimulating cooperative behavior in prey, namely cooperating reproduction, foraging capacity, etc. The cooperative behavior in prey may be generated by predation or by procedure inborn to the prey lifespan history ^[21]. Theory has pervasively payed attention to cooperative behavior in preys ^{[22,23,24,25,26,27,28,29,30]} and cooperative behavior in predators is less studied and poorly understood ^[31,32,33], in particular when space is considered explicitly. A mathematical model of prey and predator population interplay with cooperative behavior in predators through the system of nonlinear ordinary differential equations has been studied in non-spatial domain by Alves et.al. ^[34]. Motivated from their work, we modify and extend the model in a spatial domain to study its spatial dynamics.

Mechanisms of spatial pattern dynamics of predator-prey systems with cooperative behavior in predators have been comparatively new and to the best of our knowledge, not studied so far. The objective of this current investigation is to create deep intuition into methods of spatial pattern formation in predator-prey model with cooperative behavior in predators. Here, we investigate how distinct intensity of cooperation rate, basic reproduction numbers of the predator and diffusion coefficients affects the spatial patterns of predator-prey interaction.

This paper is organized as follows: In section Ⅱ, we formulate a reaction-diffusion predator-prey system with zero flux boundary conditions and non-zero initial conditions. Furthermore, we analyze the non-spatial dynamics of the model with hunting cooperation in predators and give a survey of the linear stability analysis in section Ⅲ. In section Ⅳ, we obtain the sufficient condition for Turing bifurcation with zero-flux boundary condition. In section Ⅴ, we carry out a series of numerical simulations to reveal that there is a large variety of different spatiotemporal dynamics in the spatial model for different intensities of predator's hunting cooperation, basic reproduction number of the predator and the rates of diffusion coefficients. The paper ends with a discussion.

2.   Model description

By incorporating diffusion and Holling type I functional response in the general predator-prey system with cooperative behavior in predators ^[34], we obtain the following diffusive predator-prey model as

where $U(T)$ and $V(T)$ are the densities of prey and predator population at time $T$ and location $(x_{1}, x_{2})$, respectively. Here, $r$ is the intrinsic growth rate of the prey and $K$ is its carrying capacity. The parameter $\lambda$ ($\lambda > 0$) is the attack rate per predator and prey, a ($a > 0$) describes the predator cooperation in hunting ($aV$ is cooperation term). We consider a Holling type I functional response of the form

which depends on both prey and predator densities, thereby reflecting hunting cooperation (handling-driven). The parameter $e$ is conversion efficiency and $m$ is the per capita mortality rate of predators. The non-negative constants $D_{1}$ and $D_{2}$ are the diffusion coefficients for prey and predator densities respectively.

We now non-dimensionalize the system Eqs (1) and (2) by introducing the non-dimensional variables

and non-dimensional parameters

and obtain the modified model as

where the positive constant $D$ is the ratio of diffusion coefficients of prey and predator densities and $\nabla^2\big( = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\big)$ is the usual Laplacian operator in two dimensional space $\textbf{R} = (x, y)$. To make certain that spatial patterns are governed by reaction-diffusion equations, system equations (6) and (7) is to be analyzed with the following non-zero initial conditions

and zero-flux (Neumann) boundary conditions

where $L$ denotes the size of the system in the direction of $x$ and $y$. $\nu$ is outward unit normal on the boundary $\partial \Omega$. Condition (9) implies that no individual species leave the domain.

3.   Analysis of the non-spatial model

In absence of diffusion, the equilibrium points of the system Eqs (6) and (7) are given by

The system Eqs (6) and (7) has following ecologically significant equilibrium solutions: (ⅰ) $(0, \; 0)$ (both extinct), (ⅱ) $(\omega, \; 0)$ (predator extinct) and (ⅲ) $(X^{*}, \; Y^{*})$ (coexist), where $X^{*} = \frac{1}{(1+\alpha Y^{*})}$ and $Y^{*}$ is the positive solution of

The cubic Eq (12) have one and two positive real roots when ($\omega > 1$, $\sigma \alpha > 1$ or $\sigma \alpha < 1$) and ($\omega < 1$, $\sigma \alpha > 1$) respectively.

The variational matrix about the equilibrium point $(X^{*}, \; Y^{*})$ is given by

(ⅰ) At $(0, 0)$, the variational matrix is

whose eigenvalues are $-1$ and $\sigma(>0)$. Hence, the system is unstable at the origin.

(ⅱ) At $(\omega, 0)$, the variational matrix is

whose eigenvalues are $-\sigma$ and $\omega -1$. Hence, the system is stable if $\omega < 1$.

(ⅲ) At $(X^{*}, Y^{*})$, the variational matrix is

and, the corresponding characteristic equation is

where $B_{1} = Y^{*}(1+\alpha Y^{*})+\frac{2\sigma}{\omega(1+\alpha Y^{*})}-\frac{\alpha Y^{*}}{(1+\alpha Y^{*})}-\sigma, \; B_{2} = \frac{Y^{*}(1+\alpha Y^{*})\left(\alpha \sigma + (1+\alpha Y^{*})^{2}\right)\omega-2\alpha \sigma Y^{*}}{(1+\alpha Y^{*})^{2}\omega}$.

The equilibrium solution $(X^{*}, \; Y^{*})$ is locally asymptotically stable if and only if

$\left(\sigma+Y^{*}(-1+\alpha(1-Y^{*}(2+Y^{*}\alpha)+\sigma))\right)\omega-2\sigma < 0$

and

$Y^{*}(1+Y^{*}\alpha)\left((1+Y^{*}\alpha)^{2}+\alpha\sigma\right)\omega-2Y^{*}\alpha\sigma > 0$.

4.   Analysis of the spatial model

Interior equilibrium point $(X^{*}, \; Y^{*})$ of non-spatial system is spatially homogenous steady state, that is, constant in space and time for the reaction-diffusion system Eqs (6) and (7) (diffusive model). We assume that $(X^{*}, \; Y^{*})$ is stable in non-spatial system of (6–7) which means the spatially homogenous steady state is stable with respect to spatially homogenous perturbations. Though the diffusion is often considered as a stabilizing process, it is a well known fact that diffusion can make a spatially homogenous steady state linearly unstable (Turing instability) with respect to heterogenous perturbations in a system of two interacting species ^[3,11]. The condition for Turing instability may be obtained by introducing a small heterogenous perturbation of the homogenous steady state as follows:

where $\epsilon_{1}$ and $\epsilon_{2}$ are two non-zero real numbers and $k = \left(k_{x}, k_{y}\right)$ such that $k^{2} = \left(k_{x}^{2}+k_{y}^{2}\right)$ is the wave number.

Substituting (16–17) into (6–7) and then linearizing it about interior equilibrium point $(X^{*}, \; Y^{*})$, we obtain the variational matrix as

The corresponding characteristic equation is

where

$C_{1}(k^{2}) = \left(1+\delta \right)k^{2}+\frac{2\sigma}{\omega(1+\alpha Y^{*})}-\sigma-Y^{*}\Big(\frac{\alpha}{1+\alpha Y^{*}}-1-\alpha Y^{*}\Big)$,

and

$C_{2}(k^{2}) = \delta k^{4}+\Big(-\delta \omega \alpha Y^{*}+ \big(2\sigma+(1+\alpha Y^{*})(Y^{*}+\alpha Y^{*2}-\sigma)\omega\big)\Big)(1+\alpha Y^{*})k^{2}+\Big((1+\alpha Y^{*})^{2}+\sigma \alpha\Big)\omega Y^{*}(1+\alpha Y^{*})-2\sigma \alpha Y^{*}$.

By Routh-Hurwitz criterion, the system (6–7) will be stable about $(X^{*}, \; Y^{*})$ if $C_{1}(k^{2}) > 0$ and $C_{2}(k^{2}) > 0$. As the parameters, $\delta\; \Big( = \frac{d_{1}}{d_{2}}\Big)$ and $k^{2}$ are all positive and $\frac{2\sigma}{\omega(1+\alpha Y^{*})}-\sigma-Y^{*}\left(\frac{\alpha}{1+\alpha Y^{*}}-1-\alpha Y^{*}\right) > 0$ (by the stability of the non-spatial model of (6–7)), hence, $C_{1}(k^{2}) >$ always positive. Hence, the condition for diffusive instability is $C_{2}(k^2) < 0$.

The polynomial function $C_{2}(k^{2})$ has a minimum for some value of $k^{2}$, say $k_{min}^{2}$, where

Hence, the minimum value of $k$ for which Turing instability will occur, is, $C_{2}(k_{min}^{2}) < 0$. Therefore, substitute $k_{min}^{2}$ in $C_{2}(k^{2})$, we get the sufficient condition for Turing instability

The interval of the wave number for which Turing instability take place is $(k_{-}, k_{+})$ and in this interval, we have $C_{2}(k^{2}) < 0$, where

and

5.   Numerical simulations

We will now investigate the numerical results of spatiotemporal models, namely (6–7). For numerical simulation, we set $\sigma$ and $\omega$ as $\sigma = 10.0$, $\omega = 0.8$, and consider $\alpha$ and $D$, as controlling parameters. For these values of parameters, the positive equilibrium points are (0, 0), (0.8, 0), (0.6159, 1.4171) and (0.6431, 1.2614). The steady state (0.6159, 1.4171) is stable and (0.6431, 1.2614) is unstable. Throughout our study in the spatiotemporal domain, we have considered the stable steady state (0.6159, 1.4171). Please note that the non-dimensional parameter $\omega = \frac{e\lambda}{m}K$, comprising of the dimensional carrying capacity, attack rate, per capita mortality rate of predators and the conversion efficiency, can also be interpreted as the basic reproduction number of the predator, which is defined as the average number of offspring produced by a single predator during its life time, when introduced into the prey population at carrying capacity.

We now simulate the spatial model in two dimensional space with the help of finite difference scheme for spatial derivatives. The forward Euler's numerical method is used for the non-spatial part of model (6–7) and general finite difference scheme of five point is used for the spatial part. The reaction-diffusion partial differential equations, given by (6–7), is numerically solved by using splitting method ^[36,37]. The numerical values for the step sizes of time and space have been selected adequately small for avoiding the numerical artifacts. In this study, we have employed statistically uncorrelated Gaussian white noise perturbation in space, which is mathematically denoted in two dimensional case as

where $\epsilon_{ij}$ and $\eta_{ij}$ are statistically uncorrelated Gaussian white noise perturbations with zero mean and fixed variance in two dimensional space.

For spatial model, we perform all the numerical simulations of the system (6–7) over the non-zero initial condition and zero-flux boundary conditions, in two dimensional spatial domain. The domain size is $70 \times 70$ with time-step $\Delta t = 0.001$ and space-step $\Delta x = \Delta y = 1.0$. The parameter values of $\sigma$ remain same and $\omega$, $\alpha$ are used as the controlling parameter (just like the temporal case).

Note: The Neumann zero-flux conditions are placed at boundary of the numerical domain in two dimensional problems. The size of the domain is chosen large enough so that the impact of the boundaries has been kept as small as possible during the simulation time.

We now demonstrate diffusive induced instability (Turing instability) and the corresponding patterns formation for the system (6–7). Although, the sufficient conditions for Turing instability were obtained analytically in previous section, whether they are satisfied with our corresponding set of parameter values, is yet to be tested. In order to do so, we sketch the polynomial function $C_{2} (k^{2})$ for distinct values of $D_{1}$, $\omega$ and $\alpha$ (other parameter values are fixed, namely, $\sigma = 10$, $D_{2} = 0.07$). Figures 1, 2 and 3 shows the plot of $C_{2} (k^{2})$ against the wave number (k) for different values of $D_{1}$, $\omega$ and $\alpha$ respectively. We observe that the sufficient condition of the diffusive instability min $C_{2} (k^{2}) < 0$ holds, when $D_{1}$ is adequately large, beginning from $D_{1} = 2.2$ (see Figure 1). The polynomial function $C_{2} (k^{2}) < 0$ holds, the wave number (k) fit in the interval $(k_{-}, k_{+})$ and the length of the interval increases with an increase in the diffusion coefficient $D_{1}$. We next plot $C_{2} (k^{2})$ for different values of $\omega$ (other parameters remain fixed) and we observe that the conditions for the Turing instability holds for $\omega \leq 0.804$ approximately, and the length of the interval $(k_{-}, k_{+})$ decreases with an increase in $\omega$ (see Figure 2). We next plot $C_{2} (k^{2})$ for different values of $\alpha$ (other parameters remain fixed) and we observe that the conditions for the Turing instability holds for $\alpha \leq 0.45$ approximately, and the length of the interval $(k_{-}, k_{+})$ decreases with an increase in $\alpha$ (see Figure 3).

In Figure 4, we illustrate the density distributions of prey and predator which covers two types of pattern namely spots and stripes. Figure 5 and Figure 6, shows the two dimensional stationary diffusive patterns of the model (6–7) at time t = 500 with $\omega$ = 0.800 and diffusion coefficients $D_{1} = 3.1$ and $D_{2} = 0.07$ for the prey and predator population respectively. In these figures hexagonal patterns (spots) prevail over the entire habitat eventually. In Figures 4, 5 and 6, it is observe that blue spots (minimum density of X) are distributed on a reddish background (maximum density of X), that is, the preys are segregated with low population density.

6.   Discussions and conclusions

In theoretical ecology, intensive studies of the mechanisms and scenarios of pattern formation in models of interacting populations have always been an attraction as their perception help to enhance the understanding of real-world ecological systems. In this paper, we have considered a diffusive predator-prey model with hunting cooperation in predators and type I functional response, under non-zero initial conditions and zero-flux boundary conditions. We have provided the elaborate analysis of both temporal and spatiotemporal models and studied possible scenarios of pattern formation in the diffusive predator-prey model with hunting cooperation in predators. While studying the spatiotemporal model, we first obtain the condition for diffusive instability and identified the corresponding domain in the space of controlling parameters. The hunting cooperation coefficient $(\alpha)$, the basic reproduction number of the predator $(\omega)$ and the diffusion coefficient of the prey $(D_{1})$ are the controlling parameters in our study. Using the parameter values from Turing domain, we investigate the properties of the system using extensive numerical simulations.

By varying the values of cooperation coefficient, we get dissimilar types of diffusive patterns, namely, patchy pattern (spots), stripe pattern and mixed pattern (spot-stripe). From the point of view of population dynamics, one can observe that the spot formation for preys imply that the preys are scattered and isolated with low density and the remainder region is high density, which means that the preys may break out in the area and are safe. Similarly, spot formation in predators convey that with hunting cooperation, the predators are scattered and isolated but still survives. The methods and consequences in this study may amplify the systematic investigation of spatial pattern formation in the predator-prey systems, and may nicely enforce in some different research dimensions. Further analyze are important to study the patterns dynamics of some more diffusive ecological. It would be interesting to study the traveling waves in the spatial predator-prey models with hunting cooperation in predators with type Ⅱ, Ⅲ or type Ⅳ functional responses ^{[35,36,37,38]}. This article highlights a number of research areas for future consideration in spatial pattern formation.

Acknowledgments

Teekam Singh gratefully acknowledges the School of Allied Sciences, Department of Mathematics, Graphic Era (deemed to be) University and Graphic Era Hill University Dehradun-248 002, Uttarakhand, India.

Conflict of interest

The authors declare no conflict of interest.