Dynamical system of quokka population depicting Fennecaphobia by Vulpes vulpes

1.   Introduction

Modeling the relationship between prey and predators among species plays an important role in mathematical biology. In this paper, the interaction between the quokka (Setonix brachyrus) and red fox (Vulpes vulpes) populations is deciphered. Quokkas are small marsupials (herbivores) that are endemic to the mainland of South-Western Australia and two islands, and are mainly distributed in coastal and high-rainfall areas ^[1]. The decline in the quokka population was due to agricultural development and the expansion of housing. Due to these factors, they shifted to the island and forests. Their roaming is restricted to a very small range with a small number of scattered populations ^[2]. Quokkas store fat in their tails and use it when food is scarce ^[3]. The abundance of quokkas decreased due to red fox predation and habitat alteration ^[4]. This led to the collapse of the metapopulation structure in a state of non-equilibrium such that the conservation status of the quokka is listed as vulnerable according to the International Union for Conservation of Nature ^[5]. The quokka population is extensively studied on Rottnest Island, where the ecology is different from the mainland population. Quokkas experience seasonal mortality in the summer that is attributable to protein starvation from feeding on nutrient deficient succulents to obtain enough moisture to survive ^[5]. The Quokka population should be well managed in order to protect this species.

Mathematical models provide a detailed understanding of the conservation and management of different species. Roy et al. ^[6] proposed a tritrophic food chain model to discuss the extinction dynamics of the quokka population and its reduction by red foxes. Belew et al. ^[7] analyzed a predator-prey system with hunting cooperation between predators, where the prey population is linearly harvested and is affected by fear. The authors ^[7] concluded that the fear effect, hunting cooperation, prey harvesting, and predator harvesting effectively changed the dynamics of the system. Arancibia et al. ^[8] proposed a predator-prey model with a Holling type Ⅱ functional response, the Allee effect, and a generalist predator, and concluded about the weak and strong Allee effect on the system dynamics. Recently, the complex dynamics of the discrete prey-predator model along with the Allee effect ^[9], fear effect ^[10], and refuge effect ^[11] has been demonstrated with the help of stability and bifurcation analyses.

The study of fear dynamics in a predator-prey model formulation has received increasing interest in recent years. The "ecology of fear" refers to the total impact of predators on prey communities ^[12]. A prey-predator system that introduces the cost of fear into prey reproduction with a Holling type-Ⅱ functional response was analyzed, and it was concluded that strong antipredator responses can stabilize the prey-predator interaction by neglecting the existence of periodic behaviors ^[13]. Tian and Li ^[14] proposed a prey-predator fishery system that introduced the cost of fear into prey reproduction with a Holling type Ⅱ functional response and prey-dependent harvesting. The authors ^[14] concluded that the level of fear altered the stability of the positive equilibrium in the absence of harvesting. Various continuous predator-prey models with the Allee effect, fear effect, and prey refuge have been proposed ^[15,16], and it was concluded that as the refuge of the prey or Allee effect increased, the dynamics of the system exhibited a complex behavior ^[17].

Several researchers proposed and investigated the fear effect of the predator-prey model with cross-diffusion. Zhang et al. ^[18] analyzed a delayed diffusive predator-prey model that incorporated spatial memory and a nonlocal fear effect. Sarkar and Khajanchi ^[19] discovered that the fear effect in a diffusive predator-prey system, which includes the mutual interference between predators, can result in more complex dynamics. Sun ^[20] found that the nonlocal fear effect can alter the stability of a constant steady state, ultimately leading to spatially nonhomogeneous steady states. Moreover, high levels of fear can stabilize the model by preventing the existence of periodic solutions. An increase in fear is associated with a decrease in the growth rate of prey as a result of predator anxiety ^[22]. Prey can relocate their grazing location to a safer area due to the fear of predators, thereby sacrificing the areas with the highest intake rates. Therefore, exploring the spatial dynamics of the prey-predator model is one of the most important factors. A spatio-temporal model of the predator-prey interaction integrates both spatial and temporal systems to investigate how the prey-predator populations vary over time and space ^[23,24]. Cross-diffusion plays a crucial role in the formation of patterns. The concept of cross-diffusion was first introduced by Shigesada et al. ^[25] in their study of competitive species systems, where the population pressure of one species is exerted on another. Since then, many researchers have employed this technique to investigate the Turing instability driven by cross-diffusion and demonstrated the existence of nonhomogeneous steady states that result in spatial patterns ^[26,27]. For instance, a study revealed a reverse correlation between the prey and predator populations due to the effect of cross-diffusion ^[26]. Turing patterns are generated from a diffusion-driven instability, known as the Turing instability. These patterns represent stationary solutions of the reaction-diffusion systems. Turing instability occurs when the stable, homogeneous steady state of the system becomes unstable due to a small heterogeneous perturbation around this steady state ^[27]. As a result, it produces stationary patterns such as spots, stripes, or a combination of spots and stripes ^[28,29].

Our work is motivated by the above research work ^[8,12] that can be used to describe the dynamics of the quokka population. Mainly, we focused on the fear effect and cross-diffusion along with a Holling type Ⅱ functional response and a generalist predator for the proposed system. The predator population control is essential for the survival of the quokka population, so that the fear effect due to predation can be reduced. The paper is organized as follows. A predator-prey model with the fear effect is formulated and converted into a topologically equivalent model in Section 2. In Section 3, the existence of all equilibrium points and its stability analysis is discussed. The appearance of different bifurcations that cause instability under certain conditions is demonstrated in Section 4. Cross-diffusion is incorporated into the equivalent topological model system, and the instability criterion for the Turing pattern is computed in Section 5. In Section 6, the numerical analysis and discussions are carried out to validate the theoretical results. Finally, the paper is concluded in Section 7.

2.   Formulation of predator-prey model with fear effect

In this section, we will develop a modified Leslie-Gower prey predator model by taking the Holling type Ⅱ functional response and the logistic growth rate with fear effect into account. The model formulation is based on the following assumptions:

(ⅰ) The quokka population increases with the intrinsic growth rate $r_{1}Q \left(1-\frac{Q}{K_{1}} \right)$, where the growth rate is represented by $r_{1}$ and the environmental carrying capacity is represented by $K_{1}$.

(ⅱ) Quokka sacrifice themselves in order to save their younger ones from predators, which is considered as the fear effect. When the fear of predation by red foxes increases, the growth rate of the quokka population decreases; therefore, a nonlinear term for fear is assumed in the form ${\frac{1}{1+f F}}$, where $f$ is the fear effect due to predation.

(ⅲ) When considering red foxes as generalist predators and quokka species as their favorite food, in the predator-prey model, the Holling type Ⅱ functional response is used in the form ${\frac{b Q F}{ Q+a }}$, where $b$ is the search rate for prey of $F$ and ${\frac{Q}{Q+a}}$ is the Holling type Ⅱ functional response, where $a$ is the measure of the extent to which the environment provides refuge for the quokka population.

(ⅳ) The red fox population grows with the logistic growth rate $r_{2}$. The Leslie-Gower formulation ${\frac{F}{K_{2}Q+c}}$ is based on the assumption that a reduction in a predator population has a reciprocal relationship with the per capita availability of its preferred food. Since the red fox is a generalist predator, a constant $c$ is added to the prey-dependent carrying capacity to measure the loss of the predator population in the absence of its favorite food.

Now, considering all the above assumptions, the predator-prey model along with fear effect of predation is given by the following:

with $Q\geq 0, F\geq 0$. The above system is defined by $\Omega = \{(Q, F) \in \mathbb{R}^2 \mid Q \geq 0, F \geq 0 \}$. All parameters are non-negative and their descriptions are listed in Table 1.

The above system (2.1) is converted to a topologically equivalent model (2.2) to simplify the analysis ^[8,30,31]. We introduce the change of variable and time re-scaling given by the function $\varphi: \hat{\Omega} \times \mathbb{R} \rightarrow \Omega \times \mathbb{R}$, where $\varphi(u, v, \tau) = (Q, F, t)$ and $\hat{\Omega} = \{(u, v) \in \mathbb{R}^2 \mid u \geq 0, v \geq 0 \}$ so that

The Jacobian matrix of $\varphi$ is as follows:

Then,

Therefore, $\varphi$ is a diffeomorphism for which the system (2.1) is a topologically equivalent system (2.2) and the associated system of differential equations is as follows:

where

3.   Existence and stability analysis of the equilibrium points

The existence of equilibrium points and its stability analysis are discussed in this section. The equilibrium points are essential to carry out a system's stability. A system reaches equilibrium when its state variables remain fixed and unchanged. The following are the equilibria of the system (2.2):

(ⅰ) The trivial equilibrium point $\epsilon_{0}(0, 0)$ exists;

(ⅱ) The predator-free equilibrium point $\epsilon_{1}(1, 0)$ exists;

(ⅲ) The prey-free equilibrium point $\epsilon_{2}(0, C)$ exists; and

(ⅳ) There are two endemic equilibrium points $\epsilon^{*}(u^{*}, v^{*})$ co-exists and is given by

where

Two endemic equilibrium points, $\epsilon^*(u_1^*, v_1^*)$ and $\epsilon^*(u_2^*, v_2^*)$, exist if the following conditions (3.1) and (3.2), hold:

and

We need to compute a Jacobian matrix $J$ that corresponds to each equilibrium point of the system (2.2) to determine its stability.

Theorem 1. In the system (2.2), the equilibrium point $\epsilon_0 (0, 0)$ is unstable, while the equilibrium point $\epsilon_1 (1, 0)$ is a saddle point and therefore also unstable. The equilibrium point $\epsilon_2 (0, C)$ is stable provided that the condition $A < B C \left(\frac{C}{D} + 1 \right)$ holds. Furthermore, there are two endemic equilibrium points, one is stable $\epsilon^* (u_{1}^{*}, v_{1}^{*})$ while the other one $\epsilon^*(u_{2}^{*}, v_{2}^{*})$, is unstable.

The variational matrix for the system (2.2) at $\epsilon_{0}(0, 0)$ is given by the following:

Here, both the eigenvalues $ACDE$ and $ACDR$ are clearly positive at the equilibrium point $\epsilon_{0}$. Therefore, the system dynamics (2.2) is unstable at $\epsilon_{0}(0, 0)$.

The variational matrix for the system (2.2) at $\epsilon_{1}(1, 0)$ is given by the following:

Here, one eigenvalue $(1+A)(1+C)D R$ is positive and the other one $-(1+A)(1+C)DE$ is clearly negative. Therefore, the equilibrium point $\epsilon_{1}(1, 0)$ is a saddle point and the system dynamics is unstable.

The variational matrix for the system (2.2) at the point $\epsilon_{2}(0, C)$ is given by the following:

Here, one eigenvalue is clearly negative $-AC(C+D)R$; however, the other eigenvalue $-C(BC(C+D)-AD)E$ is negative if $A < B C \left(\frac{C}{D}+1 \right)$ holds. Thus, the system dynamics is stable if $A < B C \left(\frac{C}{D}+1 \right)$ hold; otherwise, this is unstable at $\epsilon_{2}(0, C)$.

The variational matrix for the system (2.2) at the point $\epsilon^{*}(u^*, v^*)$ is given by the following:

The dynamics of the system will be stable if the above $\det(J(u^*, v^*))$ and $\text{tr}(J(u^*, v^*))$ are, respectively, positive and negative at $\epsilon^*$. Therefore, the system is stable at the endemic equilibrium point if $u^* \ge \frac{1}{2}$ and $v^* \le A + 2 u^*$ satisfy; otherwise, the system will be unstable.

4.   Bifurcation analysis

A bifurcation is defined as a sudden qualitative change in the nature of the dynamics that occurs for the most sensitive parameter.

4.1.   Existence of Hopf bifurcation

A Hopf bifurcation is a point where the equilibrium changes its stability through a pair of purely imaginary eigenvalues, giving rise to a limit cycle from an equilibrium.

Theorem 2. The necessary and sufficient conditions for the occurrence of a Hopf-bifurcation at the point $\epsilon^{*}(u^{*}, v^{*})$ when the critical value $B = B_{c} = {\frac{R(1-u^{*})(u^{*}+A)^2}{Eu^{*}(u^{*}+C)(1-2u^{*}-A)}}$ are

(a) ${\operatorname{tr}(J) \big|_{\epsilon^*} = 0}$,

(b) ${\det(J) \big|_{\epsilon^*} > 0}$,

(c) ${\frac{d}{dB} \operatorname{tr}(J) \big|_{\epsilon^*} \neq 0}$.

Proof. The characteristic equation of (3.3) at the endemic point is given by the following:

The trace in the above characteristic equation at $\epsilon^{*}(u^{*}, v^{*})$ is zero when $B = B_{c}$. Since $\text{tr}(J(u^{*}, v^{*})) = 0$, the characteristic equation at $B_c$ reduces to the following:

where $\det(J) = E R u^* (u^*+A) (u^*+C) v^* (v^*+D) \left(D (2u^*-1 + A) + \frac{R (u^*-1) (u^*+A) (2u^*+A+C-v^*) (2v^*+D)} {u^* (1+E (u^*-1) (u^*+C))+A} \right),$ which is greater than $0$ at $\epsilon^*(u^*, v^*)$. Next, to confirm the occurrence of a Hopf bifurcation, we must check the following transversality condition at $B = B_c$:

Thus, the occurrence of a Hopf bifurcation is confirmed at the critical point $B_c$. □

Example 4.1. The occurrence of a Hopf bifurcation is numerically confirmed around the point $\epsilon^*(0.390721, 0.640721)$ at $B_c = 0.071387$ for a set of parameter values considered in Eq (6.1). Here, all the conditions for the existence of the Hopf bifurcation $\operatorname{tr}(J) \big|_{\epsilon^*} = 0$, $\det(J) \big|_{\epsilon^*} = 0.0000352064 > 0$, and $\frac{d}{dB} \operatorname{tr}(J) \big|_{\epsilon^*} = 0.671263 \neq 0$ is satisfied at this critical $B_c$.

4.2.   Existence of saddle-node bifurcation

For system (2.2), the $v$ nullclines are $v = 0$ and $v = u+C$, and the $u$ nullclines are $u = 0$ and $v = \frac{-BD \pm \sqrt{(B D)^2 + 4B D (1 - u) (u+A)}}{2B}$. The endemic equilibrium point $u^*$ is the solution of $\frac{-BD \pm \sqrt{(B D)^2 + 4B D (1 - u) (u+A)}}{2B} = u+C$, which is equivalently as follows:

The roots of (4.1) are given by the following:

A saddle-node bifurcation or fold bifurcation is a local bifurcation in which two fixed points of a dynamical system collide and annihilate each other. Therefore, when $\Delta = 0$, the system has only one positive equilibrium point, that is, $u_{1} = u_{2} = G = \frac{-2BC+D(1-A-B)}{2(B+D)}$. The corresponding Jacobian matrix will be $J(G, G+C)$:

Theorem 3. The system (2.2) undergoes a saddle-node bifurcation at the equilibrium point $(G, G+C)$, where $B_{sn} = {\frac{D(A+2G-1)}{2C+D+2G}}$.

Proof. We will prove that the system (2.2) undergoes a saddle-node bifurcation at the critical point $B = B_{sn}$ based on Sotomayor's theorem ^[32]. If $\Delta = 0$, then the system possesses a single positive equilibrium point, denoted as $(G, G+C)$. The determinant of the Jacobian matrix at this equilibrium point is given by the following:

For the occurrence of a saddle-node bifurcation, the determinant of the Jacobian matrix $\det(J(G, G+C))$ must be zero. Hence, the saddle-node bifurcation exhibits at the critical point $B_{sn} = \frac{D(A+2G-1)}{2C+D+2G}$.

The Jacobian matrix (4.2) at the critical point $B_{sn}$ is given by the following:

The vector form of system (2.2) is given by the following:

Let

be the eigenvectors corresponding to the eigenvalue $\Delta = 0$ of $J(G, G+C)$ and $J(G, G+C)^\top$, respectively.

Now, we differentiate the vector function (4.4) with respect to the bifurcation parameter $B$; then, we obtain the following:

Therefore,

Now, we compute $W^\top[D^{2}f(u, v;B)(V, V)]$, where $V = (v_{1}, v_{2})$ and $D^{2}f(u, v;B)(V, V)$ is given by the following:

Hence, the system (2.2) has a saddle-node bifurcation as $0 < A < 1$ at $(G, G+C)$ by Sotomayor's theorem ^[32]. □

This theorem implies that when the predator search rate ($B$) varies, it significantly affects the predation pressure on the prey. If search efficiency exceeds a critical threshold, then the prey mortality dramatically increases, thus compromising their survival. Initially, one endemic equilibrium point is stable where the prey and predator co-exist, and the other endemic point is unstable. As the search rate increases, these equilibrium points merge at a critical point, thus causing a saddle-node bifurcation. Once the threshold is crossed, stable co-existence is no longer possible, thus leading to the collapse of the prey population.

4.3.   Existence of Bogdanov-Takens bifurcation

A double zero eigenvalue of the Jacobian matrix is obtained when $\text{tr}(G, G+C) = 0$, even if it is not a zero matrix. Therefore, there is the possibility of a two co-dimension bifurcation occurring. Let $A$ and $B$ be the bifurcation parameters. Now, we will prove the non-degeneracy condition of the Bogdanov-Takens ($BT$) bifurcation. Let $(A, B) = (A^{[BT]}+\lambda_{1}, B^{[BT]}+\lambda_{2})$ be the neighboring point of the $BT$ point $(A^{[BT]}, B^{[BT]})$ for a sufficiently small $\lambda_{i}$ with $i = 1, 2$. Then, the system (2.2) at $(A, B) = (A^{[BT]} + \lambda_{1}, B^{[BT]} + \lambda_{2})$ is given by the following:

The system (4.5) is $C^{\infty}$ smooth with respect to the variables $u, v$ in small neighborhood of $(A^{[BT]}, B^{[BT]})$.

We define $u_{1} = u-u^{*}$ and $u_{2} = v-v^{*}$; then, the above system (4.5) reduces to the following:

where

and $P_{1}, P_{2}$ are the power series in $(u_{1}, u_{2})$ with powers $u^{i}_{1} u^{j}_{2}$ satisfying $i+j\geq3$. The non-degeneracy conditions of the $BT$ bifurcation ^[33] are as follows:

Thus, we have the following theorem.

Theorem 4. Let the system parameter be such that $D = -\frac{(C + G) (A R + G R)}{ A E G-E G + 2 E G^2 + A R + G R}$, where $G = \frac{-2BC+D(1-A-B)}{2(B+D)}$, $\Delta = 0$; then, the system (2.2) experiences a $BT$ bifurcation.

Proof. If the system (2.2) exhibits a BT bifurcation at $D = D_{BT}$, then the eigenvalues corresponding to the bifurcation points will be zero. The eigenvalues are $\lambda_1 = 0$ and $\lambda_2 = (C + G) (D E G - A D E G - 2 D E G^2 - A C R - A D R - A G R - C G R - D G R - G^2 R)$. At $D_{BT} = -\frac{(C + G) (A R + G R)}{ A E G-E G + 2 E G^2 + A R + G R}$, $\lambda_2$ becomes zero. Therefore, the Jacobian matrix at the equilibrium point $(G, G+C)$ is simplified to the following:

Now, we find the Jordan normal form of $J(G, G+C)$ that has equal eigenvalues and a unique eigenvector $\begin{pmatrix} 1 \\ 1 \end{pmatrix}$. This vector will form the first column of the transformation matrix $\Gamma$. To obtain the second column of $\Gamma$, we choose a vector $\begin{pmatrix} 1 \\ 0 \end{pmatrix}$.

Thus, $\Gamma = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}$ and $\Gamma^{-1}(J(G, G+C))\Gamma = \begin{pmatrix} 0 & \frac{E G (A + G) (C + G)^2 (-1 + A + 2G) R} {E G (-1 + A + 2G) + (A + G) R} \\ 0 & 0 \end{pmatrix}$.

Hence, we have a $BT$ bifurcation of co-dimension two ^[34]. □

5.   Stability analysis of spatio-temporal system

In this section, diffusion is incorporated into the proposed system (2.2) to decipher the spatial dynamics and the effect of cross-diffusion on the system dynamics. A spatio-temporal model is introduced in the following system to demonstrate the spatial interactions between predator and prey:

with the initial conditions

and the no-flux boundary conditions

where $u(x, y, t)$ and $v(x, y, t)$ denote the population densities of the quokka and red fox, respectively, at position $(x, y) \in \Omega$ and time $t$. Here, $d_1$ and $d_2$ represent self-diffusion for the quokka and red fox, respectively. The Laplace operator $\nabla^2$ is defined as $\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$, which describes diffusion in two dimensions. Furthermore, $d_{12}$ is the cross-diffusion coefficient which models the movement of the quokka species depending on the red fox population. Similarly, $d_{21}$ is the cross-diffusion coefficient for the red fox, which models their movement as being influenced by the presence of the quokka populations. $\Omega$ is the spatial domain and $\eta$ is the unit outward normal vector to $\partial\Omega$.

To discuss the linear stability analysis of the system (5.1) about the homogeneous steady state $\epsilon^{*} (u^{*}, v^{*})$, we linearize the system (5.1) by the following perturbations ^[35]:

where $\lambda$ is the growth rate of perturbations and $0 < \eta_1, \eta_2 \ll 1$. Additionally, ${\bf{{k}}} = ({k}_x, {k}_y)$ is the wave number, and its value is $k = |\bf{{k}}|$.

After linearizing the system (5.1) using the perturbation (5.2), the characteristic equation for the growth rate $\lambda$ is derived as follows:

where

is the diffusion matrix, ${\mathbb{I}_{2 \times 2}}$ is a $2 \times 2$ identity matrix, and $J$ is the Jacobian matrix at $\epsilon^{*}(u^{*}, v^{*})$, which is the same as given in (3.3).

The characteristic equation (5.3) can be written in the following simplified form:

The characteristic equation (5.4) has the following solutions:

Turing instability occurs when $\Re \left[\lambda(k^2)\right] > 0$ for some $k \neq 0$. Since $g(k^2) > 0, \, \, \forall k\neq 0$, instability may arise only from $h(k^2) < 0$ for some $k\neq 0$. Given that $h(k^2)$ is a quadratic polynomial in $k^2$, the minimal value of $h(k^2)$ can be obtained at the critical value $k^2 = k^2_{T}$, where

We require $\det({\bf{D}}): = d_1 d_2 - d_{12}d_{21} > 0$ for the well-posedness of the spatio-temporal system. Consequently, the Turing threshold is obtained if the following conditions hold:

$k_{T}$ is a real number. Additionally, $h(k^2_T) < 0$ provided

Theorem 5. The condition of diffusive instability for system (5.1) is given by the conditions $d_1d_2-d_{12}d_{21} > 0$, (5.5), and (5.6), provided that $\epsilon^*(x^*, y^*)$ is linearly stable in the absence of diffusion.

6.   Numerical simulations and discussions

In this section, the model systems (2.2) and (5.1) are numerically simulated in detail to investigate the dynamics of the system. The temporal system is integrated using the ode45 functions with the help of MATLAB and the MatCont software ^[37]. The spatiotemporal system (5.1) is numerically solved using the finite difference scheme. The forward difference and central difference schemes are employed for the reaction and diffusion terms, respectively, with the help of MATLAB. $\Delta h$ and $\Delta t$ correspond to the distance between consecutive spatial points and the discrete step interval of time, respectively. The spatial domain is considered of the size $200 \times 200$, the space step is $\Delta h = 0.1$, and the time step is $\Delta t = 0.001$. The following set of parameter values is fixed throughout the numerical simulation:

Nullclines of the system (2.2) are illustrated in Figure 1 using the parameter values provided in (6.1). The blue and red curves represent the nullclines for the prey(quokka) and predator(red fox), respectively. In this figure, the direction field is also plotted to understand how the curves move around in the plane so that one can identify the stability and instability of the equilibrium points. The intersections of these curves indicate the equilibrium points (represented by red dots), which include $\epsilon_{0} (0, 0)$, $\epsilon_{1}(1, 0)$, and $\epsilon_{2}(0, 0.25)$, and two endemic equilibriums, $\epsilon^*_u (0.0390094, 0.289009)$ and $\epsilon^*_s (0.544741, 0.794741)$. Here $\epsilon_{0}$, $\epsilon_{1}$, and $\epsilon^*_{u}$ are unstable, while $\epsilon_{2}$ and $\epsilon^*_{s}$ are stable under certain criteria. The determinant and trace at the endemic point $\epsilon^*(0.0390094, 0.289009)$ are $\det(J(u^{*}_{1}, v^{*}_{1})) = -6.3674 \times 10^{-08}$ (negative) and $\text{tr}(J(u^{*}_{1}, v^{*}_{1})) = 0.000663683$ (positive). Thus, the system (2.2) is unstable according to the Routh-Hurwitz criteria at this point $\epsilon^*(0.0390094, 0.289009)$. The determinant and trace at the other endemic point $\epsilon^*(0.544741, 0.794741)$ are $\det(J(u^{*}_{2}, v^{*}_{2})) = 0.000180743$ (positive) and $\text{tr}(J(u^{*}_{2}, v^{*}_{2})) = -0.0113229$ (negative). Therefore, the system (2.2) is stable at this endemic point $\epsilon^*(0.544741, 0.794741)$ according to the Routh-Hurwitz criteria.

6.1.   Uncertainty and sensitivity analysis

Here, we calculate the global sensitivity index of the endemic equilibrium point to the parameter in the model. This index shows how targeting specific parameters has an impact on the global achievement of interest. Since the relationship between the parameters and output is not pre-known, using the partial rank correlation coefficient (PRCC) method is ideal. Hence, the PRCC method is used to find which parameter has a high impact on the endemic equilibrium point, so that it should be targeted for effective and efficient intervention. We begin the uncertainty and sensitivity analyses by selecting the value and range of the parameters. If the parameter estimation is uncertain, then we can treat each parameter as a random variable with a corresponding probability density function, and a uniform probability distribution function is assigned to each parameter ^[36]. The range for all parameters is considered to be $[0, 1]$. The PRCC sensitivity results for the endemic equilibrium point $\epsilon^{*}$ to all the different parameters are illustrated in Figure 2. The PRCC sign indicates the direction of association between the input and the output factors. A value $+1$ indicates that there is a positive linear relationship, $-1$ indicates that there is a negative linear relationship, and a value near $0$ indicates that there is no relationship. The PRCC sensitivity index of the system (2.2) is given in Table 2, and the most negatively influential parameter at the endemic equilibrium point $\epsilon^{*}$ is $B$ followed by $D$. In addition, the most positively influenced parameter is $A$.

The parameter $B$ indicates that the increased predation pressure leads to a decrease in the prey population. The variable $B$ is inversely related to the intrinsic growth rate of the prey ($r_1$); therefore, higher growth rates lessen the negative effects of predation. Although a higher $D$ (lower fear) allows the prey to freely forage, it can also increase the predator encounters. Finally, the parameter $E$, which is associated with carrying capacity, negatively impacts the population stability by limiting the resources and breeding opportunities, thus leading to fewer prey.

6.2.   Bifurcation numerical discussion

In Figure 3(a), the system experiences a Hopf bifurcation at the endemic point $(u^{*}, v^{*}) = (0.390721, 0.640721)$ with $B_c = 0.071387$ and a saddle-node bifurcation at $B_{sn} = 0.091278$ as the parameter $B$ changes. The Hopf bifurcation indicates that as the predator's searching rate increases, the system transition changes from a stable state to oscillatory predator-prey cycles. This means that the prey-predator population begins to periodically fluctuate, thus resembling real-world predator-prey dynamics. The green line indicates that the system is stable, and the blue dashed line indicates that the system is unstable. The first Lyapunov coefficient is $0.01598375$ at the $H$ point. Its positive value indicates a subcritical Hopf bifurcation. At a higher value of $B = 0.091278$, the limit point (LP) indicates a sudden collapse of one of the equilibrium states. The $LP$ is reported at $(u^{*}, v^{*}) = (0.186340, 0.436340)$ with the eigenvalues $0$ and $0.002816$ and the normal form coefficient $a = -0.002411581$. In Figure 3(b), the system only experiences a Hopf bifurcation at a higher critical point value $(B_c = 0.154163)$ when $C = 0$. The first Lyapunov coefficient is $-0.01083873$ at the $H$ point in Figure 3(b). This means that a stable limit-cycle bifurcates from the equilibrium when it loses stability. Its negative value indicates a supercritical Hopf bifurcation. This suggests that without an alternative food source for predators, a greater predation pressure is necessary for the transition to an oscillatory behavior.

The subcritical Hopf bifurcation in Figure 4(a) and supercritical Hopf bifurcation in Figure 4(b) are plotted at $B_c = 0.071385$ ($C = 0.25$) and $B_c = 0.154164$ ($C = 0$), respectively. In Figure 4, the green and blue lines indicate that the limit cycle is stable and unstable respectively. This indicates that when the predator's searching rate exceeds a specific threshold, stable populations become unstable, thus resulting in periodic fluctuations. This concludes that an increase in the red fox predation pressure can destabilize the quokka populations in a real ecosystem, thus leading to alternating periods of abundance and scarcity.

In Figure 5, the BT point is plotted by varying the parameters $A$ versus $B$. The BT point is reported at the endemic point $(u^{*}, v^{*}) = (0.018609, 0.268609)$, where $(A, B) = (0.121619, 0.183585)$. The stability of the system changes with the occurrence of the $BT$ point.

6.3.   Effect of cross-diffusion

How cross-diffusion and fear effects impact the population distribution and pattern formation has been explored. We have investigated the circumstances under which cross-diffusion causes instability and Turing pattern generation using a linear stability analysis. The set of parameter values is the same as specified in Eq (6.1) with the diffusion coefficients $d_{1} = 0.1, \; d_2 = 0.009$ and $d_{21} = 0.07$ to numerically illustrate the diffusive instability phenomenon. To ensure the validity of the analytical conditions for the Turing instability, the graphs of $h(k^2)$ and $\Re(\lambda(k^2))$ versus $k^2$ are presented in Figure 6 by varying the diffusion coefficient $d_{12}$. Numerically, $d_{12} = 0.005$ is less than the threshold condition $d_{12}^T = 0.0099$, then, the spatio-temporal system is Turing stable. However, if $d_{12} = 0.0125$ is greater than the threshold condition $d_{12}^T$, then the spatio-temporal system is Turing unstable and also satisfies all the criteria of Turing instability as $b_{11} + b_{22} = - 0.0113229 < 0$, $b_{11} b_{22} - b_{12} b_{21} = 0.000180743 > 0$, $d_{1}d_{2}-d_{12}d_{21} = 0.000025 > 0,$ $b_{12}d_{21}+b_{21}d_{12}-b_{11}d_{2}-b_{22}d_{1} = -0.000361555 < 0,$ $(b_{12}d_{21}+b_{21}d_{12}-b_{11}d_{2}-b_{22}d_{1})^2-4(d_{1}d_{2}-d_{12}d_{21})(b_{11}b_{22}-b_{12}b_{21}) = 0.00349948 > 0$.

Figure 7 illustrates how cross-diffusion ($d_{12})$ impacts the movement of the quokka and red foxes in space. Cross-diffusion refers to the phenomenon in which one species moves in response to the presence of the other, rather than moving randomly. The prey and predators form small distinct patches when $d_{12}$ is small. In this scenario, red foxes tend to hunt in specific areas with minimal movement. As $d_{12}$ increases, these patches stretch into more maze-like shapes, thus indicating that predators are increasingly influenced by the distribution of the quokka. When $d_{12}$ is even higher, the patterns evolve into stripes, thus suggesting a more stable and organized distribution of populations. This explores that cross-diffusion is crucial in red fox-quokka interactions, as it contributes to the formation of patterns that may prevent local extinctions and support the long-term survival of both species.

Figure 8 illustrates how increased fear interferes with the quokka aggregation. The prey tend to form isolated patches at a low density. This behavior suggests that strong fear drives them into safe zones, thus limiting their movement to avoid predators. As the density increases, these patches become more elongated and interconnected. This indicates that reduced fear allows the quokka to spread out while still maintaining a structured pattern. At the highest density, the distribution of the prey resembles stripe-like formations. This implies that in the absence of strong fear, the quokka disperse more evenly but still exhibit spatial organization due to interactions with the red fox. This evolution in patterns highlights the important role of fear in shaping the prey movement and distribution. A stronger fear leads to localized clustering, while a weaker fear promotes a broader dispersal, thus balancing the predation risk with the resource availability.

7.   Conclusions

In this paper, the predator-prey interaction model is studied that corresponds to the quantitative measure of the attributes of the quokka and red fox. The quokka population is on the edge of extinction. Apart from natural death (from road kills, climate change, lack of moisture in the environment and food habitat, etc.), the arrival of the red fox in Australia is one of the main reasons for the declination of the quokka population. Quokkas have a fear effect due to the predation by the red fox. The introduction of Europeans into Australia and human activities led to degradation of the quokka habitat. Moreover, this led to an increase in the number of feral foxes that have a negative impact on the quokka population.

We investigated a spatiotemporal predator-prey model with the Leslie-Gower term and Holling type Ⅱ functional response along with a fear effect in the quokka population. In order to simplify the analysis, the temporal system (2.1) is transformed to a topologically equivalent system (2.2). The transformed system (2.2) is analyzed throughout the paper with and without cross-diffusion. The existence and stability analysis of the equilibrium points are performed under certain conditions. The occurrence of a subcritical Hopf bifurcation at $C = 0.25$ and a supercritical Hopf bifurcation at $C = 0$ is reported at the critical point $B_{c}$ (cf. Figure 4). In addition, the system (2.2) exhibits a saddle-node bifurcation at the equilibrium point $(G, G+C)$. Finally, the system (2.2) undergoes a Bogdanov-Takens bifurcation by varying the parameters $A$ and $B$. From the bifurcation analysis, it is concluded that the predation rate is the essential parameter that rules our model system and is useful in understanding the dynamics of the predator-prey model. Sensitivity analysis is performed to find the most influential and least influential parameters using the PRCC method. The effect of fear due to predation plays a major role in the declination of the prey population.

For the spatio-temporal system, the cross-diffusion coefficient $d_{12}$ is a very crucial parameter. The homogeneous steady state is attained when $d_{12} < d_{12}^{T}$. However, the homogeneous state becomes unstable and spatial patterns appear when $d_{12} > d_{12}^{T}$. When the cross-diffusion coefficient exceeded a critical threshold, the system exhibited Turing instability, resulting in spatially heterogeneous population distributions. cross-diffusion, representing the tendency of one species to move in response to the presence of another, is identified as a key factor of pattern formation. Our result indicates that when prey experience a high level of fear due to the presence of predators, they tend to disperse more, leading to fragmented and irregular spatial distributions. As fear reduces, the prey start to form more structured aggregations, transitioning from scattered spots to stripe-like formations. These findings suggest that fear not only impacts the overall population size, but also governs the spatial arrangement of individuals within the ecosystem. Such patterns are ecologically significant as they can indicate regions of high predation pressure.

Use of AI tools declaration

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

Acknowledgments

This work was financially supported by DST INSPIRE FELLOWSHIP, REG. NO. IF230186.

Conflict of interest

The authors declared that there is no conflict of interest.