Simulation of single-species bacterial-biofilm growth using the Glazier-Graner-Hogeweg model and the CompuCell3D modeling environment

1.   Introduction

Incorporation of spatial components into a prey-predator system can potentially lead to very rich and complex spatiotemporal dynamics than the corresponding non-spatial counterpart which is in accordance with experimental observations and theoretical studies. For instance, Gause ^[1] performed a laboratory experiment regarding the growth of paramecium and didinum and identified the significance of heterogeneous spatial distribution on the stabilization and long term survival of certain species. Working on the theoretical setup introduced by Turing ^[2], Levin and Segel ^[3] explained the patchy distribution of plankton ecosystem and Klausmeier ^[4] investigated vegetation patterns in semiarid region. Diffusive prey-predator systems are able to produce a wide class of spatiotemporal patterns such as stationary (spots, labyrinthine, and mixture of stripes and spots) and non-stationary (periodic, quasi-periodic and chaotic etc.) patterns. Another class of non-stationary patterns includes traveling wave, periodic traveling wave and wave of invasion ^[5]. Prey-predator model with specialist predator, prey dependent functional response and linear death rate for predator is unable to produce stationary patterns through Turing instability. On the other hand, this type of models are capable to produce spatiotemporal chaotic patterns for equal diffusivity of both the species ^[6]. However, models with predator dependent functional response, density dependent death rate for predator or generalist predator are able to exhibit stationary patterns generated through Turing instability ^[7,8,9]. In this direction, transition from stationary pattern to different types of non-stationary patterns is getting lots of attention recently from researchers as the real instances of non-stationary patterns in ecological community are growing significantly. For example, plankton population in the Naroch Lakes in Belarus and vole population in northern Fennoscandia exhibit chaotic patterns ^[10,11].

There have been numerous evidences of Allee effect in ecological communities due to several ecological mechanisms such as difficulty in finding mates, habitat alteration and cooperative defense etc. ^[12,13,14]. Basically, Allee effect represents a positive correlation between population density and individual fitness of a species ^[12,13,14]. It is classified as the strong when the growth rate of a species becomes negative below a certain threshold density and the weak when the growth rate stays positive and increasing at low population density ^[13,15,16]. Incorporation of Allee effect in the modeling approach for prey-predator interactions is actually able to produce very rich temporal dynamics. For example, the classical Rosenzweig-MacArthur model exhibits temporal extinction scenario through homoclinic bifurcation in presence of the strong Allee effect ^[17]. Existence of two and three limit cycles in a temporal Leslie-Gower type prey-predator model with additive Allee effect in prey growth was shown in ^[18,19]. In case of two limit cycles, sub-critical Hopf bifurcation leads to an unstable limit cycle which is surrounded by a stable limit cycle ^[18]. In case of three limit cycles, homoclinic bifurcation leads to one limit cycle and Hopf bifurcation gives rise to other two ^[19]. Recently, we have identified that the temporal extinction for both the populations occurs through heteroclinic bifurcation in a prey-predator system with additive Allee effect, monotonic functional response with saturation and density dependent death rate for predators in ^[20].

Understanding the spatiotemporal dynamics of prey-predator interactions in presence of Allee effect and the corresponding spatial pattern formation has been emerging as an intriguing research theme in ecology and attracting the attention of many researchers. Du and Shi ^[21] studied a spatially heterogenous reaction-diffusion prey-predator system and found that the prey population can become extinct or persist or even blow-up depending on the initial conditions, heterogeneity of environment and certain parametric restrictions while the predator population maintains its density about a constant level. Further, they identified that heterogeneous environment can lead to the existence of Allee effect in case of strong prey growth. Wang et al. ^[22] considered a diffusive prey-predator system with strong Allee effect in prey growth and studied bifurcations of the spatially uniform, heterogeneous periodic and non-constant equilibrium solutions. In order to estimate the parametric region for spatial pattern formation, they examined the non-existence of the non-constant positive equilibrium solutions. In ^[23], the authors investigated the local and global asymptotic stability of positive homogeneous steady state and derived the conditions for Hopf bifurcation for a prey-predator model with Holling type Ⅱ functional response and additive Allee effect in prey growth. An extension of this model by incorporating the density dependent death rate for predators was presented in ^[24]. The authors studied the dissipation and persistence property, stability of the homogeneous steady states, and derived the conditions for the Hopf bifurcation of the positive homogeneous steady state. Further, they investigated the existence and non-existence of positive non-constant solutions of the system and presented some stationary Turing patterns only for the case of weak Allee effect. They found that density dependent death rate for predators plays a key role in producing different types of stationary patterns through Turing instability. Recently, a wide class of stationary and dynamic patterns for this model was reported in ^[20] in presence of both weak and strong Allee effect. They argued that the half-saturation constant plays a prominent role for the formation of this wide variety of patterns and identified the Hopf bifurcation as a necessary component to induce spatiotemporal chaos. Apart from these, they successfully found patchy invasion scenario for this model. Rao and Kang ^[25] studied a prey-predator system with multiplicative Allee effect in prey growth and Michaelis-Menten type functional response, and identified the strength of the Allee effect as the key component in formation of different types of stationary and dynamic patterns.

Mathematical modeling of ecological interactions is an evolving area where exposed discrepancies lead to alteration or reconstruction of the models. One important observation is almost all biological processes are not instantaneous in nature and require time delays such as time involved in the process of fertilization stage to birth and gestation period etc ^[26,27]. Time delay is incorporated into a prey-predator system with the understanding that the change in a specific population density at a particular time is regulated by both past and present population densities. Generally, the time delay is included in a prey-predator model through two different assumptions. For example, one can consider time delay in the prey growth term such that the prey population in absence of the predator population follows the well-known Hutchinson's equation which is actually the delayed logistic equation ^[28,29]. Another way to incorporate delay in a prey-predator model is to consider delayed numerical response term for the predator population and this type of time delay is known as the gestation delay ^[29,30]. That is, prey consumption does not immediately reflect in predator population in terms of offsprings but a time lag is required in order to take care of the gestation period. Kuang ^[27] has mentioned that non-delayed models of prey-predator interactions act as mere approximations of the real situations and ignore the reality. Therefore, it is necessary to investigate the consequences of time delay on the spatiotemporal pattern formation.

For the past few years, a considerable amount of efforts has been dedicated to investigate the dynamics of delayed spatiotemporal system for prey-predator interactions. According to the best of our knowledge, Roy Choudhury ^[31,32] first studied the spatial structure in a prey-predator model with Volterra-type distributed delays in the inter-species interaction terms by using linear analysis and multiple-scales perturbation method. Hadeler and Ruan ^[33] investigated the combined effects of delay and diffusion for delayed reaction-diffusion equations. Sen et al. ^[34] introduced a generalized approach to investigate the delay-induced spatiotemporal instability and reported the formation of spirals through Hopf-Turing transition. In ^[35], the authors considered a realistic modeling approach to explore the complex spatiotemporal behavior of populations and accordingly studied the stability and Hopf bifurcation for a spatiotemporal model with logistic type growth, nonlocal delay effect and Dirichlet boundary condition. Tian and Zhang ^[36] studied the spatiotemporal patterns in a delayed plankton system and identified that delay-induced Hopf bifurcation can lead to irregular spatial patterns. Banerjee and Zhang ^[37] investigated how Turing and delay-induced Hopf bifurcations together influence the spatial patterns for a ratio-dependent prey-predator model. They identified that comparatively large values of time delay and the ratio of diffusivity play prominent roles in destabilizing the system's dynamics. In a subsequent study ^[38], they explored the role of time delay in producing spatiotemporal chaotic patterns for a prey-predator model with Holling type Ⅱ functional response and density dependent death rate for predators. Apart from the chaotic patterns, the existence of stationary, quasi-periodic patterns were also reported in this study. Song et al. ^[39] presented results regarding persistence, stability and delay-induced Hopf bifurcation in a delayed diffusive ratio-dependent prey-predator model. They also calculated the formulae for determining the direction and stability of the Hopf-bifurcating periodic solutions and showed the existence of spatially homogeneous and inhomogeneous periodic solutions through numerical illustrations. For other related works on delay-induced Hopf bifurcation and its stability and direction for reaction-diffusion systems of prey-predator interactions, interested readers are referred to the articles ^{[40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56]}. Recently, Jankovic et al. ^[57] showed that the delay-induced spatiotemporal chaos is possible for single species population models. They considered the spatiotemporal dynamics of a single species population using two different mathematical formulations such as delayed diffusive logistic equation and integro-differential equation incorporating a spatial convolution. Also, spatiotemporal dynamics of delayed three-species food-chain models were investigated in ^[58,59,60].

Being motivated by the works discussed above, we aim to investigate the inherent relationship between time delay and emerging spatial patterns for a delayed diffusive prey-predator system with additive Allee effect in prey growth, monotonic functional response with saturation and density dependent death rate for predators. Majority of the existing literature on the delayed spatiotemporal dynamics of prey-predator interactions have considered one-dimensional spatial components and this simplification takes the outcome away from the reality. In this paper, our goal is to provide both the analytical and numerical results for the delayed spatiotemporal model with two-dimensional spatial components. Also as far as our knowledge is concerned, work on the delayed spatiotemporal dynamics in presence of the Allee effect, which is quite commonly observed for several ecological communities, is rare in literature. These facts actually make our investigation more complicated and interesting, and also set apart our proposed works in this paper from any other existing work in literature. We aim to provide verification of the analytical results regarding the direction and stability of the Hopf-bifurcating periodic solutions through numerical simulations and construction of relevant bifurcation diagram. In this study, we are especially interested in the onset of delay-mediated spatiotemporal chaos. Also, we would like to explore when the spatial and corresponding non-spatial systems will produce similar qualitative dynamical behaviors in presence of the time delay. We will present a comparative study between the spatial and non-spatial systems.

The remaining part of this paper is organized in the following fashion. We present the model and its description in the next section. Section 3 deals with the well-posedness of the problem and discusses about the associated spatially homogeneous steady states. In Section 4, we present linear analysis and derive the analytical expression for the delay induced Hopf bifurcation threshold. The succeeding section discusses about the analytical results regarding the properties of the Hopf bifurcating periodic solutions. In order to support our theoretical findings and to investigate the effect of time delay on stationary patterns, we present several numerical illustrations in Section 6. Finally, Section 7 ends this paper with a brief discussion.

2.   The delayed spatiotemporal model

In this section, we extend the spatiotemporal prey-predator model presented in ^[20,24] with additive Allee effect in prey growth, prey-dependent monotonic functional response with saturation and density dependent death rate for predator by incorporating the gestation delay for predator population as follows:

subjected to the non-negative initial conditions:

and zero-flux boundary conditions:

Here $\Omega$ is a bounded spatial domain in $\mathbb{R}^{n}$ with smooth boundary $\partial\Omega$ and $\overline{\Omega} = \Omega\cup\partial\Omega$. The zero-flux boundary conditions indicate that no individual can move across the boundary with $\frac{\partial}{\partial\nu}$ being the outward drawn normal derivative on the boundary $\partial\Omega$. The prey and predator population densities at position ${\bf{x}}\in\Omega\subset\mathbb{R}^{n}$ and time instant $t$ are represented by $N({\bf{x}}, t)$ and $P({\bf{x}}, t)$, respectively. The parameters $r$ and $K$ stand for the intrinsic growth rate and the environmental carrying capacity in absence of the Allee effect, respectively, for the prey population. However, the effective environmental carrying capacity for prey population with additive Allee effect is actually different from $K$ ^{[61,62,63,64]}. The expression of effective carrying capacity for prey population is explicitly given below. The rate of severity of the Allee effect and the prey population size at which fitness is half of its maximum value are respectively represented by the parameters $m$ and $b$. Predator captures prey at a rate $c$ while $a$ denotes the half saturation constant. The parameter $s$ represents feed concentration for predator with $q$ and $\delta$ as natural and density dependent death rate for predator, respectively. Therefore, the quantities $sq$ and $s\delta$ represent respectively the effective natural and density dependent death rates for predator. All these parameters are positive constants. Further, $\tau$ represents the gestation delay for predator population which reflects the fact that the consumption of prey can not lead to instantaneous reproduction of predator. The fractional terms $\frac{m}{N+b}$ and $\frac{cN}{N+a}$ in the model (2.1)-(2.2) represent additive Allee effect in prey growth and prey density dependent functional response for predator, respectively. The constant diffusion coefficients of prey and predator populations are respectively denoted by $d_{1}$ and $d_{2}$, and $\Delta$ represents the Laplacian operator.

Depending on the parameter values the Allee effect can act as weak or strong in nature, which has been thoroughly investigated in ^[24]. Here, we reiterate these parameter restrictions for weak and strong Allee effects briefly as follows:

(a) If the parameters obey the inequality $\frac{b^{2}r}{K} < m < br$, then the Allee effect is weak. In this case without any predator, prey population can survive for a very tiny positive initial prey population size and eventually approaches to $N^{1} = \frac{r(K-b)+\sqrt{r^{2}(K-b)^{2}+4rK(br-m)}}{2r}$. Also, $N^{1}$ acts as the effective environmental carrying capacity in presence of the additive Allee effect for prey population ^{[61,62,63,64]}.

(b) If the parameters obey the inequality $\frac{b^{2}r}{K} < br < m < \frac{r^{2}(K-b)^{2}+4br^{2}K}{4rK}$, then the Allee effect is strong. In this case without any predator, the prey population persists and approaches to $N^{1}$ for any positive initial prey population size greater than a critical prey population threshold $N^{\text{crit}}$ and goes to extinction for any positive initial prey population size less than $N^{\text{crit}}$, where $N^{\text{crit}} = \frac{r(K-b)-\sqrt{r^{2}(K-b)^{2}+4rK(br-m)}}{2r}$. This $N^{\text{crit}}$ acts as the Allee effect threshold for strong Allee effect.

In the following sections, we will examine the system (2.1)-(2.4) both analytically and numerically. We would like to mention here that numerical simulations will be carried out in two-dimensional space though all analytical results hold for $n$-dimensional space.

3.   Well-posedness and steady states

In this study, we are dealing with prey and predator populations and for this reason it is pertinent to establish the existence, uniqueness, non-negativity and boundedness of the solutions of the concerned system (2.1)-(2.4). In this section, we discuss these issues along with the existence of possible spatially homogeneous steady states. Before proceeding further, let us introduce some useful notations. Let $\mathcal{X} = C(\overline{\Omega}, \mathbb{R}^{2})$ be a Banach space of continuous functions from $\overline{\Omega}$ to $\mathbb{R}^{2}$ and $\mathcal{C} = C([-\tau, 0], \mathcal{X})$ denote the Banach space of continuous functions from $[-\tau, 0]$ to $\mathcal{X}$ equipped with the usual supremum norm. For the sake of convenience, an element $\phi\in\mathcal{C}$ is a function from $\overline{\Omega}\times[-\tau, 0]$ into $\mathbb{R}^{2}$ and is defined by $\phi({\bf{x}}, \theta) = \phi(\theta)({\bf{x}})$. We define $\omega_{t}\in\mathcal{C}$ by $\omega_{t}(\theta) = \omega(t+\theta)$, $\theta\in[-\tau, 0]$ for any continuous function $\omega:[-\tau, T)\rightarrow\mathcal{X}$ where $T > 0$.

Theorem 3.1. For any given $\phi\in\mathcal{C}$ satisfying the initial conditions (2.3), the delayed spatiotemporal system (2.1)-(2.4) admits a unique solution defined on $[0, +\infty)$ and this solution remains non-negative and bounded for all $t\geq 0$.

Proof. For any $\phi = (\phi_{1}, \phi_{2})^{T}\in\mathcal{C}$ and ${\bf{x}}\in\overline{\Omega}$, we define $F = (F_{1}, F_{2}):\mathcal{C}\rightarrow\mathcal{X}$ by

Then we can rewrite the system (2.1)-(2.4) in terms of the following abstract functional differential equation:

where $U = (N, P)^{T}$, $\phi = (\phi_{1}, \phi_{2})^{T}$ and $AU = (d_{1}\Delta N, d_{2}\Delta P)^{T}$. One can easily observe that $F$ is locally Lipschitz in $\mathcal{X}$. It follows from ^{[65,66,67,68,69]} that the system (3.1)-(3.2) admits a unique local solution on the interval $[0, T_{max})$ where $T_{max}$ represents the maximal existential time for solution of the system (3.1)-(3.2). Also, since $\bf{0}$ = (0, 0) is a lower solution of the system (2.1)-(2.2), we can easily obtain that $N({\bf{x}}, t)\geq 0$ and $P({\bf{x}}, t)\geq 0$.

Now, it remains to prove the boundedness of the solutions of the system (2.1)-(2.4). From the equation (2.1), we can easily obtain that

Let $\widehat{N}(t)$ be a solution to the following ordinary differential equation:

Thus, we have $\widehat{N}(t)\leq \max\left\{K, \max_{x\in\overline{\Omega}}\phi_{1}({\bf{x}}, 0)\right\}$ for all $t\in[0, T_{max})$ ^[24,70]. Then by using the comparison principle for parabolic partial differential equations ^[71], we obtain that $N({\bf{x}}, t)\leq \widehat{N}(t)$. Therefore,

Now, from the equation (2.2) and using the bound for $N$, we deduce that $P({\bf{x}}, t)$ satisfies the following system:

where $\rho = \max\left\{K, \max_{x\in\overline{\Omega}}\phi_{1}({\bf{x}}, 0)\right\}$. Further, let $\widehat{P}(t)$ be a solution to the following delay differential equation:

In order to find a bound for $\widehat{P}(t)$, we recall a related result presented in ^[72,73]. This result states that for an equation of the form $\frac{du(t)}{dt} = Au(t-\tau)-Bu(t)-Cu^{2}(t)$ with $A, B, C > 0$ and $u(\theta) > 0$ for $\theta\in[-\tau, 0]$, we get $\lim_{t\rightarrow +\infty}u(t) = \frac{A-B}{C}$ when $A > B$ and $\lim_{t\rightarrow +\infty}u(t) = 0$ when $A < B$. Now, we assume that $M = \max\left\{0, \frac{(c-q)\rho-aq}{\delta(\rho+a)}\right\}$. Hence, using this result we obtain $\widehat{P}(t)\leq \max\left\{M, \max_{({\bf{x}}, \theta)\in\overline{\Omega}\times[-\tau, 0]}\phi_{2}({\bf{x}}, \theta)\right\}$ for all $t\in[0, T_{max})$. Then again using the comparison principle for parabolic partial differential equations ^[71], we obtain that $P({\bf{x}}, t)\leq \widehat{P}(t)$. Therefore,

From the above discussion, we obtain that $N({\bf{x}}, t)$ and $P({\bf{x}}, t)$ are bounded on $\bar{\Omega}\times[0, T_{max})$ and also we have $T_{max} = +\infty$ from the standard theory of semi-linear parabolic partial differential equations ^[74]. This completes the proof.

Now, we discuss about the spatially homogeneous steady states of the model (2.1)-(2.2). Any steady state of the model (2.1)-(2.2) satisfies the following system of coupled partial differential equations:

But in this study, we are only concerned with the spatially homogeneous steady states. Then, any non-negative solution of the following system of coupled algebraic equations corresponds to a spatially homogeneous steady state $E = (\widehat{N}, \widehat{P})$ of the model (2.1)-(2.2):

From the above two algebraic equations, one can easily notice that there always exists a trivial extinction steady state $E_{0} = (0, 0)$. For weak Allee effect, there exists a unique prey-only axial steady state $E_{1} = \left(\frac{r(K-b)+\sqrt{r^{2}(K-b)^{2}+4rK(br-m)}}{2r}, 0\right)\equiv\left(N^{1}, 0\right)$. But in case of strong Allee effect, apart from the steady state $E_{1}$ there exists another additional prey-only axial steady state $E_{2} = \left(\frac{r(K-b)-\sqrt{r^{2}(K-b)^{2}+4rK(br-m)}}{2r}, 0\right)\equiv\left(N^{\text{crit}}, 0\right)$, whose prey component actually indicates the threshold value for prey population below which prey species can not persist in the non-spatial counterpart of the considered model (2.1)-(2.2). In order to find a spatially homogeneous coexisting steady state $E_{*} = (N_{*}, P_{*})$, we need to solve a quartic equation in $N_{*}$ as follows:

where $A_{4} = \frac{r\delta}{K}$, $A_{3} = \frac{(2a+b)r\delta}{K}-r\delta$, $A_{2} = \frac{a(a+2b)r\delta}{K}+m\delta+c(c-q)-(2a+b)r\delta$, $A_{1} = \frac{a^{2}br\delta}{K}+2am\delta+bc(c-q)-a(a+2b)r\delta-acq$ and $A_{0} = a^{2}m\delta-a^{2}br\delta-abcq$. The corresponding component for predator population, i.e. $P_{*}$, is evaluated through $P_{*} = \frac{1}{\delta}\left(\frac{cN_{*}}{N_{*}+a}-q\right)$. In this study, we use numerical computation for obtaining coexisting steady state since it is very difficult to compute coexisting steady state analytically from a quartic equation. Also, note that depending on the parameter values we can obtain one or more than one spatially homogeneous coexisting steady states for the model (2.1)-(2.2). For detailed discussion on steady states one can go through the works presented in ^[24]. In this study, we concentrate only on the case when there exists a unique spatially homogeneous coexisting steady state for the model (2.1)-(2.2).

4.   Linear stability analysis and Hopf bifurcation

In this section, we investigate the delay induced Hopf bifurcation in the delayed spatiotemporal system (2.1)-(2.4) around a typical spatially homogeneous coexisting steady state $E_{*} = (N_{*}, P_{*})$ through the linearization technique. Before proceeding further, let us introduce some useful notations for the sake of convenience. Let us consider $0 = \lambda_{0} < \lambda_{1} < \lambda_{2} < \cdots < \lambda_{n} < \cdots$ be the eigenvalues of the operator $-\Delta$ on $\Omega$ under the zero-flux boundary conditions and let $\mathcal{F}(\lambda_{i})$ be the space of eigenfunctions corresponding to $\lambda_{i}$ in $C^{1}(\Omega)$. Let $\{\psi_{ij}:\; j = 1, 2, \dots, \dim\mathcal{F}(\lambda_{i})\}$ be an orthonormal basis of $\mathcal{F}(\lambda_{i})$, $\mathbb{X}_{ij} = \{\kappa\psi_{ij}:\; \kappa\in\mathbb{R}^{2}\}$ and $\mathbb{X} = \{{\bf{U}}\in [C^{1}(\overline{\Omega})]^{2}\mid\frac{\partial\bf{U}}{\partial\nu} = {\bf{0}}\; \text{on}\; \partial\Omega\}$. Then we can write

Linearizing the model (2.1)-(2.2) about the spatially homogeneous coexisting steady state $E_{*} = (N_{*}, P_{*})$, we obtain

where $\mathcal{D} = \text{diag}(d_{1}, d_{2})$, ${\bf{U}}({\bf{x}}, t)\equiv(U_{1}({\bf{x}}, t), U_{2}({\bf{x}}, t))$, ${\bf{U}}({\bf{x}}, t-\tau)\equiv(U_{1}({\bf{x}}, t-\tau), U_{2}({\bf{x}}, t-\tau))$ and

Let us define $\mathcal{L}{\bf{U}} = \mathcal{D}\Delta{\bf{U}}({\bf{x}}, t)+\mathfrak{J}_{1}{\bf{U}}({\bf{x}}, t)+ \mathfrak{J}_{2}{\bf{U}}({\bf{x}}, t-\tau)$. Note that $\mathbb{X}_{i}$ is invariant under the operator $\mathcal{L}$ for each $i\geq 0$. Further, $\mu$ is an eigenvalue of the operator $\mathcal{L}$ if and only if $\mu$ satisfies the equation $\text{det}(\mathfrak{J}_{i}-\mu\mathbb{I}_{2})\equiv\text{det}(-\lambda_{i}\mathcal{D}+\mathfrak{J}_{1}+ \mathfrak{J}_{2}e^{-\mu\tau}-\mu\mathbb{I}_{2}) = 0$ for some $i\geq 0$ and there exists an eigenvector in $\mathbb{X}_{i}$ for this case. Since our model is subjected to the zero-flux boundary conditions, a typical eigenfunction of (4.1) is of the form

where ${\bf{x}} = (x_{1}, x_{2}, \cdots, x_{n})\in\mathbb{R}^{n}$ and ${\bf{k}} = (k_{1}, k_{2}, \cdots, k_{n})$ represents the $n$-dimensional wave-vector. Then, the characteristic equation is given by

where $C_{1} = -(A_{11}+A_{22})+k^{2}(d_{1}+d_{2})$, $C_{0} = (A_{11}-k^{2}d_{1})(A_{22}-k^{2}d_{2})$, $D_{1} = -B_{22}$, $D_{0} = (A_{11}B_{22}-A_{12}B_{21})-k^{2}d_{1}B_{22}$ and $k^{2} = \bf{k}\cdot\bf{k}$.

In absence of the time delay (that is, when $\tau = 0$), the characteristic equation reduces to the following quadratic equation:

Working on the equation (4.3), one can obtain the conditions for Turing and spatiotemporal Hopf bifurcation for the corresponding non-delayed model which have been investigated earlier in ^[20,24]. Since the main goal of this study is to understand the destabilizing role of discrete time delay, we assume that the coexisting homogeneous steady state $E_{*} = (N_{*}, P_{*})$ is locally asymptotically stable for the non-delayed model. That is, we assume the following conditions hold: $C_{1}+D_{1} > 0$ and $C_{0}+D_{0} > 0$.

Now, we investigate the effect of time delay and deduce the conditions under which the model (2.1)-(2.2) undergoes delay induced Hopf bifurcation. In order to find the delay induced Hopf bifurcation threshold, we follow the procedure presented in ^[75]. Substituting $\mu = i\omega$ into the characteristic equation (4.2) and separating the real and imaginary parts, we arrive at the following two equations:

Eliminating the trigonometric functions from the above two equations, we arrive at the following quartic equation in $\omega$:

Now, note that one can visualize the equation (4.6) as a quadratic equation of $\omega^{2}$ and if both the values of $\omega^{2}$ are negative or imaginary, the time delay can not able to alter the stability property of the homogeneous coexisting steady state. But if we have at least one positive real value of $\omega^{2}$, then destabilization of the homogeneous coexisting steady state is possible through $\tau$. In order to find the Hopf bifurcation threshold in terms of the delay parameter $\tau$, we assume $\omega^{*}$ denotes a positive root of the equation (4.6). The simplest possible assumptions that would give us a unique positive root $\omega^{*}$ of the equation (4.6) is $C_{1}^{2}-2C_{0}-D_{1}^{2} > 0$ and $C_{0}^{2}-D_{0}^{2} < 0$. The expression for the unique positive root $\omega^{*}$ is given by

Then eliminating $\sin \omega\tau$ from the equations (4.4)-(4.5) and solving for $\tau$, we obtain

for $n = 0, 1, 2, 3, \cdots$. Now, setting

and by using the Butler's lemma ^[76], we obtain that the homogeneous coexisting steady state $E_{*} = (N_{*}, P_{*})$ is unstable for $\tau > \tau^{*}$. On the other hand, the characteristic equation (4.2) does not have any root on the imaginary axis for $\tau\in[0, \tau^{*})$. In conclusion, we obtain that $E_{*} = (N_{*}, P_{*})$ is locally asymptotically stable for $\tau\in[0, \tau^{*})$ and becomes unstable for $\tau > \tau^{*}$. But in order to ensure that the occurrence of instability is due to delay driven Hopf bifurcation, we need to verify the associated transversality condition:

Differentiating the both sides of the characteristic equation (4.2) with respect to $\tau$, we get

which eventually leads to

Thus, we obtain

Now, from equations (4.4)-(4.5) we have

Hence, equation (4.10) leads to

Thus, from the assumptions for the occurrence of the unique positive root $\omega^{*}$ of the equation (4.6), we obtain $\text{Re}\left(\frac{d\mu}{d\tau}\, \middle\vert\, _{\tau = \tau^{*}}\right)^{-1} > 0$. Further, we know that $\text{sign}\left[\frac{d}{d\tau}(\text{Re}(\mu(\tau)))\, \middle\vert\, _{\tau = \tau^{*}}\right] = \text{sign}\left[\text{Re}\left(\frac{d\mu}{d\tau}\, \middle\vert\, _{\tau = \tau^{*}}\right)^{-1}\right]$ and therefore, we conclude that the inequality (4.9) is satisfied.

Therefore, summarizing the above detailed analysis we arrive at the following theorem:

Theorem 4.1. The delayed spatiotemporal system (2.1)-(2.4) undergoes Hopf bifurcation around the spatially homogeneous coexisting steady state $E_{*} = (N_{*}, P_{*})$ at $\tau = \tau^{*}$ when $C_{1}+D_{1} > 0$, $C_{0}+D_{0} > 0$ and the assumptions for the unique positive $\omega^{*}$ hold.

Remark 4.1. From the expression (4.8) for $\tau^{*}$, we can notice that $\tau^{*}$ depends on the value of $k$ for the considered spatiotemporal system (2.1)-(2.4). In other words, for the considered system (2.1)-(2.4), there will exist a non-negative $k$ for which we will have our theoretical Hopf bifurcation threshold $\tau^{*}$. Also, it is important to note that for different values of system parameters the system (2.1)-(2.4) will need different non-negative $k$-values in order to achieve the Hopf bifurcation threshold. This fact will be demonstrated more clearly in Section 6. Here, we would like to mention that by considering $k = 0$ we can eventually achieve the Hopf bifurcation threshold for the corresponding non-spatial counterpart of the spatial model (2.1)-(2.2). We denote this threshold by $\tau^{*}_{0}$, where the subscript \enquote{$0$} signifies the fact $k = 0$. Hence, we conclude that the corresponding non-spatial model undergoes the Hopf bifurcation at the coexisting steady state $E_{*} = (N_{*}, P_{*})$ at $\tau = \tau^{*}_{0}$ provided the conditions and assumptions given in Theorem 4.1 are true when $k = 0$. Also, as per our knowledge there is no literature where the dynamics of the non-spatial counterpart of the model (2.1)-(2.2) have been investigated. Therefore, our findings in this section also provide key insights about the dynamics of the corresponding non-spatial model.

Now we briefly corroborate the conclusions made about the dynamics of the non-spatial model corresponding to the model (2.1)-(2.2) in the above remark through numerical illustration. For numerical illustration, we restrict ourselves to the case of strong Allee effect. The three figures presented in the Figure 1 have been obtained by varying the value of the gestation delay parameter $\tau$ in an increasing manner while all other system parameters were kept fixed (all the parameter values are mentioned in the caption of the Figure 1). Figure 1(a) represents a bistable scenario for a lower value of $\tau$ (for example, $\tau = 1.0$), where the total extinction steady state $E_{0} = (0, 0)$ is a stable node and the coexisting steady state $E_{*} = (3.665, 3.015)$ is a stable focus. In this case, the basins of attraction of these two steady states are kept apart by the stable manifold of the prey-only steady state $E_{2} = (0.559, 0)$ whose prey component acts as a threshold prey population density characterizing the strong Allee effect. For the considered parameter values, we have computed the Hopf bifurcation threshold as $\tau^{*}_{0} = 1.275$. Also, we have obtained $\frac{d\mu}{d\tau} = 0.0321-0.0674i$ at $\tau = \tau^{*}_{0}$ which eventually confirms the occurrence of the Hopf bifurcation. Therefore, we can expect the existence of Hopf-bifurcating limit cycles by setting the value of $\tau$ to be lower or greater than $1.275$ and Figure 1(b) shows the existence of the Hopf-bifurcating limit cycles for $\tau > 1.275$. From this figure we can observe that the size of the limit cycle is gradually increasing with the increment in the value of $\tau$. Here, we would like to mention that the emerging limit cycles around the coexisting steady state $E_{*} = (3.665, 3.015)$ is stable in nature. Therefore, the Hopf bifurcation can be classified as the supercritical Hopf bifurcation. With this increasing limit cycles another interesting feature such as the heteroclinic orbit emerges and this orbit connects the two prey-only steady states $E_{2} = (0.559, 0)$ and $E_{1} = (8.941, 0)$ when $\tau = 2.0623$. We have estimated numerically the threshold value of $\tau$ for the heteroclinic bifurcation to be approximately $\tau_{HO} = 2.0623$ (see the green colored orbit presented in Figure 1(b)). Through this bifurcation the basin of attraction for the stable limit cycles around the coexisting steady state $E_{*} = (3.665, 3.015)$ is demolished and $E_{0} = (0, 0)$ becomes the sole attractor. This scenario can be easily noticeable from the Figure 1(c). These scenarios are presented in a summarized fashion in Figure 2 by a one parameter ($\tau$) bifurcation diagram. In this figure, the blue and black vertical straight lines represent the Hopf and heteroclinic bifurcation thresholds, respectively. Two completely red horizontal straight lines respectively stand for the prey-only steady states $E_{1}$ and $E_{2}$ which are always saddle points and hence unstable. The completely green horizontal straight line along the $\tau$-axis stands for the extinction steady state $E_{0} = (0, 0)$ which is always a stable equilibrium point. The partially green and partially red straight line in the middle denotes the coexisting steady state $E_{*} = (3.665, 3.015)$ which is a stable focus for $\tau$ less than $\tau^{*}_{0}$ and an unstable focus for $\tau$ greater than $\tau^{*}_{0}$. The green curve, which lies between the blue and black vertical lines, represents the existence of stable limit cycle that arises through Hopf bifurcation and disappears through heteroclinic bifurcation. Note that the green and red colors signify stability and instability properties of the steady states and limit cycles, respectively. In the next section, we derive the analytical conditions for the stability of the Hopf-bifurcating limit cycle and its direction for the system (2.1)-(2.4).

5.   Direction and stability of the Hopf bifurcation

In this section, we investigate direction and stability of the Hopf bifurcating limit cycles by using center manifold reduction and normal form theory for partial functional differential equation ^[69,75,77].

Let us consider $\widetilde{N}(\cdot, t) = N(\cdot, t)-N_{*}$, $\widetilde{P}(\cdot, t) = P(\cdot, t)-P_{*}$, $\widetilde{U}(t)\equiv(\widetilde{N}(\cdot, t), \widetilde{P}(\cdot, t))$ and $\tau = \tau^{*}+\varsigma$ for $\varsigma\in\mathbb{R}$. Then, $\varsigma = 0$ corresponds to the Hopf bifurcation threshold of the system (2.1)-(2.4). In order to normalize the delay we rescale the time $t$ by $\frac{t}{\tau}$ and also recall the Banach space decomposition $\mathbb{X} = \bigoplus_{i = 1}^{\infty}\mathbb{X}_{i}$. Now, applying these new variables in the system (2.1)-(2.2) and then dropping the tildes for the simplicity of notation, we obtain the following transformed system in space ${\bf{C}} = C([-1, 0], \mathbb{R}^{2})$:

where for $\varphi = (\varphi_{1}, \varphi_{2})^{T}\in\bf{C}$, $L_{\varsigma}:\bf{C}\rightarrow\mathbb{R}^{2}$ and $f:\bf{C}\times\mathbb{R}\rightarrow\mathbb{R}^{2}$ are respectively given by

where

and

Now by the Riesz representation theorem ^[75], there exists a $2\times2$ matrix function $\eta(\theta, \varsigma)$ for $-1\leq\theta\leq0$, whose elements are of bounded variation such that

As a matter of fact, one can choose

where $\delta$ represents a Dirac delta function. For $\varphi\in C^{1}([-1, 0], \mathbb{R}^{2})$, let us define

and

Thus, the system (5.1) is equivalent to the following system:

where $U_{t}(\theta) = U(t+\theta)$ for $\theta\in[-1, 0]$. Now for $\psi\in C^{1}([0, 1], (\mathbb{R}^{2})^{*})$, we define

and a bilinear inner product

where $\eta(\theta) = \eta(\theta, 0)$. Then, we have $\mathcal{A}(0)$ and $\mathcal{A}^{*}$ as adjoint operators. From the discussion in the previous section, we notice that $\pm i\tau^{*}\omega^{*}$ are eigenvalues of $\mathcal{A}(0)$ and hence, they are also eigenvalues of $\mathcal{A}^{*}$.

Let us assume that $p(\theta) = (p_{1}, p_{2})^{T}e^{i\tau^{*}\omega^{*}\theta}$ be the eigenvector of $\mathcal{A}(0)$ corresponding to the eigenvalue $i\tau^{*}\omega^{*}$. Then, we have $\mathcal{A}(0)p(\theta) = i\tau^{*}\omega^{*}p(\theta)$. Now, from the definition of $\mathcal{A}(0)$ we obtain

Solving the above matrix equation, we have

Similarly, let us assume that $p^{*}(\gamma) = D(p_{1}^{*}, p_{2}^{*})e^{i\tau^{*}\omega^{*}\gamma}$ be the eigenvector of $\mathcal{A}^{*}$ corresponding to the eigenvalue $-i\tau^{*}\omega^{*}$. Then, we have $\mathcal{A}^{*}p^{*}(\gamma) = -i\tau^{*}\omega^{*}p^{*}(\gamma)$. From the definition of $\mathcal{A}^{*}$, we obtain

Therefore, from the above matrix equation, we obtain

In order to ensure $\langle p^{*}(\gamma), p(\theta)\rangle = 1$, we need to determine $D$. From the equation (5.8), we deduce that

Therefore, we can take $D$ as

such that $\langle p^{*}(\gamma), p(\theta)\rangle = 1$ and $\langle p^{*}(\gamma), \overline{p}(\theta)\rangle = 0$.

For the rest of this section, we will closely follow the notations as described in ^[75]. Now, we compute the coordinates to describe the center manifold $\bf{C}_{0}$ at $\varsigma = 0$. For this purpose, we define

On the center manifold $\bf{C}_{0}$, we have

where

and $z$ and $\overline{z}$ represent local coordinates for the center manifold $\bf{C}_{0}$ in the direction of $p^{*}$ and $\overline{p}^{*}$, respectively. We are only concerned with the real solutions and note that $W$ is real if $U_{t}$ is real. For the solution $U_{t}\in\bf{C}_{0}$ of the equation (5.6), since $\varsigma = 0$ and using (5.11), we obtain

Now, we rewrite the above equation (5.13) as follows:

where

Also from the equations (5.11) and (5.12), we get

Then, it follows that

where

Now, comparing the coefficients of the equations (5.15) and (5.17), we obtain

We need to compute the values of $W_{11}(\theta)$ and $W_{20}(\theta)$ as the expression for $g_{21}$ involves both $W_{11}(\theta)$ and $W_{20}(\theta)$ for $\theta\in[-1, 0]$. It follows from the equations (5.6) and (5.11) that

where

On the other hand, from the equation (5.12) on $\bf{C}_{0}$, we obtain

Also, using the expressions (5.12) and (5.20) in the equation (5.19), we obtain

Now, comparing the coefficients of $z^{2}$ and $z\overline{z}$ from the equations (5.21) and (5.22), we obtain

For $\theta\in[-1, 0)$, it follows from the equations (5.15) and (5.19) that

Again comparing the coefficients of $z^{2}$ and $z\overline{z}$ from the equations (5.20) and (5.24), for $\theta\in[-1, 0)$ we obtain

From the equations (5.23) and (5.25) and using the definition of $\mathcal{A}(\theta)$, we deduce that

The relation $p(\theta) = (p_{1}, p_{2})^{T}e^{i\tau^{*}\omega^{*}\theta} = p(0)e^{i\tau^{*}\omega^{*}\theta}$ leads the equation (5.27) to the following first order linear ordinary differential equation in $W_{20}(\theta)$:

Solving the above differential equation (5.28) for $W_{20}(\theta)$, we obtain

where $V_{1} = (V_{1}^{(1)}, V_{1}^{(2)})\in\mathbb{R}^{2}$ denotes a two-dimensional constant vector. In a similar manner, from the equations (5.23) and (5.26) and using the definition of $\mathcal{A}(\theta)$, we deduce that

where $V_{2} = (V_{2}^{(1)}, V_{2}^{(2)})\in\mathbb{R}^{2}$ also denotes a two-dimensional constant vector. Next, we need to find the appropriate values of the constant vectors $V_{1}$ and $V_{2}$. From the definition of $\mathcal{A}(0)$ and (5.23), we derive

and

where $\eta(\theta) = \eta(\theta, 0)$. Taking $\theta = 0$, we obtain from the equation (5.19) that

Then, we have

Comparing the coefficients of $z^{2}$ and $z\overline{z}$ from the both sides of the above equation (5.34), we obtain

where $(\rho_{1}, \rho_{2})^{T}$ and $(\rho_{3}, \rho_{4})^{T}$ represent the coefficients of $z^{2}$ and $z\overline{z}$ of $f_{0}(z, \overline{z})$. Then, we have

Since $i\tau^{*}\omega^{*}$ represents the eigenvalue of $\mathcal{A}(0)$ and $p(0)$ denotes the corresponding eigenvector, we have

Now, substituting the equations (5.29) and (5.35) in equation (5.31) yields

which eventually leads to

Solving the above matrix equation (5.40), we obtain

and

where

In a similar manner, substituting the equations (5.30) and (5.36) in equation (5.32) we obtain

Solving the above matrix equation (5.41), we obtain

where

Thus, we can determine $W_{20}(\theta)$ and $W_{11}(\theta)$ from the equations (5.29) and (5.30). Furthermore, we can also compute $g_{21}$. In order to determine properties of Hopf-bifurcating periodic solution of the system (2.1)-(2.4) at the critical value $\tau^{*}$ on the center manifold, we compute the following quantities:

Here, $\zeta_{2}$ determines the direction of Hopf bifurcation, $\beta_{2}$ determines the stability of the bifurcating periodic solution and $T_{2}$ determines the period of the Hopf-bifurcating periodic solution. Therefore, following the results presented in ^[75], we summarize the properties of the Hopf bifurcation at the critical value $\tau^{*}$ in the following theorem:

Theorem 5.1. Let us assume that the system (2.1)-(2.4) undergoes Hopf bifurcation at the coexisting steady state $E_{*} = (N_{*}, P_{*})$ when $\tau = \tau^{*}$. Then the following results hold:

(i) If $\zeta_{2} > 0\; (\zeta_{2} < 0)$, then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for $\tau > \tau^{*}\; (\tau < \tau^{*})$.

(ii) If $\beta_{2} < 0\; (\beta_{2} > 0)$, then the bifurcating periodic solutions are stable (unstable).

(iii) If $T_{2} > 0\; (T_{2} < 0)$, then the period of the periodic solutions increases (decreases).

Remark 5.1. In order to numerically compute the formulae derived in this section regarding the properties of the Hopf-bifurcating periodic solutions, we have to make use of the non-negative $k$-value for which we will get the Hopf bifurcation threshold $\tau^{*}$ for a particular parameter set. Further, to match the numerical results presented in the preceding section, we have calculated the numerical values for these formulae. Our computations give $C_{1}(0) = -0.8435-0.3956i$, $\zeta_{2} = 26.3009 > 0$, $\beta_{2} = -1.6870 < 0$ and $T_{2} = 6.0453 > 0$. Detailed computations of these quantities are presented in Appendix Ⅱ. Therefore, we can expect that the Hopf bifurcation is forward in direction, i.e., the bifurcating periodic solutions arise for $\tau > \tau^{*}_{0}$; emerging periodic solutions are stable and their period increase. These conclusions can be easily verified from the figures presented in the preceding section.

6.   Numerical simulations

In this section, we carry out numerical simulations in order to demonstrate the impact of gestation delay on the spatial structures exhibited by the prey and predator populations and to corroborate the theoretical results obtained. Numerical simulations were performed over a square spatial domain $\Omega = [0,200]\times[0,200]$ along with grid sizes $\triangle x = 1.0$ and $\triangle y = 1.0$. We have considered time increment as $\triangle t = 0.01$ which is adequate in taking care of the numerical artifacts. For each numerical result, we have used small random perturbations about the homogeneous coexistence steady state as initial population distribution by using the usual Gaussian white noise terms. For the numerical integration of the reaction part, we employ the first order Euler scheme ^[38,59]. We discretize the Laplacian term at the lattice site $(i, j)$ by the following formula

where the matrix elements of $C_{l}$, $C_{r}$, $C_{d}$ and $C_{u}$ take the value $1$ except at the boundary ^[59]. Further in order to guarantee the zero-flux boundary conditions, we define $C_{l}(i, j)N(i-1, j)\equiv N(i+1, j)$ and we use similar definitions for $C_{r}$, $C_{d}$ and $C_{u}$ ^[59]. As we have observed that both the prey and predator populations exhibit patterns with the similar qualitative property and a change in scale, we chose to present the patterns only for the prey population. Also, in this section we restrict ourselves only to the case of strong Allee effect in prey growth. For this reason, we chose the following parameter values: $r = 1.0$, $K = 10.0$, $m = 1.0$, $b = 0.5$, $c = 1.0$, $q = 0.35$, $\delta = 0.0425$, $s = 2.0$, $d_{1} = 0.15$ and $d_{2} = 10.0$ ^[20,24]. Here we would like to mention that we have chosen the value of the parameter $b$ different from the value taken in ^[20,24] in order to make $N^{\text{crit}}$ visible better. For these parameter values, we obtain $\frac{b^{2}r}{K} = 0.025$, $br = 0.5$ and $\frac{r^{2}(K-b)^{2}+4br^{2}K}{4rK} = 2.756$. Therefore, our choice of the parameter value $m = 1.0$ gives rise to the strong Allee effect in prey growth.

Firstly, we aim to provide numerical results which support our theoretical findings on the Hopf bifurcation and the associated properties derived in the preceding sections. For this purpose, we have considered $a = 4.0$ along with the other parameter values mentioned above. In this case, the corresponding non-delayed spatial system exhibits spatially homogeneous distribution of both the populations at the homogeneous coexistence steady state $E_{*} = (3.665, 3.015)$. For the sake of brevity, we restrict ourselves from displaying these numerical results. Thus, this choice of the parameter values acts as a desired parameter set to determine the Hopf bifurcation threshold induced by $\tau$ numerically. Our numerical investigations suggest that for lower values of $\tau$, the system (2.1)-(2.4) retains the homogeneous distribution of both the populations as it was the case for the non-delayed system and we chose not to display them here. Gradually increasing the value of $\tau$, we have achieved the space independent periodic in time population distribution through Hopf bifurcation. In order to check this space independence of the periodic distributions, we present the plots of the limit cycles for $\tau = 1.7$ in case of spatially averaged densities and densities at a particular spatial location $(100,100)$ in Figure 3. From these two figures, we can easily observe that the limit cycles are exactly the same in both the cases. Also, one can get similar results for the plots of the density of prey against that of predator at other spatial locations instead of $(100,100)$. This eventually confirms the space independence of the periodic patterns. Through numerical simulations we have arrived at the Hopf bifurcation threshold which is approximately $\tau^{*} = 1.275$ in this case and it is actually the same as the threshold $\tau^{*}_{0}$ for the corresponding non-spatial system. Therefore, in this case we can directly compute $\tau^{*}$ by taking $k = 0$ and any positive $k$ does not have any role in determining $\tau^{*}$. For the chosen parameter set, non-delayed counterpart of the considered spatial model and its non-spatial version exhibit that the coexistence steady state is stable in both the cases. Therefore, introduction of diffusion in non-delayed temporal model does not play any prominent role in the resulting dynamics for this set of parameter values. This is exactly the reason why any positive $k$ does not contribute in determining the threshold $\tau^{*}$ for the considered spatial model. Therefore, computations of the quantities regarding the properties of the Hopf bifurcation given in Section 5 follow directly from the computed values given in Remark 2 in this case. Hence, we can conclude from the computed values of the relevant quantities that the Hopf bifurcation is forward, i.e., the bifurcating periodic solutions emerge when $\tau > 1.275$, and these solutions are stable for our delayed spatial model.

In order to compare these periodic solutions emerging for the spatial model with that of the corresponding non-spatial model, we present Figure 4. Figure 4(a) shows that the limit cycles in both the cases are almost the same whereas Figures 4(b) and (c) exhibits the small deviations between these two limit cycles in different regions. We remark that this slight differences in these two limit cycles arise due to the presence of the spatial components. Even bigger value of $\tau$ eventually leads to the spatiotemporal chaotic pattern and these scenarios are depicted in Figure 5. Figure 5(a) exhibits the chaotic pattern for the prey population and Figures 5(b) and (c) give us assurance regarding the chaotic nature of this pattern. Also, we have cross-checked the chaotic patterns by examining the sensitivity to initial distributions, estimation of Lyapunov exponent and power spectrum ^{[78,79,80,81]} but did not present the results here for the sake of briefness. All these transitions are clearly visible from the bifurcation diagram presented in Figure 6. This bifurcation diagram (that is, Figure 6) has been prepared by plotting the global maximum and minimum of the average prey density for stable homogeneous and periodic solutions and the local maxima and minima of the average prey density for chaotic solutions after discarding the initial transients ^[57]. Also, from this figure we can observe that if $\tau$ is sufficiently large beyond the Hopf bifurcation threshold $\tau^{*}$ (or $\tau^{*}_{0}$) but less than the heteroclinic bifurcation threshold $\tau_{HO}$ then the periodic solution for the corresponding non-spatial system turns out to be chaotic for the spatial counterpart. For this parameter set, chaotic distribution for both the populations emerges when $\tau > 1.96$. Comparing Figure 2 with Figure 6, we can easily notice that incorporation of diffusion enhances the persistence of both the species as in the $\tau$-range where the non-spatial system exhibits extinction scenario for both the species, the corresponding spatial system leads to the chaotic coexistence.

Now we investigate numerically the impact of gestation delay on the stationary Turing patterns. For this reason, we chose the parameter value $a = 3.5$ along with all other parameter values the same as above. From the conditions for Turing instability presented in ^[20,24], one can easily compute that these parameter values lie inside the Turing domain. A brief discussion about the analytical conditions for Turing instability and the justification for the fact that these parameter values lie inside the Turing domain of the non-delayed spatial system are presented in Appendix Ⅰ. For these parameter values, the non-delayed version of the system (2.1)-(2.4) exhibits stationary cold spot patterns for both the species. But for the sake of brevity of this paper, we did not incorporate these patterns here. Figure 7 represents the effect of delay on this stationary cold spot pattern. From Figure 7(a), we can observe that for lower values of $\tau$ the considered delayed system (2.1)-(2.4) preserves the cold spot pattern. If we increase the value of $\tau$ this cold spot pattern changes into mixture of stripes and cold spots pattern which is exhibited in Figure 7(b). Both these patterns are stationary in nature which can be noticed in Figures 8(a) and (b) for the corresponding temporal evolution of the spatial average densities of both the species. Also, we would like to mention here that for the stationary patterns of both the populations we have noticed almost one-to-one correspondence and accordingly we exhibit only the patterns for the prey population. Another interesting dynamical feature arises in this transition of one stationary pattern to another which is termed as the "loss of monotonicity". This feature is evident from the Figures 8(a) and (b). In Figure 8(a), we can see that the spatial average densities of both the populations approaches the stationary state monotonically (monotonically increasing in case of prey population and monotonically decreasing in case of predator population). Whereas in case of the stationary mixture pattern, there are initial fluctuations in the temporal evolution of the spatial average densities of both the populations before settling to the stationary state and this can be observed in the Figure 8(b). Therefore, we conclude that the transition of stationary patterns for our considered system (2.1)-(2.4) happens through the loss of monotonicity of the spatial average densities. Further gradual increments in $\tau$ leads to quasi-periodic and spatiotemporal chaotic patterns for both the populations respectively. These scenarios can be observed in Figures 7(c) and (d) for the prey population. The corresponding spatial average densities of both the populations are exhibited in Figures 8(c) and (d). In order to characterize these patterns efficiently we have incorporated the corresponding phase diagrams of the spatial average densities in Figure 9. Since the phase diagram presented in Figure 9(a) shows that the spatially averaged solution trajectory stays confined in an annular region, we have termed the corresponding pattern as the quasi-periodic pattern though it looked quite similar to a stationary pattern. Also, the annular domain will not change with the increase in the length of time interval. At this end, we would like to specify that all the patterns reported in this paper are self-organized as we have studied the considered system in homogeneous environment and there are no other external interferences. For the other choices of the parameter values in the Turing domain for the corresponding non-delayed spatial system, we have observed similar kind of transition of stationary patterns to spatiotemporal chaotic patterns due to the effect of the gestation delay. But we did not incorporate them here deliberately in order to make this paper concise.

7.   Discussion

Various research works were carried out to derive analytical expressions for the relevant quantities to determine the stability of the Hopf-bifurcating periodic solutions for delayed temporal and spatiotemporal prey-predator models ^{[30,39,40,41,46,48,50,51,54,55,56,58,59,60]}. In some of them, preliminary numerical simulations were carried out to explore the change in pattern formation scenario due to the presence of time delay. So far as our knowledge goes, there is no literature where the analytical results for the stability of Hopf-bifurcating periodic solutions for the non-spatial and spatial models have been verified numerically and the results are compared with the help of bifurcation diagrams for these two types of models.

In this paper, we have considered a delayed spatial prey-predator system with monotonic functional response with saturation to describe grazing phenomena and density dependent death rate for predator along with the additive Allee effect in prey growth. The consideration of density dependent death rate for predator is a realistic assumption and it is crucial for pattern formation as it allows the system to satisfy Turing instability conditions for certain choice of parameter values while models without it cannot produce Turing patterns ^[82]. The non-delayed version of our model cannot produce Turing patterns in absence of density dependent death rate for predators (i.e., $\delta = 0$) ^[24]. Also, one can go through the Note presented in Appendix Ⅰ for a brief justification of this observation. But it produces various kinds of Turing patterns in presence of density dependent death rate for predators (i.e., $\delta > 0$) ^[20,24]. We have incorporated the discrete time delay in the predator's growth to take care of gestation time period for predator species and studied the effect of time delay on the spatiotemporal pattern formation. The well-posedness of the problem in terms of existence, uniqueness, non-negativity and boundedness of the solutions has been established to ensure the mathematical model is appropriate. Theoretically we have derived the Hopf bifurcation threshold for the model under suitable assumptions and computed the analytical expressions in order to determine the properties of Hopf bifurcating periodic solutions with the help of center manifold reduction and normal form theory with discrete time delay as bifurcation parameter. Although the computed results are for spatial system, we have discussed on how one can actually obtain the associated results for the corresponding non-spatial system by taking $k$ appropriately (i.e., $k = 0$). Accordingly, we have discussed the Hopf bifurcation threshold and the corresponding properties of the limit cycle generated through Hopf bifurcation. Also, we have observed another bifurcation, which is a global bifurcation, the heteroclinic bifurcation for larger value of delay through which the periodic solution disappears and we have found the extinction scenario for both the species in the non-spatial model. Further, numerical illustrations have been provided in order to corroborate our theoretical results and to investigate the impact of time delay on the stationary patterns.

Our present study confirms that the stable Hopf bifurcating periodic solution changes to homogeneous in space and oscillatory in time solution and remains stable in the presence of spatial component. Further, we observe the emergence of the spatiotemporal chaotic distribution for both the species in presence of sufficiently large value of the gestation delay. Time delay changes the spatially uniform distribution of the populations to chaotic distribution through the periodic distribution. This periodic distribution emerges through Hopf bifurcation with time delay as the bifurcation parameter. Delay induced spatiotemporal pattern and enhancement of parametric domain for spatiotemporal chaotic pattern are available in literature ^{[36,37,38,59]} but its emergence from Hopf-bifurcating periodic solution is not well explained in the existing literature. For the model under consideration and for chosen parameter values, we have observed Hopf bifurcation threshold for the spatial system is the same as that for the corresponding non-spatial system. This phenomenon occurs since for the chosen parameter values the corresponding non-delayed temporal system admits asymptotically stable unique coexistence steady state and the corresponding non-delayed spatial system exhibits uniform distribution of populations at this stable steady state. Therefore, incorporation of diffusion does not play a prominent role on the destabilization of the dynamics when the parameter values are away from the Hopf bifurcation threshold for the temporal model and Turing bifurcation threshold for the spatial model. Hence, the Hopf bifurcation threshold (which has been calculated by considering $k = 0$) is the same in both non-spatial and spatial models. In this regard, we remark that positive $k$-values come into picture in determining the Hopf bifurcation threshold when the unstable positive steady state for the non-delayed temporal system is stabilized by the incorporation of spatial components into the model. Also, our numerical results have indicated that incorporation of spatial components enhances the coexistence for both the species and this claim can be easily verified by comparing the results presented in the Figures 2 and 6. From these two figures, one can easily understand that the chaotic coexistence results in the spatial model where the corresponding non-spatial counterpart exhibits extinction of both the species. This chaotic coexistence appears in the spatial model due to the localized disappearance and restoration at other nearby favorable patches in the natural habitat. Another interesting fact is that we have obtained the chaotic pattern for much higher value of the ratio of predator's diffusion coefficient to that of the prey population. And we can safely associate the reason for this phenomenon to the presence of the time delay as it is generally a common practice in literature to take this ratio to be equal to unity (especially $d_{1} = d_{2} = 1$) in order to produce the chaotic distribution of populations ^{[6,8,20,78,83]}.

Another interesting finding in this study is that, apart from the emergence of the chaotic patterns, time delay can actually alter a stationary pattern to another stationary pattern. For smaller values of the time delay the considered system retains the stationary pattern exhibited by the corresponding non-delayed system in its Turing domain. With the increment in time delay the stationary pattern gets converted into another one which eventually turns into chaotic pattern for sufficiently large time delay. In our study, we have shown this transition where cold spot pattern turns into a stationary mixture pattern (mixture of cold spots and stripes). In this case, the mixture pattern eventually settles into a chaotic pattern through the quasi-periodic one with the increase in the magnitude of time delay. Further, we have identified the loss of monotonicity property for the spatially averaged densities in this transition from one stationary pattern to another. Generally, it is believed that both the local and external factors acts together in shaping the spatial pattern formation for ecological community ^[84,85]. By chemical-physical limitations, environmental heterogeneity plays a crucial role in shaping population dynamics as one of the most common external component ^[86]. However, there exist sufficient amount of empirical manifestations where it is often very difficult to justify the observed spatial patterns through environmental heterogeneity solely ^[57,87,88]. This fact eventually leads to the well-known concept of "self-organized" spatiotemporal pattern formation. This type of patterns emerges due to internal factors alone in an absolutely homogeneous environment ^{[3,83,89,90,91]}. Therefore, all the patterns reported in this study has been classified as the self-organized patterns since we have considered our system in a uniform environment and in absence of any other external factors like periodic and/or fluctuating environmental conditions and stochasticity.

The works presented in this paper give rise to various relevant and interesting problems to investigate regarding prey and predator interactions. An immediate future direction is to study the problem with predator density dependent functional responses, namely, Beddington-DeAngelis and Michaelis-Menten types of functional responses. Another important and interesting problem will be to consider the system in presence of both time delay and external influencing factors (especially the environmental heterogeneity) and study the corresponding dynamical behaviors. Also in reality, predator's movement is generally guided by some specific ingredients passed off by prey species such as pheromones and odour etc. ^[92]. Therefore, it will be a problem of immense potential to study the problem presented in this paper in presence of prey-taxis. We will take up and explore these problems as our future studies.

Acknowledgments

The first author gratefully acknowledges the financial support provided by Science and Engineering Research Board, Government of India for pursuing his post-doctoral research. The authors express their gratitude to the learned reviewers for their insightful comments and suggestions.

Conflict of interest

The authors declare there is no conflict of interest.