Fear effect in a three-species food chain model with generalist predator

1.   Introduction

Interactions among species, such as host-pathogen, plant-herbivore, herbivore-carnivore, phytoplankton-zooplankton, and host-parasitoid, play a pivotal role in driving species evolution. These intricate connections maintain biodiversity within ecosystems. Predator-prey interactions, in particular, serve as fundamental components of ecological food chains, illustrating the interdependence of organisms through feeding relationships. Each level within a food chain represents a distinct trophic level, and energy transfers from lower trophic levels to higher ones. The study of predator-prey interactions has been a focal point of scientific inquiry since the seminal works of Lotka ^[1] and Volterra ^[2]. A wealth of research articles exploring the complex dynamics of predator-prey interactions is available in the literature ^{[3,4,5,6,7,8]}. These investigations delve into the intricate dynamics governing predator-prey relationships, shedding light on the mechanisms shaping ecosystem stability and species coexistence.

Predators exhibit varying degrees of specificity in their prey choices, with some displaying a specialized diet while others are more flexible in their food preferences. Certain species, like certain carabid and staphylinid beetles as well as certain types of ants, demonstrate a broad range of prey choices. However, the size of potential prey often plays a crucial role in determining the suitability of a target for predation ^[9]. For instance, the larvae of the green lacewings Chrysoperla carnea prefer aphids as their primary food source but can also consume other small insects and their eggs ^[10]. Ladybirds, such as Coccinella septempunctata, primarily rely on various species of aphids but supplement their diet with different insects, particularly Thysanoptera, as well as fungal spores and pollen grains ^[11]. While some coccinellids, like Coleomegilla maculata, exhibit more flexibility in their dietary choices ^[12,13], others like Rodolia (Vedelia) cardinalis are less adaptable and tend to have narrower food preferences. It is worth noting that many predators labeled as generalists might still display varying degrees of specialization in their prey selection. Additionally, predators with a broad range of prey options typically do not exhibit significant increases in population or consumption in response to a single prey species unless that specific prey constitutes a substantial portion of the available food sources. Overall, the feeding behavior of predators is often shaped by a combination of factors including prey size, availability, and the predator's adaptability, leading to varying degrees of specialization or generalization in their diet preferences.

Theoretical advancements in ecological modeling have traditionally emphasized specialist predators due to their relatively simpler parameterization. However, efforts have been made to incorporate and understand the dynamics of generalist predators in relation to prey populations through two distinct modeling approaches. The first approach involves studying scenarios where a single predator species interacts with multiple prey species, with each prey serving as a potential food source for the predator ^[8]. Alternatively, the second approach assumes a single prey species that the predator consumes while also obtaining supplemental nutrition from alternative food sources. Researchers such as Erbach et al. ^[4] investigated predator-prey systems incorporating a generalist predator, focusing on phenomena like bistability and the occurrence of limit cycles within the dynamics of these systems. Similarly, Magal et al. ^[5] explored a model involving pests and their natural enemies, including a generalist predator, to understand the dynamics of such systems. Moreover, studies by Mondal et al. ^[14] delved into predator-prey systems featuring a generalist predator, specifically analyzing the influence of cooperative hunting strategies on system dynamics. Their research investigated bistabilities and various forms of local and global bifurcations exhibited by these complex ecological systems. In summary, theoretical research has increasingly acknowledged the significance of generalist predators in ecological dynamics, exploring their interactions within predator-prey systems through sophisticated mathematical models to uncover intricate patterns of stability, oscillations, and bifurcations that characterize these ecological relationships.

Predators play a fundamental role in regulating prey populations by hunting for sustenance. However, recent scientific interest has focused on understanding the broader impacts of predation beyond direct consumption, known as non-consumptive or fear effects, on prey demography. These effects encompass changes in the reproductive and survival rates of prey species due to the perceived threat of predation ^[15,16,17]. Experimental evidence has demonstrated that the mere presence of predators can significantly alter prey demography in addition to direct killing. For instance, Zanette et al. ^[18] conducted a comprehensive field experiment involving Song sparrows over an entire breeding season. They shielded a specific group of the sparrow population from any predation risks but exposed them to the perceived fear of predators, such as predator sounds, while preventing direct killing of the entire population. The results revealed a substantial 40% reduction in the offspring production of Song sparrows due to predation fear. Similarly, Sheriff et al. ^[19] observed a 30% decline in the survival rates of adult female Snowshoe hares attributed to predation fear. Additionally, Suraci et al. ^[20] found that the presence of carnivorous animals led to reduced hunting rates and population sizes of raccoons, subsequently causing a notable increase in the populations of certain crab species and fishes that are preyed upon by raccoons. These studies collectively highlight the significant influence of predation fear on prey demography, showcasing that the mere presence or perception of predator risk can lead to substantial changes in reproductive success, survival rates, and population dynamics of various prey species.

Predator-prey interactions prompt significant alterations in the foraging behavior and habitat selection of animals and birds. These behavioral and physiological adaptations in response to the presence of predators often outweigh the direct impact of predation itself. Prey species employ diverse anti-predation tactics, such as heightened vigilance, alterations in habitat utilization, and physiological adjustments ^[21,22]. Studies have demonstrated that mule deer, when faced with the threat of lion attacks, adjust their foraging activities as a defensive measure ^[23]. Likewise, elk within the Yellowstone park ecosystem exhibit physiological changes in response to wolf predation ^[24]. While relocating to safer habitats may initially enhance survival prospects, such shifts can potentially compromise long-term prey fitness. Habitats deemed less energy-efficient may exert adverse effects on prey demographics ^[21]. Furthermore, prey animals perceiving predation threat demonstrate reduced foraging behavior ^{[20,24,25,26]}. In avian species, fear of predation manifests in decreased reproductive success, including diminished egg production, lower hatching rates, heightened nestling mortality, thereby influencing demographic shifts in bird populations.

The seminal work addressing the impact of fear within a predator-prey system was initiated by Wang et al. ^[7], where the fear of predation was modeled to account for a reduction in the birth rate of prey. Subsequently, a plethora of studies have delved extensively into the effects of predation-induced fear, exploring various traits inherent in both prey and predator species ^{[27,28,29,30,31,32,33,34]}. Zhang et al. ^[34] investigated predator-prey dynamics by incorporating the fear phenomenon alongside prey refuge. Their findings suggested that the fear factor contributes to system stabilization by eliminating persistent periodic oscillations and concurrently reducing predator density. However, higher fear strength consistently led to a decline in the predator population towards extinction, irrespective of refuge size. In a separate analysis, Wang et al. ^[33] explored the effect of fear on the dynamics of a modified Leslie-Gower predator-prey system with prey refuge, inferring that predation fear and prey refuge have divergent effects on species persistence. Mukherjee ^[35] introduced the concept of the cost of fear, which elevates the death rate while reducing the reproduction rate of prey. Additionally, Mondal et al. ^[36] conducted a comprehensive study encompassing hunting cooperation, prey refuge, harvesting, and predation-induced fear in both autonomous and nonautonomous predator-prey systems. They highlighted that the fear of predation diminishes the birth rate and intensifies intraspecies competition among prey.

In natural ecosystems, it is widely observed that species prey upon one while also serving as prey for other species. These interlinked interactions form intricate food chains, creating complex networks within ecosystems. This complexity raises the fundamental question of how the fear experienced by one or more species within this chain can impact the demography of the entire population cascade. The experimental investigation by Suraci et al. ^[20] explored a trophic cascade by manipulating the playback of large carnivores. Their findings showcased a decline in the population of meso-carnivores and an increase in the population of their prey. This suggests that the fear induced by top predators within the food chain explicitly or implicitly influences each species within the ecosystem. Cong et al. ^[3] delved into a three-species food chain system enriched with the fear of predation, analyzing its dynamics. Additionally, Panday et al. ^[37] examined a three-species food chain, incorporating the fear effect in the predation rate of the intermediate predator along with growth rates. Their study observed trophic cascading triggered by fear. Their results highlighted that the fear factor stabilizes the system's dynamics, and the cost associated with fear in the intermediate predator enhances the persistence of species within the ecosystem.

In this study, we investigate a three-species food chain model comprising a prey species, an intermediate specialist predator, and a top predator with a generalist feeding behavior. We incorporate the impact of predation induced fear, accounting for diminished birth rates and heightened intraspecies competition in the model formulation. By employing a comprehensive array of mathematical and numerical techniques, our primary objective is to elucidate the complex dynamics and interactions inherent within this food chain. Specifically, we aim to delve into the intricate effects arising from fear responses and the consequential costs associated with predation, seeking a deeper understanding of their influence on the dynamics of this ecological system.

The subsequent sections of the paper are structured as follows: we describe our proposed model in Sections 2 and 3 delves into a comprehensive exploration of our model and lays out the foundational mathematical concepts. In particular, Section 3.1 rigorously examines the positivity and boundedness of the solution within the proposed system. Following this, Section 3.2 elucidates the identification and characterization of feasible equilibrium points, along with their corresponding stability conditions. Furthering our analysis, Section 3.3 scrutinizes the potential occurrence of local bifurcations under specific parametric configurations, providing an in-depth assessment of the system's behavior. Moving forward to Section 4, we employ numerical examples and employ graphical representations to illustrate and apprehend the dynamic nature exhibited by the proposed model system. Conclusively, Section 5 encapsulates the paper with insightful concluding remarks, synthesizing the findings and implications drawn from the preceding sections.

2.   Mathematical model

Here, we propose a mathematical model for a food chain system of an ecological community. We consider three species; prey, intermediate predator and top predator, and denote their densities by $x(t)$, $y(t)$ and $z(t)$, respectively, at any instant of time $t > 0$. We formulate our model based on the following ecological assumptions.

1) In the absence of predators, the prey species grow logistically with intrinsic growth rate $r$ and carrying capacity $r/d$; $d$ is the density-dependent death rate.

2) Predation of prey species by intermediate predators is quantified by Holling type-Ⅱ functional response, which is given by $\frac{\alpha_1xy}{m+x}$, where $\alpha_1$ is the capture rate and $m$ is the half saturation constant.

3) Intermediate predators are of specialist type, i.e., only the considered prey serves as the food source for them. Therefore, their growth is proportional to the consumed food through predation. Here, we consider $\theta_1$ as the respective gain in the density of intermediate predators; the value of $\theta_1$ lies between $0$ and $1$ for ecological reasons.

4) The density of intermediate predators diminishes due to natural mortality as well as intraspecies competition. We denote the natural mortality and intraspecies competition coefficient of the intermediate predators by $\delta_1$ and $\delta_2$, respectively. Their density also depletes due to predation by top predators that is modeled by the bilinear functional response, $\alpha_2yz$.

5) Predation has positive feedback on the growth of top predators. Therefore, the growth rate is proportional to the functional response with proportionality constant $\theta_2$. The density of top predators depletes due to natural mortality at a rate $\delta_3$.

6) Unlike intermediate predators, the top predators are generalist in nature, i.e., even in the absence of intermediate predators, the top predator population increase due to availability of food from other sources. This growth is captured by Beverton-Holt like function ^[4,14], and is given by $\frac{\eta z}{1+\eta_0z},$ where $\eta$ is the per capita reproduction rate, and $\eta_0$ is the density-dependent strength.

With these assumptions, a three-species food chain model is obtained as follows:

where the initial population densities are taken as positive, i.e., $x(0) > 0, \ y(0) > 0 \ \text{and} \ z(0) > 0$. In system (2.1), the value of parameter $\eta$ should be more than that of $\delta_3$, i.e., the natural mortality rate of top predators. That is to say, the top predators can persist in the ecosystem even in the absence of intermediate predators.

Predation-induced fear in prey species is a multifaceted phenomenon that significantly impacts their population dynamics beyond direct predation. Studies such as Zanette et al. ^[18] have highlighted its role in diminishing birth rates, while research by Allen et al. ^[38], Clinchy et al. ^[39], and Sheriff et al. ^[40] emphasizes its influence on foraging behavior and adult survival. This fear alters the intrinsic growth rate of prey species, influencing their overall population dynamics. Furthermore, the perception of predation risk prompts prey species to adopt anti-predation strategies, leading to reduced foraging, seeking safer habitats, and allocating more energy and time to vigilance. These behavioral adaptations often result in the congregation of prey in specific areas, forming high-density patches. Consequently, this aggregation intensifies intraspecific competition among prey due to limited resources in these localized areas, compounded by reduced foraging over larger territories. To mathematically model this, we have incorporated the cost of fear by adjusting the parameters governing prey population growth and intraspecific competition. Specifically, we multiply the growth rate ($r$) of prey species by a decreasing function $1/(1+\omega_1 y)$, which accounts for the reduction in reproductive capacity due to fear-induced changes in behavior. Simultaneously, we consider the impact of fear on intraspecific competition by applying an increasing function ($1+\omega_2 y$) to the intraspecific competition coefficient ($d$). Here, $\omega_1$ and $\omega_2$ represent fear parameters, enabling us to quantify and integrate the effects of fear on both reproductive capacity and intraspecific competition within the mathematical model governing the dynamics of prey populations.

Predation-induced fear not only impacts the growth rate of prey species but also exerts a substantial influence on the behavior and predation activities of intermediate predators. Experimental study by Gordon et al. ^[41] demonstrated that the presence of apex predators, such as dingoes, suppresses the foraging and predation activities of intermediate predators like feral cats. This suppression alleviates perceived predation risk for smaller prey species like desert rodents. Areas with a higher abundance of dingoes exhibited reduced feral cat density, allowing the proliferation of rodent populations. Conversely, regions with fewer dingoes showed lower rodent numbers alongside elevated feral cat densities. Similarly, in a study by Suraci et al. ^[20], they have induced fear responses in mesopredators, like raccoons, through auditory stimuli mimicking large carnivores. This induced fear led to a significant reduction in foraging activity and increased vigilance among raccoons. The resultant decline in raccoon foraging had cascading effects on the ecosystem, benefiting the prey species of mesopredators such as intertidal crabs, fish, worms, and red rock crabs. Moreover, in the Yellowstone National Park ecosystem, the reintroduction of wolves, as apex predators, significantly altered the behavior and grazing patterns of elk, illustrating the far-reaching effects of apex predators on the ecosystem's trophic dynamics ^[42]. In our modeled system, both the growth term and predation rate of intermediate predators are modified by a factor $\frac{1}{1+\omega_3 z}$, where $\omega_3$ is the fear parameter and $z$ represents the density of top predators. This adjustment accounts for the fear-induced changes in the behavior of mesopredators due to the presence of top predators. As the fear factor or density of top predators increases, the foraging and predation activities of mesopredators decrease. Additionally, the food obtained and consumed by the mesopredator population contributes to their growth. Therefore, the growth term in the mesopredator population is proportional to a modified functional response with a proportionality constant $\theta_1$, incorporating the same factor $\frac{1}{1+\omega_3 z}.$ Furthermore, the intraspecific competition among intermediate predators is influenced by predation-induced fear of top predators. Hence, the intraspecific competition coefficient ($\delta_2$) is adjusted by an increasing function $(1+\omega_4z)$, where $\omega_4$ is the fear parameter. This modification of the model formulation captures the intricate interplay between fear-induced behavioral changes in intermediate predators, their predation activities, growth, and intraspecific competition, elucidating their impacts on trophic interactions within the ecosystem. The schematic representation of these modifications is visually depicted in Figure 1, and the model equations are obtained as follows:

with initial condition $x(0) > 0, \ y(0) > 0 \ \text{and} \ z(0) > 0$.

3.   Analysis of system (2.2)

3.1.   Positivity and boundedness of solutions

In ecological modeling, it is crucial to ensure that the formulated model not only captures the interactions among species but also reflects real-world ecological principles, such as population stability, boundedness, and the absence of infinite growth or extinction.

Theorem 1. System (2.2) exhibits a unique solution that is positive for all $t \geq 0$.

Proof. As the functions on the right-hand side of system (2.2) are continuous as well as locally Lipschitz in the positive octant $\mathbb{R}^3_+$, solution of system (2.2) exists uniquely in the interval $[0, T)$, where $0 < T \leq \infty$ ^[43].

System (2.2) can also be written as

Therefore,

where

Thus, $x(t), \ y(t), \ z(t) > 0$, for all $t\geq0$, whenever $x(0)$, $y(0)$ and $z(0)$ are positive. Hence, the solution trajectory of the system (2.2) that originates within the positive octant lies there indefinitely.

Theorem 2. Solutions of system (2.2) are uniformly bounded.

Proof. Analyzing the first equation of system (2.2), we get

Second equation of system (2.2) implies

Define $u = \theta_2 y+z$. Therefore,

Let $\gamma = \min\{\delta_1, \delta_3\}$, then we have

Above differential inequality implies

Therefore, $0 < u(y(t), z(t)) \leq \left(\frac{\theta_2\theta_1^2\alpha_1^2}{4\delta_2}+\frac{\eta}{\eta_0}\right)+\epsilon$, for any $\epsilon > 0$ as $t\rightarrow \infty$. Thus, each solution of (2.2) is eventually attracted to the region:

Therefore, all the dynamical variables presented in the system (2.2) are bounded.

3.2.   Equilibrium and stability analysis

3.2.1.   Equilibria of the system

Due to the nonlinear nature of the predator-prey model represented by system (2.2), exact analytical solutions are not feasible. Therefore, to understand the long-term behavior of the system, we focus on identifying equilibrium points. These points are obtained by setting the growth rates, as described by the differential equations of the model system, equal to zero. Equilibrium points represent steady states where population changes cease. It can be easily seen that the system (2.2) exhibits three axial equilibria; $E_0(0, 0, 0)$, $E_1(\hat{x}, 0, 0)$, and $E_2(0, 0, \hat{z})$; a planar equilibrium $E_3(\hat{x}, 0, \hat{z})$, where $\hat{x} = r/d$ and $\hat{z} = (\eta-\delta_3)/\eta_0\delta_3 > 0$, as $\eta > \delta_3$. Another planar equilibrium with $z$ component zero can be obtained by putting $z = 0$ in the equilibrium equations of system (2.2), which leads to the following two algebraic equations whose intersection point in the feasible region gives the top predator-free equilibrium $\tilde{E}$:

● The curve (3.1) passes through $(x, y) = (r/d, 0)$ and for $x = 0$, it cuts the positive $y$-axis at some point. Further, the slope of the curve at $(r/d, 0)$ is negative.

● The curve (3.2) cuts the $x$-axis at $x = m\delta_1/(\theta_1\alpha_1-\delta_1)$ and the $y$-axis at $y = -\delta_1/\delta_2$. Also, it has a horizontal asymptote $y = (\theta_1\alpha_1-\delta_1)/\delta_2$ and a vertical asymptote $x = -m$. Moreover, the slope is always positive on the right side of the vertical asymptote.

Therefore, it can be concluded that the curves (3.1) and (3.2) intersect uniquely in the feasible region if $\frac{m\delta_1}{\theta_1\alpha_1-\delta_1} < \frac{r}{d}$; the intersecting point gives the $x$ and $y$-components of the equilibrium $\tilde{E}(\tilde{x}, \tilde{y}, 0)$.

Now, in order to obtain a feasible interior equilibrium, we first write the variable $y$ in terms of $z$ from the third equilibrium equation of system (2.2) to get,

Clearly, $y$ will be positive if $z > \frac{\eta-\delta_3}{\eta_0\delta_3} = \hat{z}$. Notably, $g(\hat{z}) = 0$ and $g'(z) > 0$. Putting this value of $y$ in the first two equilibrium equations of the system (2.2), we get the following two isoclines

● For curve (3.4), at $z = \hat{z}$, $x = \frac{m(\delta_1+\alpha_2\hat{z})}{\frac{\theta_1\alpha_1}{1+\omega_3\hat{z}}-(\delta_1+\alpha_2\hat{z})}$, and the slope $\frac{dx}{dz} > 0$ in the positive quadrant.

● Putting $z = \hat{z}$ in curve (3.5), we obtain $x = \frac{r}{d} = \hat{x}$. At the point $(\hat{x}, \hat{z})$, the slope $\frac{dx}{dz}$ of the curve (3.5) is negative.

Note that, if $\frac{m(\delta_1+\alpha_2\hat{z})}{\frac{\theta_1\alpha_1}{1+\omega_3\hat{z}}-(\delta_1+\alpha_2\hat{z})} < \frac{r}{d}$ and the slope of the isocline (3.5) remains negative in the positive quadrant, then the two isoclines intersect uniquely, and the intersecting point gives the $x$ and $z$ components of the interior equilibrium $E^*(x^*, y^*, z^*)$. By using $x^*$ and $z^*$ in Eq (3.3), the $y$ component of $E^*$ can be obtained.

3.2.2.   Local stability

In this section, we conduct a local stability analysis to evaluate the behavior of the system (2.2) near its equilibrium points. This analysis helps characterize whether the system tends towards or moves away from an equilibrium when initiated close to, but not precisely at, that particular point. Specifically, we investigate whether the equilibria are locally asymptotically stable, meaning that nearby initial conditions lead the system to approach the equilibrium point as time ($t$) progresses towards infinity. To perform the local stability analysis of the feasible equilibria in the system (2.2), we compute the Jacobian matrix, a matrix of partial derivatives, which provides insights into the system's behavior around these points. The Jacobian matrix corresponding to the system (2.2) is obtained as follows:

where

Theorem 3. The axial equilibria $E_0$, $E_1$ and $E_2$, and the planar equilibrium $\tilde{E}$ are always unstable.

Proof. Eigenvalues of the Jacobian matrix $\boldsymbol{J}$ evaluated at the equilibrium $E_0$ are $r$, $-\delta_1$ and $\eta-\delta_3$; positivity of one eigenvalue makes the equilibrium $E_0$ unstable.

The matrix $\boldsymbol{J}$ after evaluation at the equilibrium $E_1(\hat{x}, 0, 0)$ gives the eigenvalues as $-r$, $\frac{\theta_1\alpha_1\frac{r}{d}}{m+\frac{r}{d}}-\delta_1$ and $\eta-\delta_3$. As $\eta > \delta_3$, the equilibrium $E_1$ is unstable.

Eigenvalues of the matrix $\boldsymbol{J}$ corresponding to the equilibrium $E_2(0, 0, \hat{z})$ are $r$, $-\delta_1-\frac{\alpha_2(\eta-\delta_3)}{\eta_0\delta_3}$ and $-\frac{\delta_3}{\eta}(\eta-\delta_3)$. Again, this equilibrium is unstable as one eigenvalue is always positive.

Now, for the equilibrium $\tilde{E}(\tilde{x}, \tilde{y}, 0)$, one eigenvalue of the matrix $\boldsymbol{J}$ comes out to be $\eta-\delta_3+\theta_2\alpha_2\tilde{y}$, which is always positive; the other two eigenvalues are the roots of the following quadratic equation:

where

In view of the positivity of one eigenvalue, the equilibrium $\tilde{E}$ is unconditionally unstable.

Theorem 4. The equilibrium $E_3(\hat{x}, 0, \hat{z})$ is locally asymptotically stable if the following condition is satisfied

Proof. Eigenvalues of the matrix $\boldsymbol{J}$ evaluated at the equilibrium $E_3(\hat{x}, 0, \hat{z})$ are obtained as $-r$, $-\frac{\delta_3}{\eta}(\eta-\delta_3)$, and $\frac{1}{1+\omega_3 \hat{z}}\frac{\theta_1\alpha_1\hat{x}}{m+\hat{x}}-\delta_1-\alpha_2\hat{z}$. Clearly, two eigenvalues are negative, therefore sign of the third eigenvalue will determine the stability of this equilibrium. Hence, $E_3$ is stable if $\frac{\theta_1\alpha_1\hat{x}}{m+\hat{x}} < (1+\omega_3\hat{z})(\delta_1+\alpha_2\hat{z})$.

Theorem 5. The equilibrium $E^*(x^*, y^*, z^*)$, if exists, is locally asymptotically stable if and only if the following conditions are satisfied:

where $J_i$'s ($i = 1, 2, 3$) are defined in the proof.

Proof. Evaluating the Jacobian $\boldsymbol{J}$ at the point $E^*$, we get the following matrix:

where

Characteristic equation corresponding to matrix $\boldsymbol{J^*}$ is obtained as follows:

where

By the Routh-Hurwitz criterion, one can infer that the interior equilibrium $E^*$ is locally asymptotically stable if and only if $J_1 > 0$, $J_3 > 0$ and $J_1J_2-J_3 > 0$.

All the feasible equilibria exhibited by model system (2.2) with the condition for existence and their stability are mentioned in the following table.

3.3.   Bifurcation analysis

3.3.1.   Transcritical bifurcation

In a transcritical bifurcation, a dynamic system undergoes a critical transformation where two equilibria within the system interchange their local stability characteristics as a specific parameter is systematically altered. This transition marks a pivotal point where the stability properties of these equilibria switch, leading to significant alterations in the system's behavior. Evaluating the Jacobian matrix $\boldsymbol{J}$ at the equilibrium $E_3(\hat{x}, 0, \hat{z})$, we get

The equation $\frac{1}{1+\omega_3\hat{z}}\frac{\theta_1\alpha_1r}{r+md}-\delta_1- \frac{\alpha_2(\eta-\delta_3)}{\eta_0\delta_3} = 0$ will determine a critical value of $\eta$ (say $\eta^*$), at which the matrix $\boldsymbol{J}_{E_3}$ has a zero eigenvalue. Let $\mathcal{U} = [u_1 \; u_2 \; u_3]^T$ and $\mathcal{V} = [v_1 \; v_2 \; v_3]^T$ be the eigenvectors corresponding to the zero eigenvalue of the matrices $\boldsymbol{J}_{E_3}$ and $\boldsymbol{J}_{E_3}^T$, respectively, where

Consider $G = [g^1 \; g^2 \; g^3]^T$, where

Now, the transversality conditions are given by

where all the notations are the same as in Theorem $1$ of Section $4.2$ of ^[44]. Thus, system (2.2) exhibits degenerate transcritical bifurcation ^[45] around the equilibrium $E_3$ at $\eta = \eta^*$. Note that for non-degenerate transcritical bifurcation, $\mathcal{V}^T D G_{\eta}(E_3;\eta^*)\mathcal{U}$ must be nonzero.

3.3.2.   Hopf bifurcation

In ecological modeling, nonlinear interactions among populations, even in systems with a low level of complexity involving two or three species, can result in intricate dynamical behaviors. Oscillatory patterns, commonly observed in population dynamics, are indicative of complex dynamics within ecological systems. The occurrence of oscillations or the manifestation of a limit cycle typically corresponds to a Hopf bifurcation in the system. The Hopf bifurcation is characterized by the appearance or disappearance of a periodic orbit due to a local change in the stability properties of an equilibrium point. In our analysis, we aim to explore the potential occurrence of a Hopf bifurcation and determine its direction concerning the coexistence equilibrium $E^*$ concerning the parameter $\omega_1$ that induces bifurcation. Specifically, we investigate how changes in the parameter $\omega_1$ influence the stability of the equilibrium point $E^*$ and whether these changes lead to the emergence or disappearance of a periodic orbit. The direction of the Hopf bifurcation signifies whether an equilibrium point becomes stable or unstable, leading to the creation or elimination of periodic behavior within the system.

Suppose there is a critical value of $\omega_1$ $(\omega_1^* \ say)$ at which the conditions $J_1 > 0$ and $J_3 > 0$ hold but $J_1J_2-J_3 = 0$. In this case, the characteristic equation (3.9) becomes

The above equation has two purely imaginary roots, say $\lambda_{1, 2} = \pm\iota\sqrt{J_2}$, and a negative real root, say $\lambda_3 = -J_1$. Assume that at any point $\omega_1$ in the $\epsilon$-neighborhood of $\omega_1^*$, $\lambda_{1, 2} = \gamma_1\pm\iota\gamma_2$. Putting this in Eq (3.9) and separating real and imaginary parts, we get

As $\gamma_2\neq0$, from Eq (3.13), we have

Using this value of $\gamma_2$ in Eq (3.12), we get

Differentiating the above equation with respect to $\omega_1$ and using the fact that $\gamma_1(\omega_1^*) = 0$, we get

Clearly, the left side of the above equation will be nonzero if $\frac{d}{d\omega_1}(J_1J_2-J_3)|_{\omega_1 = \omega_1^*}\neq 0.$ Therefore, we can say that the system (2.2) exhibits Hopf bifurcation around the equilibrium $E^*$ as $\omega_1$ crosses the critical value $\omega^*_1$.

For a clear understanding of the instability behavior, it is needed to obtain the amplitude and the initial period of the periodic solutions. For this, we set $J_3 = \kappa J_1J_2$ in the characteristics equation (3.9). Assuming $\lambda$ as a continuous function of $\kappa$, we can rewrite Eq (3.9) as

At $\kappa = \kappa^* = 1$, $J_3 = J_1J_2$ and Eq (3.14) factorizes into $(\lambda+J_1)(\lambda^2+J_2) = 0$, that gives the roots $\lambda(\kappa^*) = -J_1$ and $\lambda(\kappa^*) = \pm i\sqrt{J_2}$. This assures the occurrence of Hopf bifurcation in system (2.2). Notably, this new parametrization, $0\leq\kappa\leq\kappa^*$ corresponds to $0\leq\omega_1\leq\omega_1^*$, $0\leq\kappa = \kappa^*$ corresponds to $\omega_1 = \omega_1^*$ and $\kappa\geq\kappa^*$ is identical to $\omega_1\geq\omega_1^*$. Setting $\kappa = \kappa^*+\mu^2\xi$, where $\|\mu\| < 1$ and $\xi = \pm 1$, we get $\lambda(\kappa) = \lambda(\kappa^*+\mu^2\xi)$. Expanding by Taylor series about $\kappa^*$ gives

where $'$ represent derivative with respect to $\kappa$. Differentiating both side of Eq (3.15) and simplifying, we get

As $\Re(\lambda(\kappa^*)) = 0$, we have $\Re(\lambda'(\kappa^*)) > \frac{J_1J_2}{2(J_1^2+J_2)} > 0$. Substituting the values of $\lambda(\kappa^*)$ and $\lambda'(\kappa^*)$ in Eq (3.15), we get the following approximation:

From the above equation, we obtain the amplitude and initial period of the oscillatory solutions that occurred along with the loss of stability when $\kappa > \kappa^*$ as $\exp\left(\frac{J_1J_2\mu^2\xi}{2(J_1^2+J_2)}\right)$ and $\frac{2\pi}{\sqrt{J_2}\left(1+\frac{J_1^2 \mu^2\xi}{2(J_1^2+J_2)}\right)}$, respectively, where $\mu = \sqrt{\frac{|\kappa-\kappa^*|}{|\xi|}}$.

3.3.3.   Direction and stability of Hopf bifurcation

By using normal form theory ^[46], we discuss the direction of Hopf bifurcation with stability properties of bifurcating oscillatory solutions of system (2.2). The eigenvectors $W_1$ and $W_2$ corresponding to the eigenvalues $\lambda_1 = i\phi$ and $\lambda_3 = -J_1$, at $\omega_1 = \omega_1^*$, where $\phi = \sqrt{J_2}$, are respectively given by:

with

In view of the following transformations

system (2.2) transforms to

where

with

Thus, the interior equilibrium $E^*$ shifts to the origin $(0, 0, 0)$ for system (3.19), and the corresponding Jacobian matrix is obtained as

where $\frac{\partial L^1}{\partial \mathfrak{X}} = \frac{\partial L^2}{\partial \mathfrak{Y}} = \frac{\partial L^1}{\partial \mathfrak{Z}} = \frac{\partial L^3}{\partial \mathfrak{X}} = \frac{\partial L^3}{\partial \mathfrak{Y}} = \frac{\partial L^2}{\partial \mathfrak{Z}} = 0, -\frac{\partial L^1}{\partial \mathfrak{Y}} = \frac{\partial L^2}{\partial \mathfrak{X}} = \phi, \ \text{and} \ \frac{\partial L^3}{\partial \mathfrak{Z}} = D_1$.

The values of $g_{11}$, $g_{02}$, $g_{20}$, $G_{101}$, $G_{110}$, $G_{21}$, $h_{11}$, $h_{20}$, $\phi$, $\phi_{20}$, $\phi_{11}$, and $g_{21}$ are computed by using the following relations:

Solutions of the following equations will give the values of $\omega_{11}$ and $\omega_{20}$:

Then, we compute

Now, using the above values, we calculate the following quantities:

where $\alpha'(0) = \frac{d}{d\omega_1}(\Re\{\lambda_1(\omega_1)\})|_{\omega_1 = \omega_1^*}$ and $\phi'(0) = \frac{d}{d\omega_1}(\Im\{\lambda_1(\omega_1)\})|_{\omega_1 = \omega_1^*}$. The direction of Hopf bifurcation, supercritical or subcritical, is determined by the sign of $\mu_2$. A positive ($\mu_2 > 0$) value indicates a supercritical Hopf, while a negative ($\mu_2 < 0$) value indicates a subcritical Hopf bifurcation. Further, the stability of the oscillatory solutions is characterized by the parameter $\beta$; a negative ($\beta_2 < 0$) value indicates stability and a positive ($\beta_2 > 0$) value indicates instability. Additionally, the period of the oscillatory solutions increases or decreases when $T_2 > 0$ or $T_2 < 0$.

4.   Numerical simulations

In this section, we aim to computationally investigate the dynamics of system (2.2) using MATLAB and MATCONT. Initially, we focus on analyzing the behavior of system (2.2) in the absence of predation-induced fear. Subsequently, we aim to assess the influence of fear on the dynamical characteristics of system (2.2) and the population dynamics of the species within the ecosystem. To facilitate numerical observations, we consider a specific set of parameter values for system (2.2):

It is important to note that unless specifically mentioned, we consistently employ the aforementioned set of parameter values for simulations. To explore diverse dynamics manifested by system (2.2), we systematically vary certain parameters within biologically plausible ranges. Through these numerical analysis, we aim to glean insights into the various behaviors and interactions exhibited by system (2.2) under different parameter configurations. This exploration will provide valuable perspectives on the ecosystem's dynamics and the impact of fear-induced responses on the species' population dynamics within this framework.

4.1.   Sensitivity analysis

In order to address the inherent uncertainties in determining parameter values for system (2.2), we employ global sensitivity analysis, leveraging two statistical techniques: Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCs) ^[47,48]. LHS involves a stratified sampling approach without replacement, allowing simultaneous variation of multiple parameters in an efficient manner. Meanwhile, PRCC evaluates the strength and direction of correlation between model output and input parameters, yielding values within the interval $[-1, 1]$. Assuming a uniform distribution for the input parameters, namely $\omega_1$, $\omega_2$, $\omega_3$, $\omega_4$, $\alpha_1$, $\eta$ and $\eta_0$, we conduct 50 simulations per LHS of system (2.2). Utilizing the baseline values specified in Eq (4.1) alongside $\omega_1 = 0.5$, $\omega_2 = 0.5$, $\omega_3 = 0.1$, $\omega_4 = 0.5$, $\alpha_1 = 12$, $\eta = 0.6$ and $\eta_0 = 2$. We allow parameters to deviate within a range of $\pm 25\%$ from these nominal values. This approach facilitates a comprehensive exploration of the parameter space and aids in understanding the sensitivity of the model output to variations in these parameters. By systematically analyzing how alterations in parameter values impact the system's behavior, we aim to enhance our comprehension of the model's robustness, uncertainties, and the influence of individual parameters on the overall dynamics of the system.

Figure 2 illustrates the Partial Rank Correlation Coefficients (PRCC) values assigned to the considered input parameters within model system (2.2), utilizing the density of prey species as the output variable. Notably, parameters exhibiting higher PRCC values exert a more pronounced impact on the density of prey species. Consequently, parameters affecting the density of prey species positively (via positive PRCC values) or negatively (via negative PRCC values) are identified. From the analysis, parameters exhibiting a negative influence on the density of prey species include $\omega_1$, $\omega_2$, $\alpha_1$ and $\eta_0$, while those contributing positively to the density of prey species comprise $\omega_3$, $\omega_4$ and $\eta$. Among these parameters $\alpha_1$, $\eta$ and $\eta_0$ emerge as particularly significant. Identifying these key parameters assumes significance in devising effective control strategies essential for preserving prey species within the ecosystem. This sensitivity analysis suggests that strategies aimed at augmenting (diminishing) parameters associated with positive (negative) PRCC values would effectively enhance the density of prey species in the ecosystem. Understanding the influential parameters and their respective impacts on the population dynamics of prey species provides valuable insights for implementing targeted interventions or management practices crucial for sustaining and conserving the prey species within the ecosystem.

4.2.   Dynamics of system (2.2) in the absence of fear factors

In this section, we initiate simulation by assigning a zero value to all fear factors to investigate the influence of parameters $\eta$ and $\eta_0$ on the dynamics of the system (2.2). This investigation aims to elucidate the effects arising from the inclusion of generalist top predators rather than specialized predators. Upon setting $\eta = 1$, the proposed system demonstrates a unique stable interior equilibrium $E^* = (2.23442, 0.58000, 3.46831)$. The corresponding Jacobian matrix yields eigenvalues of $-0.10878$ and $-0.06073\pm 1.82598i$. Concurrently, all other feasible boundary and planar equilibria manifest as intrinsically unstable. Figure 3 portrays the time series solutions for all system variables and the respective phase portrait. It is notable that as the value of $\eta$ diminishes, the stability of the system is compromised, resulting in instability. When $\eta = 0.5$, Figure 4 illustrates the time series solution and the phase portrait, revealing the periodic oscillations around the co-existence equilibrium. This figure demonstrates that trajectories converging towards the limit cycle, regardless of their initiation from within or outside the cycle.

Subsequently, Figure 5 illustrates the equilibrium curve across a significant spectrum of $\eta$ values, serving to elucidate the repercussions of augmented food availability on the dynamics of the system. In this depiction, stable equilibria are denoted in blue, whereas unstable equilibria are marked in red. Notably, for lower $\eta$ values, the system manifests oscillatory behavior. However, with an increase in $\eta$, the co-existence equilibrium stabilizes via a supercritical Hopf bifurcation occurring at $\eta = 0.84127$ (designated as $H$). To visually comprehend these bifurcations, Hopf bifurcation diagrams are represented in Figure 6. As $\eta$ continues to escalate, the system experiences a transcritical bifurcation at $\eta = 3.71586$ (labeled as $BP$). At this bifurcation point, stability transitions from the internal equilibrium to the intermediate predator-free equilibrium $E_3$, rendering the former infeasible. This transition is accompanied by a noteworthy observation that as $\eta$ increases toward the $BP$ point, the density of top predators steadily rises, whereas that of intermediate predators diminishes due to their consumption by the top predators. Consequently, the density of prey species surges. Upon surpassing the threshold value at the $BP$ point, the intermediate predators face extinction, and the density of prey species saturates at $r/d$. Interestingly, the population of top predators continues to surge, even in the absence of intermediate predators, owing to the availability of alternative food resources. It is important to note that the saturated prey density, coupled with the heightened abundance of top predators, may culminate in the extinction of intermediate predators within the ecosystem.

To investigate the influence of $\eta_0$ on the system's dynamics, we examine the equilibrium curves concerning $\eta_0$, as depicted in Figure 7. In contrast to the effect observed with variations in $\eta$, we note that for lower values of $\eta_0$, the intermediate predator-free equilibrium $E_3$ exhibits stability. At $\eta_0 = 0.40985$, a transcritical bifurcation denoted as $BP$ occurs, marking the transition to a stable branch of the co-existence equilibrium. Additionally, at $\eta_0 = 2.44551$, a supercritical Hopf bifurcation takes place, leading to system instability. It is notable that the impact of $\eta_0$ on the equilibrium densities of species is opposite to the effects witnessed with alterations in $\eta$. Ecologically interpreting the dynamics portrayed in Figures 5 and 7, it becomes evident that system stability is sustained when supplementary food sources, aside from the primary prey (intermediate predators), significantly contribute to the proliferation of top predators. These findings underscore the intricate interplay between different food sources and their effects on the stability and dynamics within the ecosystem.

Thereafter, we present the Hopf curve in the $\eta$-$\eta_0$ plane (depicted in Figure 8(a)), which delineates the bi-parametric regions of stability and instability. The blue-colored region signifies stability, while the brown-colored region denotes instability. This graphical representation succinctly encapsulates the previously discussed dynamics: an increase in $\eta$ tends to stabilize the system, whereas an increase in $\eta_0$ tends to destabilize it. In addition, we have plotted the Hopf curve in the $\eta$-$\delta_3$ bi-parametric plane (Figure 8(b)). In this representation, the colored regions within the parametric plane correspond to feasibility, i.e., $\eta > \delta_3$. Notably, we observe that if the mortality rate of the top predator remains relatively low, the system tends to maintain stability, and the dynamics exhibit no significant alterations with increasing $\eta$. However, for higher values of $\delta_3$, the stability behavior around the co-existence equilibrium shifts from unstable to stable with increasing $\eta$. The critical point $GH$ delineates the generalized Hopf bifurcation point, signifying a transition in the nature of the Hopf bifurcation from supercritical to subcritical, or vice-versa. Biologically interpreting these findings, it implies that for top predator species with shorter life spans, the growth rate resulting from additional food resources, besides their natural prey (intermediate predator), must be sufficiently substantial to maintain system stability. Conversely, for species with longer life spans, no such stringent condition is necessary to ensure stability. These insights shed light on the complex interplay between top predator characteristics, food resource availability, and the resulting ecosystem stability.

4.3.   Impacts of fear factors on the dynamics of system (2.2)

To explore the influence of fear parameters on the dynamics of system (2.2), we first plot equilibrium curves concerning $\omega_1$, while assigning other fear parameters zero. In Figure 9, equilibrium curves are depicted for three distinct values of $\alpha_1$, denoted as $(1)$, $(2)$, and $(3)$ corresponding to $\alpha_1 = 8$, $\alpha_1 = 10$, and $\alpha_1 = 12$, respectively. The blue and red color signifies stable and unstable equilibria of the system. In curve $(1)$, it is evident that the system maintains stability without any alteration as $\omega_1$ varies. However, curve $(2)$ demonstrates a different behavior, showcasing two instances of Hopf bifurcation occurring at $\omega_1 = 0.24314$ and $\omega_1 = 0.47583$. Within the interval $0.24314 < \omega_1 < 0.47583$, the system exhibits oscillatory behavior around the interior equilibrium, while outside this range, it remains dynamically stable. The corresponding Hopf bifurcation diagram is presented in Figure 10(a), showcasing a phenomenon referred to as a "bubbling" effect induced by the fear parameter $\omega_1$. This observation aligns with findings reported by ^[29,49]. For curve $(3)$, a distinctive pattern emerges where the system displays oscillatory behavior at lower values of $\omega_1$, transitioning into stability through a supercritical Hopf bifurcation at larger values (occurs at $\omega_1 = 1.487185$). The associated Hopf bifurcation diagram is illustrated in Figure 10(b), revealing a declining trend in the equilibrium density of prey species with increasing values of $\omega_1$. Figure 11 delineates the region of stability and instability, demarcated by the Hopf curve, in the $\omega_1-\alpha_1$ parametric plane. Notably, in scenarios involving intermediate predators with higher predation rates, the prey species exhibit heightened anti-predation activities, perceiving greater predation risk. Consequently, all species coexist stably, albeit with lower prey density.

Subsequently, Hopf bifurcation diagrams were generated concerning the fear parameters $\omega_2$, $\omega_3$, and $\omega_4$ in Figure 12(a)–(c), respectively. These illustrations reveal a consistent trend where all fear parameters exhibit a stabilizing influence on the system's dynamics. To further analyze the stability and instability regions in the $\omega_1-\omega_2$, $\omega_1-\omega_4$, and $\omega_2-\omega_4$ parameter spaces, Figure 13(a)–(c) were constructed. Specifically, in Figure 13(a), Hopf curves are plotted for four distinct values of $\alpha_1$. It is observed that an increase in the predation rate of the intermediate predator ($\alpha_1$) constricts the stability region. Notably, for $\alpha_1 = 16$, only region Ⅰ remains stable, while for $\alpha_1 = 14$, stability expands to encompass regions Ⅰ $\bigcup$ Ⅱ. Following this trend, for $\alpha_1 = 11$, stability encompasses regions Ⅰ $\bigcup$ Ⅱ $\bigcup$ Ⅲ $\bigcup$ Ⅳ, leaving region Ⅴ as the sole region of instability. Thus, Figure 13(a) illustrates that when either the cost of intermediate predator-induced fear on growth rate or intra-species competition of prey species is high, the system exhibits stable behavior around the interior equilibrium. Conversely, when both costs are relatively low, the system becomes unstable. This implies that an excessive predation rate of the intermediate predator imposes additional costs on prey, leading to increased investment in vigilance and reduced foraging areas, thereby escalating intra-species competition. Moreover, Figure 13(b) illustrates that achieving stability depends on the interplay between the cost of intermediate predator-induced fear on growth rate and the cost of top predator-induced fear on the intra-species competition of the intermediate predator. Interestingly, Figure 13(c) demonstrates a straightforward influence of the costs associated with fear-induced intra-species competition. Higher values of either $\omega_2$ or $\omega_4$ lead to system stability, whereas lower values of both parameters result in system instability.

The alteration in the parameter $\omega_3$ yields distinct population density variations among the three species, as illustrated by the bar diagrams in Figure 14. The graphical representation demonstrates a positive correlation between the increment in the value of $\omega_3$ and the equilibrium density of the prey species, while concurrently observing a decrease in the equilibrium densities of intermediate and top predators. This ecological trend suggests that heightened apprehension induced by top predators triggers behavioral adjustments in intermediate predators, compelling them to limit their foraging activities, consequently leading to a decline in their population density. The reduced foraging efforts of intermediate predators, in turn, alleviate predation pressure on the prey species, culminating in an eventual surge in their population density. The decline in the density of top predators is primarily linked to the scarcity of intermediate predators, as these constitute a favored food source for the top predators. This qualitative behavioral pattern is in accordance with empirical findings observed by Suraci et al. ^[20], wherein the fear elicited by large carnivorous animals initiates cascading effects on lower trophic levels. In summary, alterations in $\omega_3$ exhibit a pronounced impact on predator-prey dynamics, eliciting shifts in population densities among trophic levels. The observed changes underscore the intricate interplay of fear-driven behavioral adaptations and predator-prey interactions, aligning with empirical evidence highlighting the cascading effects of fear-induced responses in ecological systems.

5.   Conclusions

Various factors, including resource availability, body size, and ecological niches, play pivotal roles in shaping the prey preferences exhibited by predator species. While some predators exhibit highly selective feeding behaviors, others display a more generalized approach to their dietary habits. Theoretical investigations into predator-prey dynamics and food web structures have predominantly centered on specialized predators. However, it's crucial to recognize the significance of generalist predators within ecosystems, as they contribute substantially to biodiversity maintenance and serve as effective bio-control agents. In this study, we have constructed a food chain model comprising three distinct species: prey, intermediate predators (of specialist nature), and top predators (generalists). The growth dynamics of the top predators, when devoid of their primary prey (intermediate predators), are described using a Beverton-Holt like function. Beyond the direct impact of predation on species demography, we have also accounted for the non-consumptive effects of predation. Specifically, we incorporated the "cost of fear", introducing heightened intraspecies competition among both the primary prey and intermediate predators, leading to a reduction in their reproduction rates. Empirical evidence supports the notion that the presence of top predators suppresses foraging and predation activities among intermediate predators. Hence, we have adjusted the predation rate of intermediate predators by introducing a decreasing function, which is influenced by the fear parameter and the density of top predators. Our analysis extensively delves into the qualitative behaviors exhibited by the proposed model. We have scrutinized the potential existence of Hopf-bifurcation, investigating the direction and stability of the resultant periodic solution arising from this bifurcation phenomenon. This study contributes to the deeper understanding of predator-prey dynamics, shedding light on the intricate interplay between species interactions, predator behaviors, and the broader ecological implications within food chain systems.

Our findings highlight the stability of the intermediate predator-free equilibrium under conditions where the maximum growth rate falls below a critical threshold. In such instances, the intermediate predator faces extinction within the system. This phenomenon can be elucidated by the significant contribution of supplementary food sources to the exponential growth of top predators, thereby elevating their population density, consequently leading to a decline in the density of intermediate predators. Although this benefits the prey species initially, their density reaches saturation due to limited resources. As a result, a persistent increase in top predator density, coupled with constrained prey availability, ultimately drives the intermediate predators to extinction. However, in scenarios where there is minimal fear of predation and higher natural mortality rates among top predators (resulting in shorter life spans), the role of additional food sources in fostering the growth of top predators must be substantially significant for the ecosystem's dynamical stability. Our research also illustrates the influence of fear parameters on the equilibrium densities of prey and predators, as well as on the dynamics within the considered food chain system. Notably, all fear parameters exhibit stabilizing effects, robustly contributing to the system's stability. Heightened predation pressure triggers increased anti-predation activities among prey species, which are imperative for the stable coexistence of all species within the ecosystem. Fear induced by intermediate predators, whether manifesting as a reduced birth rate or heightened intraspecies competition, diminishes the abundance of prey species. Conversely, the fear induced by top predators on intermediate predators enhances prey density. This fear suppresses the foraging and predation activities of intermediate predators, setting off a trophic cascading effect within lower trophic levels. Reduced foraging activities consequently diminish intermediate predator abundance, ultimately amplifying the abundance of prey species. Consequently, the deliberate manipulation or propagation of fear can emerge as a valuable tool for augmenting endangered species and ensuring a balanced ecosystem. Additionally, the introduction of top predators could prove advantageous for biodiversity conservation efforts.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors wish to extend their heartfelt appreciation to the reviewers for their invaluable comments and insightful suggestions, which have significantly contributed in enhancing the quality of this paper. Soumitra Pal is thankful to the Council of Scientific and Industrial Research (CSIR), Government of India for providing financial support in the form of a senior research fellowship (No. 09/013(0915)/2019-EMR-I). The research of Hao Wang is partially supported by the Natural Sciences and Engineering Research Council of Canada (Discovery Grant RGPIN-2020-03911 and Accelerator Award RGPAS-2020-00090).

Data availability

The data is provided within the article that supports the results of the study.

Conflict of interest

There are no conflicts of interest disclosed by the authors. Hao Wang is an editor-in-chief for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article.